3. Results and discussion

We conduct simulations of two-dimensional Rayleigh–Bénard convection in porous media over the ranges of Rayleigh number $10^4\leq R\leq 10^{11}$ and Darcy number $5\times 10^{-8}\leq \xi \leq 10^6$. All the results reported here are for porous media with inertia coefficient $c_F=0.145$ and porosity $\epsilon =0.99$ for brevity. A comparison of the Nusselt number versus the Rayleigh number for different Darcy numbers is shown in figure 1(a). The spatio-temporally averaged Nusselt number is calculated by

(3.1)\begin{equation} N=1+(\sigma R)^{1/2}\langle\, u_{\bot}T \,\rangle, \end{equation}

which is obtained from the average in the whole computational domain over 1000 separate snapshots sampled in the final statistically steady states. The convective heat transfer increases with Rayleigh number, and the relationship is strongly modified by the Darcy number. The Nusselt number becomes significantly smaller with decreasing Darcy number; however, the difference diminishes for high Rayleigh numbers. The dependence of the Nusselt number $N$ on the Rayleigh number $R$ can be well represented by $N=0.146\times R^{0.286}$ for $R>2\times 10^7$ at the highest Darcy number $\xi =10^6$, which approaches the relationship for the Rayleigh–Bénard turbulent convection of free fluids (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009). According to (2.4), the Forchheimer term and the Darcy term decrease with increasing Darcy number, and become negligible in comparison with the buoyancy force term when the Darcy number is high enough, though it often falls in the range $\xi \leq 10^{-2}$ for most porous media in previous studies (Guo & Zhao Reference Guo and Zhao2002; Nield & Bejan Reference Nield and Bejan2013). The generalized momentum equation (2.1) reduces to the Navier–Stokes equation for free fluid flows when the porosity is approximately unity, which allows one to study the turbulent convection in both porous media and free fluids by considering a single equation. A fitting line with a slope of 2/7 is provided for reference in figure 1(a), wherein all the data for various Darcy numbers gradually collapse onto the unique line for the highest Rayleigh numbers that approach $R\simeq 10^{11}$. A preasymptotic regime with a slope of unity can be roughly discerned for those data with Darcy number $\xi \leq 2\times 10^{-6}$ and Rayleigh number $R\geq 10^7$, for which a fitting line is included for reference only. The observation is consistent with the results of the Darcy model (Hewitt et al. Reference Hewitt, Neufeld and Lister2012). However, the convective flow departs from the near-Darcian regime with increasing Rayleigh number until it eventually attains the asymptotic regime which exhibits an attenuated slope. The Nusselt number $N$ versus the Rayleigh–Darcy number $R_{\xi }$ is shown in figure 1(b). A collapse of the data occurs when the Darcy number $\xi$ and the Rayleigh–Darcy number $R_{\xi }$ become small. A line with unit slope is included for reference, which indicates the Darcy limit for porous media. A fitting line with a slope of 2/7 can also be obtained for Darcy number $\xi \geq 10^{-2}$, which corresponds to the relationship of the free fluid limit for Rayleigh–Bénard turbulent convection. Note that the data for $\xi =10^6$ have been shifted leftwards by a factor of $10^8$ in figure 1(b) so as to facilitate a comparison with the data for $\xi =10^{-2}$. A gradual transition of the flow state towards the free fluid limit can be observed with increasing Darcy number $\xi$ and Rayleigh–Darcy number $R_{\xi }$.

Figure 1.

Comparison of the Nusselt number $N$ versus the Rayleigh number $R$ (a) and the Rayleigh–Darcy number $R_{\xi }$ (b) for different Darcy numbers $\xi =10^6$ (circle), $10^{-2}$ (cross), $10^{-3}$ (up triangle), $10^{-4}$ (down triangle), $10^{-5}$ (square), $2\times 10^{-6}$ (plus) and $10^{-6}$ (diamond). For comparison, the data for $\xi =10^6$ have been shifted leftwards by a factor of $10^8$ in (b).

A scaling analysis can be conducted based on the generalized Navier–Stokes equations and the heat transport equation in order to identify when the Darcy limit applies. The inertia coefficient $c_F$ and the porosity $\epsilon$ are two constants that can be eliminated, which do not affect the result. The characteristic velocity scale $u$ can be determined by balancing the Darcy term with the buoyancy force term in (2.4) and the pressure gradient term in (2.1), whose scales are of the order of unity:

(3.2)\begin{equation} u \sim \xi\left(\frac{R}{\sigma}\right)^{1/2}. \end{equation}

The convective flow departs from the Darcian regime when the Forchheimer quadratic inertia term becomes important in comparison with the Darcy term:

(3.3)\begin{equation} \xi^{{-}1/2}u \sim \xi^{{-}1}\left(\frac{\sigma}{R}\right)^{1/2}, \end{equation}

where the velocity scale can be approximated by the characteristic velocity. The boundary for Darcy's law to hold is obtained by substituting equation (3.2) into (3.3):

(3.4)\begin{equation} R\,\xi^{3/2} \sim 1, \end{equation}

where the Prandtl number $\sigma$ is assumed to be unity. The Darcy term can also be compared with the Brinkman term, which results in

(3.5)\begin{equation} \xi^{{-}1} \sim l^{{-}2}. \end{equation}

The parameter $l$ is the typical length scale of the flow, which can be estimated by a balance of advection and diffusion in the heat transport equation (2.3):

(3.6)\begin{equation} ul^{{-}1} \sim \frac{1}{(\sigma R)^{1/2}}\,l^{{-}2}. \end{equation}

The typical length scale can be deduced from (3.2) and (3.6):

(3.7)\begin{equation} l \sim \frac{1}{R\,\xi}. \end{equation}

The constraint of (3.4) can be reproduced by substituting equation (3.7) into (3.5). Compare the Darcy term with the convective inertia term

(3.8)\begin{equation} l^{{-}1}u \sim \xi^{{-}1}\left(\frac{\sigma}{R}\right)^{1/2}, \end{equation}

and substitute (3.2) and (3.7) into (3.8), which still leads to the same constraint as in (3.4). The flow is dominated by the Darcian effect in the regime of $R\,\xi ^{3/2} \ll 1$; however, the non-Darcian effects become notable in the regime of $R\,\xi ^{3/2} \gg 1$, where the flow state is dominated by those force terms of the Forchheimer quadratic inertia term, the Brinkman term and the convective inertia term. Darcy's law is more likely to break down for larger $R$, unless the Darcy number $\xi$ is made sufficiently small. The range of the Darcy regime can therefore be roughly located by (3.4), which corresponds to the range $R_{\xi }<10^2$ in figure 1(b).

The Reynolds number $\varLambda$ and the Nusselt number $N$ versus the Darcy number $\xi$ for different Rayleigh numbers are illustrated in figure 2. The spatio-temporally averaged Reynolds number is defined in terms of the velocity fluctuation (Lohse & Toschi Reference Lohse and Toschi2003):

(3.9)\begin{equation} \varLambda =\vartheta h/\nu, \end{equation}

where $\vartheta$ is the root mean square velocity fluctuation. The Reynolds number $\varLambda$ (figure 2a) and the Nusselt number $N$ (figure 2b) both increase with Darcy number above the onset of convection, and a unique critical Darcy number can be roughly identified for any given Rayleigh number $R$ by $\varLambda \geq 0$ and $N\geq 1$, respectively. The critical Darcy numbers are observed to decrease with increasing Rayleigh number in both panels, which indicates that fluid-saturated porous media with a higher Rayleigh number are more prone to buoyancy instability. The phenomenon qualitatively coincides with the result of the Darcy model, in which the critical value of the product of the Darcy number $\xi$ and the Rayleigh number $R$ remains constant at about $4{\rm \pi} ^2$, though the generalized momentum equation is probably not suitable for the analysis near the onset of convection where the Reynolds number is infinitesimally small. We now focus on the asymptotic behaviour of heat transfer in turbulent convection with high-porosity fluid-saturated media, which is well above the onset of convection with Reynolds number far greater than unity.

Figure 2.

Darcy number $\xi$ dependence of Reynolds number $\varLambda$ (a) and Nusselt number $N$ (b) for various Rayleigh numbers $R=10^{11}$ (diamond), $10^{10}$ (cross), $10^9$ (up triangle), $10^8$ (down triangle), $10^7$ (square), $10^6$ (plus) and $10^5$ (circle).

The Reynolds number and the Nusselt number both arrive at a saturation regime with increasing Darcy number, as illustrated in figure 2. A phase transition towards the saturation regime of the Reynolds number can be observed; however, the saturation point is not readily identified with increasing Rayleigh number, though it does exist where the generalized Navier–Stokes equations gradually reduce to the equations of free fluids. A simple algorithm is adopted which is able to roughly determine the saturation point $\xi _{s}$ based on the estimation of curvature from the discrete data $\varphi (\xi _i)\ (i=1,2,\ldots )$ that approach a saturation regime:

(3.10)\begin{equation} \xi_{s}=\max\{\xi_{i} \mid C[{\rm lg}\, \varphi({\rm lg}\,\xi_{i})]\}, \end{equation}

where $C[\cdot ]$ is the usual curvature function and $\varphi$ represents the Reynolds number $\varLambda$ and the Nusselt number $N$. The saturation point $\xi _{s}$ is attained at the maximum value of $\xi _i$ that meets the criterion of (3.10). Consequently, the data are divided into two parts by $\xi _{s}$, above which $\varphi (\xi _i)$ are considered to fall into the saturation regime. Note that $\lg \varphi (\xi _i)$ is not directly obtained from the raw data, which have to be smoothed through a least-squares curve-fitting process. The operation removes the fluctuation of the raw data towards the saturation regime in order to obtain a monotonically increasing curve, which avoids the possibility of an unexpected incorrect saturation point. The saturated Darcy number $\xi _{s}$ for Reynolds number versus Rayleigh number $R$ is illustrated in figure 3, where a fitting line with a slope of $-2/3$ is included for reference. A scaling analysis can also be conducted based on the generalized Navier–Stokes equations and the heat transport equation in order to roughly estimate the saturation point of the Reynolds number. The characteristic velocity scale $u$ can be determined by balancing the Forchheimer quadratic inertia term with the buoyancy force term in (2.4) and the pressure gradient term in (2.1) in consideration of the fact that the nonlinear drag plays a major role when the flow state approaches the saturation regime:

(3.11)\begin{equation} u \sim \xi^{1/4}. \end{equation}

The transition occurs where the Forchheimer quadratic inertia term becomes less significant in comparison with the Brinkman term and the convective inertia term. The Brinkman term is used to obtain the constraint, which provides a criterion to identify the lower bound of the saturation regime:

(3.12)\begin{equation} \xi^{{-}1/2}u \sim \left(\frac{\sigma}{R}\right)^{1/2}l^{{-}2}. \end{equation}

Here $l$ is the typical length scale of the convective flow, which can be obtained by substituting equation (3.11) into (3.6):

(3.13)\begin{equation} l \sim R^{{-}1/2}\xi^{{-}1/4}, \end{equation}

where the Prandtl number $\sigma$ is assumed to be unity. The lower boundary of the saturation regime is determined by substituting equations (3.11) and (3.13) into (3.12):

(3.14)\begin{equation} R^{1/2}\xi^{3/4} \sim 1. \end{equation}

The relationship can also be obtained by comparing the Forchheimer quadratic inertia term with the convective inertia term. The Forchheimer quadratic inertia term and the Darcy term gradually become negligible with increasing Rayleigh number $R$ and Darcy number $\xi$, and the flow state eventually attains the saturation regime with $R^{1/2}\xi ^{3/4}\gg 1$. The constraint in (3.14) can be written in the form $\xi \sim R^{-2/3}$, which is consistent with the fitting line in figure 3. An increasing Rayleigh number $R$ results in a decreasing saturated Darcy number $\xi _{s}$, above which the convective state is driven away from the Forchheimer regime until it is practically indistinguishable from that of free fluids. Note that (3.14) and (3.4) actually represent the same scaling relationship, although they correspond to two parallel lines with a slope of $-2/3$ in the $(\xi ,R)$ phase diagram, where the former lies above the latter.

Figure 3.

Saturated Darcy number $\xi _{s}$ versus Rayleigh number $R$ for Reynolds number $\varLambda$ (plus) and Nusselt number $N$ (square).

The values of various force terms have been calculated to evaluate the effect of the Darcy number on the flow, including the convective inertia term and the Brinkman viscous term in the generalized momentum equation (2.1), the Forchheimer quadratic inertia term, the Darcy term and the buoyancy force term in (2.4). The peak values of these force terms in the fluid field of porous media are recorded, which are averaged over $100$ separate snapshots sampled in the final statistically steady states. The assessment is based on the maximum value in the fluid field of each force term (Lage Reference Lage1992), as these are more distinguishable in comparison with the spatio-temporally averaged values. For clarity, the peak values of various force terms $F_{\upsilon }$ in the generalized Navier–Stokes equation versus Darcy number are illustrated in figure 4 only for Rayleigh number $R=10^{11}$. The buoyancy force term remains constant over the whole range of Darcy numbers, which provides a baseline for evaluating the data of the other force terms. The Brinkman viscous term increases with Darcy number from the onset of convection and reaches the same order as the buoyancy force term for $\xi \geq 10^{-6}$. The convective inertia term is smaller than the buoyancy force term for $\xi \leq 10^{-6}$, while it increases with the Darcy number, which exceeds the value of the buoyancy force term at about $\xi =10^{-5}$ and becomes two orders of magnitude larger than the latter for $\xi \geq 10^{6}$. The Darcy term and the Forchheimer quadratic inertia term counteract the buoyancy term and retard the occurrence of saturation of the Reynolds number, both of which increase with the Darcy number from the onset of convection until reaching their maxima at $2\times 10^{-7}$ and $10^{-5}$, respectively. The Darcy term dominates the Brinkman viscous term and the two inertia terms for $\xi \leq 2\times 10^{-7}$, which indicates a regime that can be approximated by the Darcy model, whereas the deviation from the near-Darcian regime is apparent with increasing Darcy number. The Darcy term and the Forchheimer quadratic inertia term both decrease with increasing Darcy number; however, the former declines even more significantly in comparison with the latter for $\xi \geq 10^{-6}$. The Forchheimer quadratic inertia term becomes more than three orders of magnitude larger than the Darcy term even for $\xi \geq 10^{-2}$. The ratio of the Forchheimer drag to the Darcy drag can be expressed as

(3.15)\begin{equation} \lambda_{F}=c_F |\boldsymbol{u}|\left(\frac{\xi R}{\sigma}\right)^{1/2}, \end{equation}

wherein the absolute value of fluid velocity $|\boldsymbol {u}|$ increases with Rayleigh number (Bejan Reference Bejan1987). The nonlinear drag of the Forchheimer quadratic inertia term can be neglected in the regime of the Darcy model; however, the effect of the term has to be considered for the cases in which the Rayleigh number or the Darcy number is high enough. The retardation effect in porous media with high Darcy number is largely the result of the Forchheimer inertia term instead of the Darcy term, as the former increases with Rayleigh number. Porous media behave like free fluids in the saturated regime of Reynolds number, where the Forchheimer term and the Darcy term in (2.4) are negligible in comparison with the buoyancy force term.

Figure 4.

Peak values of the force terms $F_{\upsilon }$ in the generalized Navier–Stokes equation versus Darcy number $\xi$ at Rayleigh number $R=10^{11}$ for the buoyancy force term (cross), the Brinkman term (circle), the Darcy term (plus), the Forchheimer quadratic inertia term (square) and the convective inertia term (diamond).

The flow states in fluid-saturated porous media are displayed in the $(\xi ,R)$ phase diagram in figure 5. The Rayleigh number and the Darcy number are two collaborative factors which determine the intensity of convection. The figure is sampled by collecting the data points from figure 2. The lower boundaries are delimited with the minimal value of the Darcy number for different regimes; hence the data points correspond to the phase boundaries. A theoretically predicted line is plotted, $\xi =4{\rm \pi} ^2/R$, which provides the lower boundary for the regime with Nusselt numbers above unity, since the generalized Navier–Stokes model is not suitable for the analysis near the onset of convection where the Reynolds number is infinitesimally small. The near-Darcian regime is characterized by the linear drag of the Darcy term, which is an order of magnitude larger than the other force terms except the buoyancy force term in the generalized Navier–Stokes equation. The boundary line with a slope of $-2/3$ is plotted in accordance with (3.4) for Darcy's law to hold. The flow state comes into the Forchheimer regime in the case where the Forchheimer term is larger than the Darcy term, where the nonlinear drag effect dominates the flow in fluid-saturated porous media. The convective inertia regime can be located by those delimiting points, over which the convective inertia term is even larger than the buoyancy force term and the convection becomes highly unsteady and spatially inhomogeneous. The delimiting points roughly mark out the phase boundary, which exhibits the combined effect of Rayleigh number and Darcy number on the flow states.

Figure 5.

Phase diagram of the flow states in fluid-saturated porous media with lower boundaries for the near-Darcian regime (plus), the Forchheimer regime (circle) and the convective inertia regime (cross).

The asymptotic behaviour of the Nusselt number appears rather surprising, and differs markedly from that of the Reynolds number according to figure 2. A phase transition towards the saturated regime of the Nusselt number occurs at a smaller Darcy number $\xi _{s}=10^{-6}$ for Rayleigh number $R=10^{11}$ than that necessary for $R=10^{5}$, which is approximately unity. The saturation point decreases with increasing Rayleigh number, which is also illustrated in figure 3. A fitting line with a slope of $-1$ is provided for reference only, which, however, cannot be readily identified by a simple scaling analysis based on the generalized Navier–Stokes equations and the heat transport equation. The difference of the saturation point $\xi _s$ between Reynolds number $\varLambda$ and Nusselt number $N$ is rather small for $R=10^5$; nevertheless, the gap increases with Rayleigh number. This discovery is somewhat counterintuitive in consideration of the fact that the heat transfer is mainly conducted through fluid motion in porous media.

We are then concerned with the different dependence on the Darcy number of the Nusselt number and the Reynolds number. The Reynolds number is related to the extent of turbulence, while the Nusselt number is determined by the heat transfer across porous media. The saturated regime of the Nusselt number emerges prior to that of the Reynolds number, where the convection in porous media approaches that of free fluid flow. In contrast, the premature saturation of the Nusselt number can be ascribed to the counter-gradient heat transport by large-scale flows (Qiang et al. Reference Qiang, Cao, Li and Zhang2017). The Rayleigh–Bénard turbulent convection is dominated by the plumes that sprout from the thermal boundary layers, which go along with the large-scale flows in the convective layer (Hewitt et al. Reference Hewitt, Neufeld and Lister2012). Some of the plumes do not release their thermal energy completely during travel towards the opposite plates, and continue their motion in the reverse direction with the large-scale flows, which results in a counter-gradient heat flux (van der Poel, Stevens & Lohse Reference van der Poel, Stevens and Lohse2013; Huang & Zhou Reference Huang and Zhou2013). The counter-gradient heat fluxes increase not only with Darcy number, but also with Rayleigh number, which suppresses the effective global heat transport. The statistics of the velocity and temperature in the whole fluid field are illustrated by the joint probability density function of vertical velocity $u_\perp$ and temperature $T$. The logarithm of an instantaneous normalized probability distribution in the statistically steady states is plotted in figure 6. The quadrants are labelled by $+$ and $-$ with the velocity and temperature corresponding to the gradient and counter-gradient heat fluxes, respectively. The probability distribution $P(u_\perp ,T)$ is concentrated in the first and third quadrants for $R=10^{8}$ with $\xi =10^{-4}$ (figure 6a), where the ascending hot flows and descending cold flows dominate the fluid field that corresponds to the gradient heat fluxes. The probability distribution $P(u_\perp ,T)$ is negligible in the second and fourth quadrants, which correspond to the weak counter-gradient heat fluxes. Nonetheless, the counter-gradient heat fluxes become conspicuous for $R=10^{8}$ with $\xi =10^{6}$ (figure 6b), in which case $P(u_\perp ,T)$ is clearly observable in the second and fourth quadrants. The discrepancy even increases for the probability distribution with $R=10^{11}$ between $\xi =10^{-5}$ (figure 6c) and $\xi =10^{6}$ (figure 6d).

Figure 6.

Logarithm of the normalized joint probability density function of vertical velocity $u_{\perp }$ and temperature $T$ for different Rayleigh numbers and Darcy numbers: $R=10^{8}$, $\xi =10^{-4}$ (a); $R=10^{8}$, $\xi =10^6$ (b); $R=10^{11}$, $\xi =10^{-5}$ (c); $R=10^{11}$, $\xi =10^6$ (d). The quadrants with $u_{\perp }T>0$ and $u_{\perp }T<0$ are labelled by $+$ and $-$, which correspond to the gradient and counter-gradient heat fluxes, respectively.

The gradient and counter-gradient heat fluxes can be calculated via conditional averaging in the relevant area of the computational domain:

(3.16a,b)\begin{equation} N^{+}=\langle N \mid u_{{\perp}} T>0 \rangle, \quad N^{-}=\langle N \mid u_{{\perp}}T<0 \rangle. \end{equation}

The instantaneous vertical profiles of the gradient and counter-gradient heat fluxes averaged in the horizontal direction are presented in figure 7 in order to provide a direct view of their relative strengths for different Darcy numbers and Rayleigh numbers, where the distance between the two parallel conducting plates is scaled to unity. The solid line corresponds to the heat flux generated by hot plumes, while the dashed line is related to cold plumes. The positive and negative values of $N$ represent the gradient and counter-gradient heat fluxes, respectively, which fluctuate along the vertical direction until they approach zero at the thermal boundary layers. The amplitude of $N^-$ remains obviously smaller than that of $N^+$ for $R=10^{8}$ with $\xi =10^{-4}$ (see figure 7a), whereas the discrepancy diminishes for the same $R$ with $\xi =10^6$ (see figure 7b). A similar behaviour is also observed for the higher Rayleigh number $R=10^{11}$, where $N^-$ is negligible with $\xi =10^{-5}$ (see figure 7c). Nonetheless, the amplitudes of both $N^+$ and $N^-$ increase significantly with $\xi$, and the amplitude of $N^-$ becomes comparable to that of $N^+$ for $\xi =10^6$ (see figure 7d). The heat fluxes strongly fluctuate for high Rayleigh numbers, which can be ascribed to the discontinuous convective plumes with small scales (Qiang et al. Reference Qiang, Cao, Li and Zhang2017). The gradient and counter-gradient heat fluxes both increase not only with Rayleigh number, but also with Darcy number.

Figure 7.

Horizontally averaged gradient heat flux $N^+$ and counter-gradient heat flux $N^-$ for different Rayleigh numbers and Darcy numbers: (a) $R=10^{8}$, $\xi =10^{-4}$; (b) $R=10^{8}$, $\xi =10^6$; (c) $R=10^{11}$, $\xi =10^{-5}$; (d) $R=10^{11}$, $\xi =10^6$. The solid line corresponds to the heat flux generated by hot plumes and the dashed line to that generated by cold plumes.

A comparison is made between the gradient and counter-gradient heat fluxes with large-scale flows in order to shed more light on the issue of the premature saturation of the Nusselt number. The average has been conducted in the whole computational domain over 100 separate snapshots sampled in the final statistically steady states. For brevity, we present only illustrative results for $R=10^{11}$, with the two component parts $N^+$ and $N^-$ versus Darcy number $\xi$ exhibited in figure 8. The counter-gradient heat flux $N^-$ remains apparently smaller than the gradient heat flux $N^+$ over the whole range $10^{-6}\leq \xi \leq 10^6$; however, the discrepancy diminishes for high Darcy number. Component $N^-$ increases with $\xi$ more rapidly than $N^+$ for $\xi \geq 10^{-6}$, which coincides with the saturated Darcy number $\xi _{s}$ for the Nusselt number illustrated in figure 3. Generally, the gradient heat transport is compensated for by the counter-gradient heat transport to a great extent with increasing Darcy number, which attenuates the effects of large-scale flows on the global heat transfer especially for high Darcy numbers. The Nusselt number can also be obtained with the gradient and counter-gradient heat fluxes by

(3.17)\begin{equation} N=1+\gamma N^+{+}(1-\gamma)N^-, \end{equation}

where $\gamma$ and $1-\gamma$ represent the contents of the gradient and counter-gradient heat fluxes in the fluid layer and are obtained by counting the numbers of grid points in the computational domain with $u_{\perp }T>0$ and $u_{\perp }T<0$, respectively. The parameter $\gamma$ decreases with increasing $\xi$, which is illustrated in the inset of figure 8. Here $\gamma \geq 0.9$ is clearly observed for $\xi \leq 10^{-6}$, where the gradient heat flux dominates in the fluid layer, whereas $\gamma$ approaches 0.5 for $\xi =10^6$, where the content and intensity of the counter-gradient heat flux are comparable to those of the gradient heat flux. The parameter $1-\gamma$ also provides a measurement of the ratio of the quantity of plumes that do not lose their thermal energy completely after reversing their movement direction to the quantity of plumes that sprout from the thermal boundary layers, which is related to the instability mechanism operating in the thermal boundary layers (Qiang et al. Reference Qiang, Cao, Li and Zhang2017).

Figure 8.

Gradient heat flux $N^+$ (cross) and counter-gradient heat flux $N^-$ (plus) versus Darcy number for Rayleigh number $R=10^{11}$. Inset: Darcy number dependence of the content of the gradient heat flux.

The temperature fields for different Darcy numbers with $R=10^{11}$ are illustrated in figure 9, which are scaled by the height of the convective layer $h$ and the non-dimensional minimum and maximum of temperature ($-0.5\leq T\leq 0.5$). A large-scale circulation can be readily identified for Darcy number $\xi =10^{6}$ (see figure 9a), which resembles the flow pattern of Rayleigh–Bénard convection for free fluids. The hot plumes accumulate on one side near the lower boundary layer, while the cold plumes concentrate on the other side close to the upper boundary layer. The hot and cold plumes, which are laterally separated, participate in large-scale circulation towards opposite plates. Some of the plumes do not lose their thermal energy completely after reaching the opposite plates, and are dragged back towards their originating thermal boundary layers by the large-scale circulation. Nonetheless, the large-scale circulation is attenuated for Darcy number $\xi =10^{-2}$ (see figure 9b), such that the coherent structures of isolated plumes are replaced by continuous plume clusters. The large-scale quasiperiodic convective structure observed with high Darcy numbers can no longer be sustained for $\xi =10^{-5}$ (see figure 9c) and is broken down, while unsteady plumes that originate from the thermal boundary layer drive vigorous columnar exchange flows through the porous media. The flow pattern for Darcy number $\xi =10^{-6}$ is very reminiscent of the Darcy model with high Rayleigh–Darcy number (see figure 9d), although the models involved are rather different (Hewitt et al. Reference Hewitt, Neufeld and Lister2012). The wavenumber $k$ of the columnar flows is measured using a time average of the Fourier transform of the temperature field at the mid-height of the computational domain. We also obtained very similar results by applying the same method to the vertical velocity field. The columnar flows in the convection layer of porous media can be easily discerned for the near-Darcian regime, where the wavenumber can also be determined by directly counting the local maxima of the temperature and vertical velocity. The Rayleigh number $R$ dependence of the time-averaged horizontal wavenumber $k$ is presented in figure 10 for Darcy numbers $2\times 10^{-7}\leq \xi \leq 10^{-4}$. The wavenumber $k$ increases with Rayleigh number $R$, and a fitting line $k=0.005\times R^{0.387}$ is provided for $\xi =2\times 10^{-7}$. The wavenumber becomes smaller with increasing Darcy number for the same Rayleigh number, which can also be roughly observed in figure 9(c,d). The discrepancy of the wavenumber $k$ between different Darcy numbers $\xi$ becomes apparent with increasing Rayleigh number $R$, which means that the convective structures are rather different for $R\geq 5\times 10^8$. We note that the Nusselt number $N$ remains approximately invariant with Darcy number $10^{-6}\leq \xi \leq 10^{-4}$ for Rayleigh number $R=10^{11}$ in view of figure 2(b), which therefore indicates that the heat transfer is independent of the convective structures in the saturation regime. The phenomenon is consistent with the observation in free fluid flow, where heat transfer in turbulent convection is primarily carried out by thermal plumes and is not sensitive to the strength and structure of the convection, especially for high Rayleigh numbers (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009). The wavenumber $k$ is scaled by the Rayleigh–Darcy number $R_{\xi }^{0.4}$ in the inset of figure 10, which provides a comparison with the result of the Darcy model for reference only. The scaling exponent 0.387 is similar to the numerical value of $0.4$ obtained with the Darcy model (Hewitt et al. Reference Hewitt, Neufeld and Lister2012).

Figure 9.

A snapshot of the temperature field for different Darcy numbers (a) $\xi =10^{6}$, (b) $\xi =10^{-2}$, (c) $\xi =10^{-5}$ and (d) $\xi =10^{-6}$ with Rayleigh number $R=10^{11}$, illustrating the heat fluxes in fluid-saturated porous media.

Figure 10.

Rayleigh number $R$ dependence of the time-averaged horizontal wavenumber $k$ for different Darcy numbers: $\xi =10^{-4}$ (plus), $2\times 10^{-5}$ (circle), $10^{-5}$ (cross), $2\times 10^{-6}$ (diamond), $10^{-6}$ (up triangle) and $2\times 10^{-7}$ (square). Inset: time-averaged horizontal wavenumber $k$ scaled by the Rayleigh–Darcy number $R_{\xi }^{0.4}$.

Finally, we verify the relationship of the thickness of the thermal boundary layer with the Nusselt number:

(3.18)\begin{equation} \delta_T\simeq\frac{h}{2N}, \end{equation}

where the thickness of the thermal boundary layer $\delta _T$ is determined by the profile of the root mean square temperature fluctuation $\vartheta _T$, which is obtained from averaging in the horizontal direction over 100 separate snapshots in the final statistically steady states. Typical profiles of the value of $\vartheta _T$ as a function of the distance $z$ from the top plate are provided for Rayleigh number $R=10^{11}$ with different Darcy numbers in figure 11. Symmetrically, similar profiles near the bottom plate have also been obtained. The distance is normalized by the height of the convective layer $h$. The temperature fluctuation reaches a maximum at a position $\delta _T\simeq 0.0025$ from the top and bottom plates over the range of Darcy numbers $10^{-6}\leq \xi \leq 10^6$, which is consistent with the results for the Nusselt number in figure 2 with reference to (3.18). The amplitude of $\vartheta _T$ at distance $\delta _T=0.0025$ from the plates decreases with Darcy number, which can be ascribed to reduced entrainment of the plumes by the lateral flushing of the large-scale flows into the interior of the convective layer. The phenomenon also accords with the behaviour of the Reynolds number in figure 2(a), which indicates that the intensity of the turbulence decreases with Darcy number. The position of the maximum temperature fluctuation shifts to the interior of the convective layer for $\xi <10^{-6}$, where the increasing thermal boundary layer $\delta _T$ corresponds to the smaller Nusselt number $N$ outside the saturation regime in figure 2(b). Consequently, the heat transfer of the turbulent convection in fluid-saturated porous media is related to behaviour of the thermal boundary layer instead of the strength and structure of the large-scale flows, which have different effects on the Nusselt number and Reynolds number. The effect of the counter-gradient heat transport $N^-$ is as a regulator that adaptively compensates for the gradient heat transport, which results in the Nusselt number being irrelevant to the Reynolds number.

Figure 11.

Root mean square temperature profile versus the normalized distance from the plate with Rayleigh number $R=10^{11}$ for Darcy numbers $\xi =10^6$ (circle), $10^{-2}$ (cross), $10^{-4}$ (up triangle), $10^{-5}$ (down triangle), $10^{-6}$ (square), $2\times 10^{-7}$ (plus) and $10^{-7}$ (diamond).