2 Numerical model

2.1 Governing equations

We consider a Boussinesq fluid in a 2-D RB cell with unity aspect ratio, containing an obstacle array located on a square lattice, as depicted in figure 1. The fluid is heated from below and cooled from above with a temperature difference $\unicode[STIX]{x1D6E5}$ between two horizontal plates a distance $L$ apart in the vertical direction. The diameter of the circular, solid obstacles is $D$. The minimum value of the obstacle–obstacle and obstacle–wall separations is $l$. The porous structures are characterized by the porosity $\unicode[STIX]{x1D719}$, which measures the volume fraction of the fluid phase. A value of $\unicode[STIX]{x1D719}=1$ corresponds to traditional RB convection without the obstacle array. For the convection cell with unity aspect ratio, we have $l=(L-N\times D)/(N+1)$, where $N$ is the number of the obstacles along the horizontal and vertical directions. The corresponding porosity is $\unicode[STIX]{x1D719}=1-N^{2}\unicode[STIX]{x03C0}D^{2}/(4L^{2})$. For fixed $D$, only specific values of $\unicode[STIX]{x1D719}$ (or $l$) can be achieved, since $N$ is an integer. We note that, besides the case with regular obstacle arrangement on a square lattice, the effect of obstacle arrangement will also be considered in this study (see figure 11). The non-dimensional governing equations describing the flow dynamics in the pores read

(2.1)$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}{\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\sqrt{{\displaystyle \frac{Pr}{Ra}}}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+T\boldsymbol{e}_{z}+\boldsymbol{f},\\ {\displaystyle \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}T)={\displaystyle \frac{1}{\sqrt{PrRa}}}\unicode[STIX]{x1D6FB}^{2}T,\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ is the velocity vector with components $(u,v)$ along the horizontal and vertical directions $(x,z)$, $p$ is the pressure, $T$ is the temperature and $\boldsymbol{e}_{z}$ is the unit vector in the vertical direction. The immersed boundary force term $\boldsymbol{f}$ is added to account for the presence of the obstacle array. The governing equations were non-dimensionalized using $L$ for length, $\unicode[STIX]{x1D6E5}$ for temperature, the free-fall velocity $U=\sqrt{g\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}L}$ for velocity and $L/U$ for time, where $g$ is the gravitational acceleration and $\unicode[STIX]{x1D6FD}$ the isobaric expansion coefficient.

The two dimensionless parameters in the governing equations (2.1) are the Rayleigh number, $Ra=g\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}L^{3}/(\unicode[STIX]{x1D708}\unicode[STIX]{x1D705})$, and the Prandtl number, $Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$, where $\unicode[STIX]{x1D708}$ is the kinematic viscosity and $\unicode[STIX]{x1D705}$ the thermal diffusivity. The heat transfer property is measured by the Nusselt number, $Nu=\sqrt{RaPr}\langle vT\rangle _{x,t}-\langle \unicode[STIX]{x2202}_{z}T\rangle _{x,t}$, where $\langle \cdot \rangle _{x,t}$ denotes taking averages over any horizontal plane and time. In practice, the average is taken over the top and bottom plates. No-slip and no-penetration boundary conditions are imposed at all solid surfaces, including the fluid–obstacle interfaces. The horizontal top and bottom plates are isothermal and the sidewalls are thermally insulated.

Note that the heat transfer between the fluid and solid phases is considered. For simplicity, we assume the same thermal properties for the two phases, including $\unicode[STIX]{x1D70C}c_{p}$, thermal conductivity $k$ and thermal diffusivity $\unicode[STIX]{x1D705}=k/(\unicode[STIX]{x1D70C}c_{p})$, where $\unicode[STIX]{x1D70C}$ is the density and $c_{p}$ the specific heat capacity. The temperature equation for the fluid and solid phases reads

(2.2)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}_{cp}T)=\frac{1}{\sqrt{PrRa}}\unicode[STIX]{x1D6FB}^{2}T,\end{eqnarray}$$

where $\boldsymbol{u}_{cp}$ (to be defined below) is the velocity for the fluid and solid phases depending on the position in the domain (Verzicco Reference Verzicco2002, Reference Verzicco2004; Stevens, Lohse & Verzicco Reference Stevens, Lohse and Verzicco2014; Ardekani et al. Reference Ardekani, Abouali, Picano and Brandt2018a,Reference Ardekani, Al Asmar, Picano and Brandtb; Sardina et al. Reference Sardina, Brandt, Boffetta and Mazzino2018).

2.2 Numerical methods

In this subsection the essential elements of the numerical methods are presented. The governing equations are discretized using a second-order finite-difference method with a pressure-correction scheme. A fractional-step third-order Runge–Kutta scheme is employed for the time stepping of the explicit terms and a Crank–Nicolson scheme for the implicit terms. A uniform, staggered, Cartesian grid is used in this study. The pressure Poisson equation is solved efficiently using a fast-Fourier-transform-based solver. For more details of the numerical schemes of the governing equations, we refer the reader to van der Poel et al. (Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015).

To account for the obstacle array, a direct-forcing immersed boundary method (IBM) in the Euler–Lagrange framework is adopted (Uhlmann Reference Uhlmann2005; Breugem Reference Breugem2012). For the resolution of the obstacles, $N_{L}$ markers are distributed uniformly along the boundary of each obstacle, with the Lagrangian grid size being about 0.7 times the Eulerian grid size. In each Runge–Kutta substep, a first prediction velocity is obtained by advancing the momentum equations in time without considering the IBM force $\boldsymbol{f}$. The first prediction velocity is then interpolated from the Eulerian grid to the Lagrangian grid using the moving-least-squares approach (Vanella & Balaras Reference Vanella and Balaras2009; de Tullio & Pascazio Reference de Tullio and Pascazio2016; Spandan et al. Reference Spandan, Meschini, Ostilla-Mónico, Lohse, Querzoli, de Tullio and Verzicco2017, Reference Spandan, Lohse, de Tullio and Verzicco2018). The IBM force $\boldsymbol{f}$ required on each Lagrangian marker for satisfying no-slip and no-penetration conditions is computed, which is then spread back to the Eulerian grid using the moving-least-squares approach. The force $\boldsymbol{f}$ is used to update the velocity, followed by a standard pressure correction step. We note that the pressure field required to enforce the incompressibility condition has vanishing gradient at the immersed boundary. Thus, the pressure correction step does not change practically the numerical accuracy of the IBM scheme, as observed by Fadlun et al. (Reference Fadlun, Verzicco, Orlandi and Mohd-Yusof2000) and Kempe & Fröhlich (Reference Kempe and Fröhlich2012). As a validation of the implementation of the IBM, we considered a steady, axisymmetric shear flow inside a circular domain driven by the rotation of the circular boundary at a fixed rotation rate without body force. The numerical simulation showed that the azimuthal velocity increases linearly with the distance away from the rotation axis and the vorticity is uniform, in quantitative agreement with the theoretical results.

The heat transfer between the fluid and solid phases is realized by solving the temperature equation (2.2) in both phases (Ardekani et al. Reference Ardekani, Abouali, Picano and Brandt2018a,Reference Ardekani, Al Asmar, Picano and Brandtb; Sardina et al. Reference Sardina, Brandt, Boffetta and Mazzino2018). A phase indicator $\unicode[STIX]{x1D709}$ is introduced to represent the solid volume fraction and to distinguish the fluid and solid phases on each grid point of two components of velocity. The value of $\unicode[STIX]{x1D709}$ is determined from a level-set function $\unicode[STIX]{x1D701}$ given by the signed distance of four corner nodes to the obstacle surface using the formula (Ardekani et al. Reference Ardekani, Abouali, Picano and Brandt2018a; Kempe & Fröhlich Reference Kempe and Fröhlich2012)

(2.3)$$\begin{eqnarray}\unicode[STIX]{x1D709}=\frac{\displaystyle \mathop{\sum }_{n=1}^{4}-\unicode[STIX]{x1D701}_{n}{\mathcal{H}}(-\unicode[STIX]{x1D701}_{n})}{\displaystyle \mathop{\sum }_{n=1}^{4}|\unicode[STIX]{x1D701}_{n}|},\end{eqnarray}$$

where ${\mathcal{H}}$ is the Heaviside step function. Then the velocity of the combined phase is defined at each point in the computational domain as

(2.4)$$\begin{eqnarray}\boldsymbol{u}_{cp}=(1-\unicode[STIX]{x1D709})\boldsymbol{u}_{f}+\unicode[STIX]{x1D709}\boldsymbol{u}_{p},\end{eqnarray}$$

where $\boldsymbol{u}_{f}$ and $\boldsymbol{u}_{p}$ denote velocities of the fluid and solid phases, respectively. For fixed obstacles considered in this study we obviously have $\boldsymbol{u}_{p}=0$. We note that the coupling of the IBM and heat transfer between different phases has been considered in other convection problems, such as turbulent convection confined by walls with finite conductivities (Verzicco Reference Verzicco2002, Reference Verzicco2004; Stevens et al. Reference Stevens, Lohse and Verzicco2014).

Table 1.

Numerical details of typical grid resolutions. The columns from left to right denote the Rayleigh number $Ra$, the number of obstacles $N_{o}$, the porosity $\unicode[STIX]{x1D719}$, the grid resolution $N_{x}\times N_{z}$, the Nusselt number $Nu$, the number of grid nodes $N_{BL}$ in the thermal boundary layers and the number of grid nodes $N_{D}$ per obstacle diameter. The thickness of the thermal boundary layer is estimated as $\unicode[STIX]{x1D6FF}_{th}=L/(2Nu)$. We note that the table only presents some typical cases. Totally, 67 different cases were simulated, and the results are collected in figure 3.

Figure 3.

(a) Variations of $Nu$ with $Ra$ for different $\unicode[STIX]{x1D719}$. (b) The compensated plot of (a). (c) Results for different $\unicode[STIX]{x1D719}$ collapse by rescaling the thermal boundary layer thickness $\unicode[STIX]{x1D6FF}_{th}=L/(2Nu)$ with the pore scale $l$. (d) Variation of $Nu$ with the Darcy Rayleigh number $Ra^{\ast }=RaDa$.

In this study a uniform Eulerian grid is used and the resolution is chosen to fully resolve the boundary layers and the smallest scales in the bulk (Stevens, Verzicco & Lohse Reference Stevens, Verzicco and Lohse2010; Shishkina et al. Reference Shishkina, Stevens, Grossmann and Lohse2010; Zhang, Zhou & Sun Reference Zhang, Zhou and Sun2017). The numerical details of some typical grid resolutions are given in table 1. For the small-$Ra$ cases, the grid is restricted by the resolution of the obstacles. For these cases a mesh with a grid of $540\times 540$ is used, and the number of grid nodes per obstacle diameter is 22. For larger $Ra$, the mesh size chosen is adequate to achieve a full resolution of the thermal boundary layer, with the thermal boundary layer resolved by at least 10 grid points. For the highest $Ra$ ($Ra=10^{10}$), a mesh with a grid of $3072\times 3072$ is used. This ensures that the turbulent flow in the bulk and boundary layers is fully resolved. For the traditional RB convection at $Ra=10^{10}$, the typical Kolmogorov length scale $\unicode[STIX]{x1D702}$ is estimated by the global criterion $\unicode[STIX]{x1D702}=LPr^{1/2}/[Ra(Nu-1)]^{1/4}$, and the Batchelor scale $\unicode[STIX]{x1D702}_{B}$ is estimated by $\unicode[STIX]{x1D702}_{B}=\unicode[STIX]{x1D702}Pr^{-1/2}$ (Shishkina et al. Reference Shishkina, Stevens, Grossmann and Lohse2010). We find that the grid spacing $\unicode[STIX]{x1D6E5}_{g}$ satisfies $\unicode[STIX]{x1D6E5}_{g}\lesssim 0.16\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D6E5}_{g}\lesssim 0.33\unicode[STIX]{x1D702}_{B}$. The thermal boundary layers are resolved with 14 grid points, which agrees with the recommendations of Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010). Long-time averages are conducted for the calculation of $Nu$. The difference of $Nu$ obtained by time averages over the first and second halves of the simulations (both taken after initial transients) was smaller than 1 %. Besides, the relative difference between the top- and bottom-wall $Nu$ was smaller than 1 %.

Simulations were performed for varying Rayleigh numbers (from $10^{5}$ to $10^{10}$) and porosities (from 0.75 to 1). The Prandtl number and obstacle diameter were fixed at $Pr=4.3$ and $D=0.04$, respectively.

3 Heat transfer

In this section the influence of the obstacle array on the heat transfer properties is examined. The presence of the obstacle array has a significant influence on the heat transfer properties. The normalized Nusselt number $Nu(\unicode[STIX]{x1D719})$ for two fixed $Ra$ is shown in figure 2. In figure 3, $Nu(Ra)$ for five fixed porosities $\unicode[STIX]{x1D719}$ is shown. In the traditional RB convection with $\unicode[STIX]{x1D719}=1$, $Nu$ increases with $Ra$ following an effective power law $Nu\sim Ra^{0.30}$, consistent with the results in the literature for 2-D RB convection (van der Poel et al. Reference van der Poel, Stevens and Lohse2013; Zhang et al. Reference Zhang, Zhou and Sun2017; Zhu et al. Reference Zhu, Stevens, Verzicco and Lohse2017). For fixed $\unicode[STIX]{x1D719}<1$, it is found that the variation of $Nu$ with $Ra$ exhibits two scaling regimes for the parameter range studied. Compared with the results of traditional RB convection with $\unicode[STIX]{x1D719}=1$, the heat transfer is reduced for small $Ra$ even if $Nu$ increases with $Ra$ with a steep effective power law $Nu\sim Ra^{0.65}$; on the other hand, for large enough $Ra$, $Nu$ becomes larger than the corresponding value without obstacles and increases with $Ra$ with an effective power law $Nu\sim Ra^{0.30}$, similar to that of traditional RB convection. The increase of $Nu$ in the large-$Ra$ regime is more visible in the compensated plot of figure 3(b) and in figure 2. Note that the regime of heat transfer enhancement with decreasing $\unicode[STIX]{x1D719}$ displayed in figure 2 corresponds to the large-$Ra$ regime shown in figure 3. Figure 3(a,b) shows that the trends of variation of $Nu$ with $Ra$ for different values of $\unicode[STIX]{x1D719}$ are similar, and the critical $Ra$ for scaling crossover increases as $\unicode[STIX]{x1D719}$ is decreased. Considering that the values of $\unicode[STIX]{x1D719}$ and $Nu$ determine the characteristic length scales of the regular porous medium and flow structures, respectively, the scaling crossover indicates a competition of the length scale of porosity and that of the flow structures. When we rescale the thermal boundary layer thickness $\unicode[STIX]{x1D6FF}_{th}=L/(2Nu)$ with the obstacle separation $l$, it is found that the results for different $\unicode[STIX]{x1D719}$ approximately collapse, and the scaling crossover occurs when $\unicode[STIX]{x1D6FF}_{th}/l\approx 1$, as shown in figure 3(c). The data collapse indicates that the heat transfer properties are determined by the length scales of the flow structures and the pore scale of the regular porous medium.

It is interesting to examine the relevance of our results to pure Darcy-type convection. For this, it is instructive to show how $Nu$ varies with the so-called Darcy Rayleigh number $Ra^{\ast }=RaDa$, where $Da=K/L^{2}$ is the Darcy number. Here $K$ is the permeability, which depends on the geometry of the porous media and represents the resistance for the fluid to flow through the porous media. We estimate $Da$ with porosity $\unicode[STIX]{x1D719}$ and obstacle diameter $D$ based on Kozeny’s equation (Nithiarasu et al. Reference Nithiarasu, Seetharamu and Sundararajan1997; Nield & Bejan Reference Nield and Bejan2006):

(3.1)$$\begin{eqnarray}Da=\frac{\unicode[STIX]{x1D719}^{3}D^{2}}{150(1-\unicode[STIX]{x1D719})^{2}L^{2}}.\end{eqnarray}$$

The results are shown in figure 3(d). We find that in the small-$Ra^{\ast }$ regime, $Nu(Ra^{\ast })$ shows a trend to collapse as $\unicode[STIX]{x1D719}$ is decreased, which is a signature of the transition to Darcy-type flow and consistent with the observation based on the generalized non-Darcy model (Nithiarasu et al. Reference Nithiarasu, Seetharamu and Sundararajan1997). We note that in porous media there may exist multiple flow states with different heat transfer efficiencies for the same parameters. The solution may not be able to sample all the flow configurations due to the suppression of fluctuations in the presence of an obstacle array. The numerically realized steady state may be different for different initial conditions. Despite this, we expect that the differences of the statistics of multiple flow states are relatively small, particularly for the large-$Ra$ cases, and the global trend of variation of $Nu$ with $Ra$ and $\unicode[STIX]{x1D719}$ will not be qualitatively affected by the existence of multiple flow states.

Figure 4.

The joint PDFs of the vertical velocity $v$ and the temperature fluctuation$\unicode[STIX]{x1D6FF}T=T-T_{m}$ on the horizontal mid-plane $z=0.5$ for different $\unicode[STIX]{x1D719}$ at $Ra=10^{8}$: (a) $\unicode[STIX]{x1D719}=1$, (b) $\unicode[STIX]{x1D719}=0.92$ and (c) $\unicode[STIX]{x1D719}=0.82$. Corresponding values of the cross-correlation $C$ between $v$ and $\unicode[STIX]{x1D6FF}T$ are labelled at the top right-hand corner of each panel.

We now focus on the statistics of the velocity and temperature fields on the horizontal mid-plane $z=0.5$, to quantify the modifications of heat transfer properties as $\unicode[STIX]{x1D719}$ is decreased, particularly in the regime where $Nu(\unicode[STIX]{x1D719})$ varies non-monotonically. Figure 4 plots the joint probability density function (PDF) $P(v,\unicode[STIX]{x1D6FF}T)$ of vertical velocity $v$ and temperature fluctuation $\unicode[STIX]{x1D6FF}T$ for different $\unicode[STIX]{x1D719}$ at $Ra=10^{8}$, where $\unicode[STIX]{x1D6FF}T=T-T_{m}$ and $T_{m}=0.5$ is the arithmetic mean temperature. In the absence of the obstacle array (i.e. for $\unicode[STIX]{x1D719}=1$), the velocity can attain considerably larger values than that in porous media at the same $Ra$. The temperature in the bulk is well mixed with $\unicode[STIX]{x1D6FF}T\approx 0$. Thus, $P(v,\unicode[STIX]{x1D6FF}T)$ is concentrated in a slender region located on the horizontal axis, as shown in figure 4(a). From $P(v,\unicode[STIX]{x1D6FF}T)$ the correlation between $v$ and $\unicode[STIX]{x1D6FF}T$ can be identified. Due to buoyancy, flow elements with positive (negative) $\unicode[STIX]{x1D6FF}T$ are more likely to move upward (downward). Note that there is a non-negligible probability for the occurrence of counter-gradient convective heat transfer, which is manifested in the finite values of $P(v,\unicode[STIX]{x1D6FF}T)$ in the second and fourth quadrants, consistent with the observation of counter-gradient convective heat transfer in Sugiyama et al. (Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010) and Huang & Zhou (Reference Huang and Zhou2013). The properties of temperature fluctuations and convective heat transfer can also be quantified by the normalized PDFs of $\unicode[STIX]{x1D6FF}T$ and $v\cdot \unicode[STIX]{x1D6FF}T$ on the horizontal mid-plane $z=0.5$, which are given in figure 5. The strong mixing of temperature in the bulk is manifested by the peak value of $P(\unicode[STIX]{x1D6FF}T)$ at $\unicode[STIX]{x1D6FF}T=0$ in figure 5(a), and the correlation between $v$ and $\unicode[STIX]{x1D6FF}T$ is reflected by the asymmetry of the left- and right-hand branches of $P(v\cdot \unicode[STIX]{x1D6FF}T)$ in figure 5(b). The counter-gradient convective heat transfer is quantified by the left-hand branch of $P(v\cdot \unicode[STIX]{x1D6FF}T)$.

As $\unicode[STIX]{x1D719}$ is decreased, it is found that the fluctuation properties are significantly influenced, as shown in figure 4(b,c). The velocity distribution is narrowed down, showing that the convection strength is decreased for the flow through regular porous media. Compared with the results of traditional RB convection with $\unicode[STIX]{x1D719}=1$, the temperature mixing is less efficient and large temperature fluctuations are more likely to appear. The decreased mixing efficiency is also manifested by the flattening of the peak of $P(\unicode[STIX]{x1D6FF}T)$ around $\unicode[STIX]{x1D6FF}T=0$ in figure 5(a). From $P(v,\unicode[STIX]{x1D6FF}T)$ it is found that the counter-gradient convective heat flux is suppressed, which is also manifested by the rapid decrease of the left-hand branch of $P(v\cdot \unicode[STIX]{x1D6FF}T)$ as $v\cdot \unicode[STIX]{x1D6FF}T$ is decreased from 0 in figure 5(b). Note that as $\unicode[STIX]{x1D719}$ is decreased, the distribution of $P(v,\unicode[STIX]{x1D6FF}T)$ is more compact and concentrated along a line with positive slope, which suggests that the correlation between $v$ and $\unicode[STIX]{x1D6FF}T$ is enhanced in regular porous media. We here quantify the correlation between $v$ and $\unicode[STIX]{x1D6FF}T$ using the cross-correlation $C\equiv \langle g_{v}g_{T}\rangle$, where $g_{\unicode[STIX]{x1D712}}\equiv [\unicode[STIX]{x1D712}-\langle \unicode[STIX]{x1D712}\rangle ]/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D712}}$, $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D712}}\equiv \sqrt{\langle \unicode[STIX]{x1D712}^{2}\rangle -\langle \unicode[STIX]{x1D712}\rangle ^{2}}$ and $\langle \cdot \rangle$ denotes the average over the horizontal mid-plane and time. Values of $C=1$ and $-1$ correspond to perfect correlations and anti-correlations, respectively, and $C=0$ corresponds to no correlation. The values of $C$ for the three cases of figure 4 are shown in the top right-hand corner of each panel, which definitely demonstrate the enhancement of the correlation between $v$ and $\unicode[STIX]{x1D6FF}T$ with decreasing $\unicode[STIX]{x1D719}$. For both the traditional RB convection and convection in regular porous media, the cross-correlation between temperature fluctuation $\unicode[STIX]{x1D6FF}T$ and horizontal velocity $u$ on the horizontal mid-plane is small with $|C|<0.1$.

Figure 5.

The normalized PDFs of (a) $\unicode[STIX]{x1D6FF}T/\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FF}T}$ and (b) $v\cdot \unicode[STIX]{x1D6FF}T/\unicode[STIX]{x1D70E}_{v\cdot \unicode[STIX]{x1D6FF}T}$ on the horizontal mid-plane $z=0.5$ for different $\unicode[STIX]{x1D719}$ at $Ra=10^{8}$, where $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FF}T}$ and $\unicode[STIX]{x1D70E}_{v\cdot \unicode[STIX]{x1D6FF}T}$ are the standard deviations of $\unicode[STIX]{x1D6FF}T$ and $v\cdot \unicode[STIX]{x1D6FF}T$, respectively. The dashed lines indicate the Gaussian distribution.

The variation of $P(v,\unicode[STIX]{x1D6FF}T)$ with $\unicode[STIX]{x1D719}$ shown in figure 4 leads us to identify two competing effects of the obstacle array on the heat transfer. On the one hand, the flow becomes more coherent with enhanced correlation between temperature fluctuation and vertical velocity. This is beneficial for heat transfer. On the other hand, the obstacle array slows down the convection due to enhanced resistance, which is unfavourable for heat transfer. The competing effects of the obstacle array on the heat transfer result in the non-monotonic variation of $Nu$ with $\unicode[STIX]{x1D719}$ as was already shown in figure 2. As $\unicode[STIX]{x1D719}$ is slightly decreased from 1, the enhancement of the flow coherence increases the heat transfer efficiency. While, when $\unicode[STIX]{x1D719}$ is small enough, the convection is strongly suppressed, leading to heat transfer reduction. The heat transfer enhancement on increasing flow coherence is a widespread phenomenon and has been observed in various flow configurations, such as confined RB convection (Huang et al. Reference Huang, Kaczorowski, Ni and Xia2013; Chong et al. Reference Chong, Wagner, Kaczorowski, Shishkina and Xia2018) and partitioned RB convection (Bao et al. Reference Bao, Chen, Liu, She, Zhang and Zhou2015). Chong et al. (Reference Chong, Yang, Huang, Zhong, Stevens, Verzicco, Lohse and Xia2017) revealed a universal, non-monotonic behaviour of turbulent transport in systems with competing stabilizing and destabilizing forces. The non-monotonic variation of $Nu$ with the increase of the stabilizing effect of porous structure is consistent with the observation of Chong et al. (Reference Chong, Yang, Huang, Zhong, Stevens, Verzicco, Lohse and Xia2017).

Figure 6.

Nusselt number $Nu(t)$ in a time interval of 200 dimensionless units for different $\unicode[STIX]{x1D719}$ at (a) $Ra=10^{7}$ and (b) $Ra=10^{9}$.

Figure 7.

Typical snapshots of the instantaneous (a,d,g) temperature $T$, (b,e,h) the velocity magnitude $|\boldsymbol{u}|$ and (c,f,i) the local convective heat flux $v\cdot \unicode[STIX]{x1D6FF}T$ in the vertical direction at (a–c) $(Ra,\unicode[STIX]{x1D719})=(10^{9},1)$, (d–f) $(Ra,\unicode[STIX]{x1D719})=(10^{9},0.82)$ and (g–i) $(Ra,\unicode[STIX]{x1D719})=(10^{7},0.82)$. Circles in (d–i) indicate the obstacle array. Note that according to the temperature equation (2.2), the temperature in the obstacles is well defined.

Figure 6 shows sampled time records of $Nu$ for different $\unicode[STIX]{x1D719}$ at $Ra=10^{7}$ and $10^{9}$, showing the influence of the obstacle array on the temporal characteristics of fluctuations. It is found that for relatively small $Ra$, high-frequency fluctuations are significantly suppressed as $\unicode[STIX]{x1D719}$ is decreased, while for large enough $Ra$, the suppression of high-frequency fluctuations is less visible for the values of $\unicode[STIX]{x1D719}$ studied. The influence of the obstacle array on the fluctuations is dependent on the relative magnitudes of the spatial coherence length and the pore scale. At small $Ra$, the length scales of the flow structures are large compared to the pore scale, and the resistance of the obstacle array to the convection is strong. Thus flow fluctuations are significantly suppressed, and there is no developed turbulence in the pores. As $Ra$ is increased, the convection is more energetic and the length scales of the flow structures are decreased. Thus the obstacle array has smaller effects on the fluctuations, and the flow in the pores becomes chaotic or even turbulent at large enough $Ra$.

4 Flow structure

In this section we focus on the flow structures of RB convection in regular porous media. Figure 7 displays typical snapshots of the instantaneous temperature $T$, the velocity magnitude $|\boldsymbol{u}|$ and the convective heat flux $v\cdot \unicode[STIX]{x1D6FF}T$ in the vertical direction. In the traditional RB convection with $\unicode[STIX]{x1D719}=1$, the flow consists of a well-organized large-scale circulation (LSC) with two well-established counter-rotating corner rolls (Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010), as shown in figure 7(a–c). Due to the formation of the LSC, the velocity in the bulk is non-uniform. Following the LSC, plumes detaching from the thermal boundary layers mainly go up and down near the sidewalls. The temperature at the cell centre is well mixed with $T\approx T_{m}$. From the distribution of $v\cdot \unicode[STIX]{x1D6FF}T$ it is confirmed that counter-gradient convective heat flux appears locally around the LSC and corner rolls, which is due to the bulk dynamics and the competition between the corner-flow rolls and the LSC (Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010; Huang & Zhou Reference Huang and Zhou2013). The inclusion of an array of obstacles has a significant influence on the flow structures, as displayed in figure 7(d–i). Convection strength is reduced and the LSC is suppressed due to the impedance of the obstacle array. Temperature mixing at the cell centre is less efficient. Thermal plumes detaching from the thermal boundary layers can penetrate deep in the bulk, forming convection channels, namely the regions with strong flows that are formed between the obstacles. When $Ra$ is relatively large, the characteristic length scales of the flow structures are smaller than the pore scale, and the flow is chaotic and dominated by fragmented plumes, as shown in figure 7(d); while when $Ra$ is relatively small the flow is dominated by large-scale plumes, as shown in figure 7(g). Comparing figures 7(f,i) and 7(c), it is found that the counter-gradient convective heat flux is reduced in regular porous media, which is attributed to the suppression of the LSC and coherence of plume dynamics due to the impedance of the obstacle array.

Snapshots of temperature fields at $Ra=10^{7}$ and different porosities are shown in figure 8, demonstrating the enhancement of flow coherence with decreasing $\unicode[STIX]{x1D719}$. For relatively large $\unicode[STIX]{x1D719}$, plumes wander randomly in the bulk, as shown in figure 8(a). As $\unicode[STIX]{x1D719}$ is further decreased, plumes cannot move freely due to the impedance of the obstacle array and prefer moving along convection channels, as shown in figure 8(b,c).

Figure 8.

Snapshots of the instantaneous temperature field at $Ra=10^{7}$ and different $\unicode[STIX]{x1D719}$, corresponding to the filled symbols in figure 2: (a) $\unicode[STIX]{x1D719}=0.99$, (b) $\unicode[STIX]{x1D719}=0.90$ and (c) $\unicode[STIX]{x1D719}=0.85$.

To further investigate the influence of the obstacle array on the convection, we study the velocity distribution in the bulk region. Figure 9 shows the PDFs of the velocity magnitude $|\boldsymbol{u}|$ inside a circular domain of radius $r$ at the cell centre for different $\unicode[STIX]{x1D719}$ at $Ra=10^{8}$. In the traditional RB convection with $\unicode[STIX]{x1D719}=1$, the velocity magnitude is spatially non-uniform due to the formation of the LSC. The mean velocity vanishes at the cell centre and increases linearly with $r$ at leading order. Due to this non-uniformity, $P(|\boldsymbol{u}|)$ shifts along the $|\boldsymbol{u}|$ axis as $r$ is increased, and follows a linear relation for small $|\boldsymbol{u}|$, $P(|\boldsymbol{u}|)\sim |\boldsymbol{u}|$, as shown in figure 9(a). The obstacle array has a significant influence on the velocity distribution, as shown in figure 9(b,c). The maximum convection velocity is decreased due to the drag of the obstacle array, and the results for different values of $r$ collapse, indicating that the velocity distribution is spatially uniform in the bulk. The shape of $P(|\boldsymbol{u}|)$ also changes qualitatively and deviates from the power-law distribution as $\unicode[STIX]{x1D719}$ is decreased. In regular porous media the probability density is larger for $|\boldsymbol{u}|$ to take smaller values, and when $\unicode[STIX]{x1D719}$ is small enough, $P(|\boldsymbol{u}|)$ satisfies an exponential distribution, as shown in figure 9(c). The change of velocity distribution demonstrates the transition of flow organization as $\unicode[STIX]{x1D719}$ is decreased, from the formation of the LSC in the traditional RB convection to the wandering motion of plumes penetrating in the pores in the porous-media convection.

Figure 9.

Log–log plots of the PDFs of the velocity magnitude $|\boldsymbol{u}|$ inside a circular domain of radius $r$ at the cell centre for different $\unicode[STIX]{x1D719}$ at $Ra=10^{8}$: (a) $\unicode[STIX]{x1D719}=1$, (b) $\unicode[STIX]{x1D719}=0.92$ and (c) $\unicode[STIX]{x1D719}=0.82$. The insets are the same results shown in semi-log plots. In (a) and the inset of (c) guiding lines of constant slope are included, describing the respective exponential scaling behaviours of $P(|\boldsymbol{u}|)$. As $|\boldsymbol{u}|$ increases, $P(|\boldsymbol{u}|)$ increases linearly for small $|\boldsymbol{u}|$ in (a), and it decreases exponentially in the intermediate range of $|\boldsymbol{u}|$ in (c).

Figure 10.

The joint PDFs of two velocity components $(u,v)$ inside a circular domain of radius $r=0.4$ at the cell centre for different $\unicode[STIX]{x1D719}$ at $Ra=10^{8}$: (a) $\unicode[STIX]{x1D719}=1$, (b) $\unicode[STIX]{x1D719}=0.92$ and (c) $\unicode[STIX]{x1D719}=0.82$.

Figure 10 shows the joint PDFs of two velocity components ($u,v$) inside a circular domain of radius $r=0.4$ at the cell centre for different $\unicode[STIX]{x1D719}$ at $Ra=10^{8}$. It is found that both the size and shape of $P(u,v)$ change as $\unicode[STIX]{x1D719}$ is decreased from 1. The shrinking of $P(u,v)$ indicates the reduction of convection strength. The change of contour levels manifests the increased probability for $|\boldsymbol{u}|$ to take smaller values. When $\unicode[STIX]{x1D719}$ is small enough, the contour lines of $P(u,v)$ protrude along the horizontal and vertical coordinate axes, indicating the formation of convection channels.

How robust are our results with respect to the arrangement of the obstacles? To answer this question, we study how the obstacle arrangement influences the flow structure. We consider three different arrangements with similar porosities, i.e. two regular arrangements and one random arrangement. Simulations are performed at $Ra=10^{9}$. The snapshots of temperature $T$, velocity magnitude $|\boldsymbol{u}|$ and convective heat flux $v\cdot \unicode[STIX]{x1D6FF}T$ in the vertical direction are shown in figure 11. It is seen that the flow structure is significantly influenced by the obstacle arrangement. In figure 11(a–c), the porous medium I is constructed by placing the circular obstacles on a square lattice, as done in the rest of this paper. As already discussed, thermal plumes can penetrate up and down without obstruction, forming vertical convection channels with strong fluid and heat transport. In figure 11(d–f), the regular porous medium II is constructed by rotating the lattice of porous medium I by $45^{\circ }$. In such an obstacle arrangement, the vertical convection channels observed in porous medium I are obstructed by the obstacles, except the vertical channels close to the sidewalls. Despite that, we find that thermal plumes can penetrate deep into the bulk along zigzag channels aligning vertically. Besides, strong fluid and heat transport occur close to the sidewalls, where the fluid can penetrate up and down without obstruction. In figure 11(g–i), a random porous medium is constructed (porous medium III), with the pore scale $l\geqslant l_{min}=0.005$, where $l_{min}$ is the minimum pore scale. It is observed that the plume motion is less organized compared to those in the regular porous media I and II. Curved convection channels with strong flow emerge in the pores, with no preferred flowing directions. Except the localized vertical channels, the curved convection channels are less efficient at transporting heat in the vertical direction. The Nusselt numbers for the obstacle arrangements I, II and III are 51.96, 53.97 and 47.14, respectively. The lower $Nu$ in porous medium III is attributed to the less organized plume motion and less efficient vertical heat transfer.

Figure 11.

Snapshots of (a,d,g) temperature $T$, (b,e,h) velocity magnitude $|\boldsymbol{u}|$ and (c,f,i) convective heat flux $v\cdot \unicode[STIX]{x1D6FF}T$ in the vertical direction for three different obstacle arrangements with similar porosities at $Ra=10^{9}$. (a–c) Porous medium I consists of obstacles located on a square lattice, and the porosity $\unicode[STIX]{x1D719}=0.68$. (d–f) Porous medium II is constructed by rotating the obstacle lattice of I by $45^{\circ }$, with porosity $\unicode[STIX]{x1D719}=0.67$. (g–i) Porous medium III is constructed by randomly placing the obstacles in the cell, with porosity $\unicode[STIX]{x1D719}=0.68$.

5 Energy dissipation rates

Rayleigh–Bénard convection is a closed system with exact balance properties. The overall kinetic and thermal energy dissipation rates $\unicode[STIX]{x1D716}_{u,T}$ are related to the governing and response parameters by two exact relations (Shraiman & Siggia Reference Shraiman and Siggia1990; Grossmann & Lohse Reference Grossmann and Lohse2000; Ahlers et al. Reference Ahlers, Grossmann and Lohse2009):

Figure 12.

Typical snapshots of the kinetic energy dissipation rates $\text{log}_{10}\unicode[STIX]{x1D716}_{u}(\boldsymbol{x})$ on logarithmic scale for different $Ra$ and $\unicode[STIX]{x1D719}$: (a,d) $Ra=10^{6}$, (b,e) $Ra=10^{8}$ and (c,f) $Ra=10^{9}$; (a–c) $\unicode[STIX]{x1D719}=1$ and (d–f) $\unicode[STIX]{x1D719}=0.92$.

(5.1a,b)$$\begin{eqnarray}\langle \unicode[STIX]{x1D716}_{u}\rangle _{V,t}=\frac{\unicode[STIX]{x1D708}^{3}}{L^{4}}(Nu-1)RaPr^{-2},\quad \langle \unicode[STIX]{x1D716}_{T}\rangle _{V,t}=\unicode[STIX]{x1D705}\frac{\unicode[STIX]{x1D6E5}^{2}}{L^{2}}Nu,\end{eqnarray}$$

where $\langle \cdot \rangle _{V,t}$ denotes the volume and time average, and

(5.2a,b)$$\begin{eqnarray}\unicode[STIX]{x1D716}_{u}=\frac{1}{2}\unicode[STIX]{x1D708}\mathop{\sum }_{ij}\left[\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}\right]^{2},\quad \unicode[STIX]{x1D716}_{T}=\unicode[STIX]{x1D705}\mathop{\sum }_{i}\left[\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}x_{i}}\right]^{2}.\end{eqnarray}$$

We note that the exact global relations (5.1) are not affected by the presence of the obstacles. Figures 12 and 13 plot typical snapshots of the non-dimensional kinetic and thermal energy dissipation rates $\log _{10}\unicode[STIX]{x1D716}_{u,T}(\boldsymbol{x})$ on logarithmic scale for $Ra=10^{6},10^{8},10^{9}$ and $\unicode[STIX]{x1D719}=1,0.92$. In the traditional RB convection with $\unicode[STIX]{x1D719}=1$, it is found that the characteristic length scales of the flow structures become smaller for larger $Ra$, as shown in figures 12(a–c) and 13(a–c). A well-organized LSC develops at sufficiently large $Ra$, and velocity boundary layers appear close to both the horizontal plates and vertical sidewalls, which dominate the dissipation of kinetic energy. There are no thermal boundary layers along the sidewalls due to the imposed adiabatic condition for temperature, and the dissipation of thermal energy is dominated by plumes in the bulk and the thermal boundary layers near the horizontal top and bottom plates. The bulk contributions of both the kinetic and thermal energy dissipation rates are small (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001, Reference Grossmann and Lohse2004; Zhang et al. Reference Zhang, Zhou and Sun2017). Results of convection in a regular porous medium are displayed in figures 12(d–f) and 13(d–f). When $Ra$ is small, convection is dominated by large-scale hot and cold plumes, penetrating to the top and bottom plates, respectively, along convection channels with fast velocities. Along these channels intense dissipation of kinetic energy occurs with length scale characterized by the pore scale. As $Ra$ is increased, the length scales of the flow structures become smaller, and plumes wander in the bulk, leading to intense dissipation of kinetic energy in the bulk, in distinct contrast to the results of traditional RB convection with $\unicode[STIX]{x1D719}=1$. For sufficiently large $Ra$, the kinetic energy dissipation in the bulk is concentrated around the obstacle surfaces due to the interaction between plumes and the obstacle array, and the pore scale is no longer a proper length scale for the kinetic energy dissipation. Due to the wandering motion of plumes in regular porous media, intense dissipation of thermal energy can occur locally in the bulk. However, the averaged bulk contribution of thermal energy dissipation remains small as compared to the boundary layer contribution.

Figure 13.

Typical snapshots of the thermal energy dissipation rates $\text{log}_{10}\unicode[STIX]{x1D716}_{T}(\boldsymbol{x})$ on logarithmic scale for different $Ra$ and $\unicode[STIX]{x1D719}$: (a,d) $Ra=10^{6}$, (b,e) $Ra=10^{8}$ and (c,f) $Ra=10^{9}$; (a–c) $\unicode[STIX]{x1D719}=1$ and (d–f) $\unicode[STIX]{x1D719}=0.92$. In (e,f), we show the obstacles as thin white circles. Note that the temperature field in them is well defined according to (2.2), and correspondingly the thermal energy dissipation rate $\unicode[STIX]{x1D716}_{T}(\boldsymbol{x},t)$ is non-zero at these locations.

To further quantify the spatial distributions of the kinetic and thermal energy dissipation rates, we plot the profiles of averaged energy dissipation rates $\langle \unicode[STIX]{x1D716}_{u,T}\rangle _{x,t}$ along the vertical direction $z$ in figure 15, showing the contributions of the bulk and boundary layer regions. Using (5.1) $Nu$ can be obtained by integrating these profiles of $\langle \unicode[STIX]{x1D716}_{u,T}\rangle _{x,t}$. Figure 15(a,b) plots the results of traditional RB convection with $\unicode[STIX]{x1D719}=1$. The contributions of the boundary layer region dominate for both the kinetic and thermal energy dissipation rates, consistent with the observations in figures 12 and 13. From figure 15(a) it seems to be less obvious that the kinetic energy dissipation is dominated by regions close to the walls. We note that $\langle \unicode[STIX]{x1D716}_{u}\rangle _{x,t}$ contains the kinetic energy dissipation from the velocity boundary layers near the sidewalls. Contrary to the velocity, the temperature field has no boundary layers close to the sidewalls, which are adiabatic. The results of convection in regular porous media are shown in figure 15(c,d). For these cases the obstacle array has a significant influence on the distribution of $\langle \unicode[STIX]{x1D716}_{u}\rangle _{x,t}$ due to the no-slip condition at the pores, while $\langle \unicode[STIX]{x1D716}_{T}\rangle _{x,t}$ is only mildly affected. The relative contribution of the bulk region of $\unicode[STIX]{x1D716}_{u}$ is significantly enhanced, while $\unicode[STIX]{x1D716}_{T}$ remains dominated by the boundary layer. For relatively large $Ra$, intense dissipation of kinetic energy occurs around the obstacle surfaces. These observations are also consistent with those seen in figures 12 and 13.

Figure 14.

Variations of the averaged (a) kinetic and (b) thermal energy dissipation rates $\langle \unicode[STIX]{x1D716}_{u,T}\rangle _{x,t}$ with $z$ for different $\unicode[STIX]{x1D719}$ at $Ra=10^{6}$.

Figure 15.

Variations of the averaged (a,c) kinetic and (b,d) thermal energy dissipation rates $\langle \unicode[STIX]{x1D716}_{u,T}\rangle _{x,t}$ with $z$ for different $Ra$ at (a,b) $\unicode[STIX]{x1D719}=1$ and (c,d) $\unicode[STIX]{x1D719}=0.92$. In (c,d) the grey areas indicate obstacle positions.

Figure 16.

(a) Plots of $Re(Ra)$ in the small-$Ra$ regime. (b) Comparison of $Nu$ obtained based on the heat flux averaging $Nu=-\langle \unicode[STIX]{x2202}_{z}T\rangle _{x,t}$ over the horizontal plates (dashed lines) and (5.3) with $c=8.0$ (solid lines with symbols) for the same four cases for the porosity $\unicode[STIX]{x1D719}$ shown in (a).

In figure 2 it is shown that the heat transfer is enhanced when $\unicode[STIX]{x1D719}$ is slightly decreased from 1 at fixed $Ra$, and further decrease of $\unicode[STIX]{x1D719}$ reduces the heat transfer compared with the $\unicode[STIX]{x1D719}=1$ case. To show the influence of the obstacle array on the heat transfer from a local perspective, we plot the profiles of $\langle \unicode[STIX]{x1D716}_{u,T}\rangle _{x,t}$ for various $\unicode[STIX]{x1D719}$ at $Ra=10^{6}$ in figure 14. As $\unicode[STIX]{x1D719}$ is slightly decreased from 1, the kinetic energy dissipation is enhanced in the bulk, which is the local manifestation of heat transfer enhancement due to the obstacle array. When $\unicode[STIX]{x1D719}$ is further decreased, the pore scale decreases accordingly and convection is significantly slowed down due to the drag of the obstacle array. Correspondingly, $\langle \unicode[STIX]{x1D716}_{u}\rangle _{x,t}$ is decreased in both the boundary layer and the bulk regions, and $\langle \unicode[STIX]{x1D716}_{T}\rangle _{x,t}$ is also decreased, particularly in the boundary layer region. Due to the smallness of the spatial coherence length and strong convection at large $Ra$, smaller porosity (or pore scale) is needed to significantly slow down the convection and reduce the overall energy dissipation rates.

For fixed $\unicode[STIX]{x1D719}$, $Nu$ increases with $Ra$ with an effective power law $Nu\sim Ra^{0.65}$ when $Ra$ is relatively small, and it changes to $Nu\sim Ra^{0.30}$ for large enough $Ra$. For small $Ra$ with the spatial coherence length larger than the pore scale, turbulence is suppressed in the pores. In the presence of the obstacle array one basically gets additional laminar-type boundary layers in the bulk region, namely laminar shear flows in the pores, and viscosity dominates such that the flow becomes of Prandtl–Blasius–Pohlhausen type. Intense dissipation of kinetic energy occurs along the convection channels. Thus, the overall kinetic energy dissipation rate can be estimated as $\langle \unicode[STIX]{x1D716}_{u}\rangle _{V,t}\sim \unicode[STIX]{x1D719}\unicode[STIX]{x1D708}U_{rms}^{2}/l^{2}$, corresponding to the so-called $\infty$-regime in Grossmann & Lohse (Reference Grossmann and Lohse2001) for large $Pr$ and small $Ra$, but with $l$ replacing the box size $L$. Based on the exact relation between $\langle \unicode[STIX]{x1D716}_{u}\rangle _{V,t}$ and $Nu$, we can further estimate $Nu$ as

(5.3)$$\begin{eqnarray}Nu\approx c\cdot \unicode[STIX]{x1D719}\left(\frac{L}{l}\right)^{4}Pr^{2}Re^{2}Ra^{-1}+1,\end{eqnarray}$$

where $Re=U_{rms}l/\unicode[STIX]{x1D708}$ and $c$ is an undetermined constant (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001, Reference Grossmann and Lohse2004). Plots of $Re(Ra)$ in the small-$Ra$ regime are shown in figure 16(a) for four fixed $\unicode[STIX]{x1D719}$, and then $Nu$ can be estimated using (5.3). With $c=8.0$, it is found that $Nu$ obtained using (5.3) is consistent with the results based on the definition of $Nu$ for different $Ra$ and $\unicode[STIX]{x1D719}$, as shown in figure 16(b), which provides additional confirmation of our observations regarding the modifications of the heat transfer and flow structure by the porous media. As $Ra$ is increased, the characteristic scales of the flow structures become smaller and intense dissipation of kinetic energy tends to occur around the obstacle surfaces. For large enough $Ra$, convection in the pores becomes chaotic or even turbulent, the pore scale is no longer a proper length scale for the kinetic energy dissipation rate and the estimate of $Nu$ through (5.3) is no longer applicable.

6 Summary

We have studied the heat transfer and flow structure of 2-D RB convection in regular porous media using pore-scale modelling. The porous medium is comprised of circular, solid obstacles located on a square lattice. The heat transfer between the fluid and solid phases is considered.

The obstacle array has two competing effects on the heat transfer. On the one hand, the flow becomes more coherent with the correlation between temperature fluctuation and vertical velocity enhanced and the counter-gradient convective heat transfer suppressed, leading to heat transfer enhancement. On the other hand, the convection strength is reduced due the impedance of the obstacle array, leading to heat transfer reduction. The coexistence of the two competing effects leads to the non-monotonic behaviour of $Nu(\unicode[STIX]{x1D719})$ as $\unicode[STIX]{x1D719}$ is decreased from 1, as shown in figure 2. The heat transfer enhancement is consistent with the counterintuitive observation that an appropriate strength of a stabilizing force can enhance heat transfer by increasing flow coherence (Chong et al. Reference Chong, Yang, Huang, Zhong, Stevens, Verzicco, Lohse and Xia2017). Significant enhancement of the heat transfer due to the increased flow coherence was also observed in the confined inclined convection in low-$Pr$ fluids (Zwirner & Shishkina Reference Zwirner and Shishkina2018). Due to the emergence of system-sized plume columns and the interaction of these plume columns with the opposed boundary layers in inclined convection, an increase of the heat transfer by a factor of approximately 2.3 can be realized. The observations of heat transfer enhancement in these distinct systems are manifestations of a universal mechanism to enhance turbulent transfer by increasing flow coherence.

The influence of porosity on flow properties is dependent on $Ra$, and two different heat transfer regimes are observed at fixed $\unicode[STIX]{x1D719}$. In the small-$Ra$ regime where viscosity dominates, $Nu$ is decreased as compared to the $\unicode[STIX]{x1D719}=1$ case with a steep effective power law $Nu\sim Ra^{0.65}$, while in the large-$Ra$ regime $Nu$ is increased as compared to $\unicode[STIX]{x1D719}=1$ with the classical power law $Nu\sim Ra^{0.30}$. The scaling crossover occurs when the thickness of the thermal boundary layer $\unicode[STIX]{x1D6FF}_{th}$ is comparable to the pore scale $l$.

The influence of the obstacle array on the heat transfer is also analysed from the local perspective, namely the energy dissipation rates. It is found that the bulk contribution of kinetic energy dissipation rate is enhanced as $\unicode[STIX]{x1D719}$ is slightly decreased from 1; while when $\unicode[STIX]{x1D719}$ is small enough, convection is significantly slowed down by the obstacle array, and the kinetic energy dissipation rate is decreased in the whole cell. For small $Ra$, due to the large coherence length of the flow structures and the impedance of the obstacle array, turbulence is suppressed in the pores and additional laminar-type boundary layers appear in the bulk. It is shown that also for porous-media convection $Nu$ can be estimated from $Re$ through (5.3) similarly to that suggested by Grossmann & Lohse (Reference Grossmann and Lohse2001, Reference Grossmann and Lohse2004) for the viscosity-dominated so-called $\infty$-regime. Regarding the turbulent modulations, the results of this study suggest that it is possible to modulate the heat transfer with porous structure and that the transitional $Ra$ between different heat transfer regimes can be changed by manipulating the pore scale.

The influence of the obstacle arrangement on the flow structure is also studied. Three different arrangements with similar porosities are considered. The flow structure is found to be significantly influenced by the obstacle arrangement. In regular porous media, thermal plumes can penetrate deep into the bulk along convection channels aligning vertically, resulting in strong fluid and heat transport. When the obstacles are randomly located in the cell, the plume motion is less organized, and curved convection channels with strong flow emerge in the pores, with reduced efficiency of vertical heat transfer.

In this study we have mainly focused on the heat transfer and flow structure of 2-D RB convection in regular porous media. In the future, it will be of interest to extend the study to the 3-D case, to allow for a one-to-one comparison with experiment. Following the comparative study of Chong et al. (Reference Chong, Yang, Huang, Zhong, Stevens, Verzicco, Lohse and Xia2017) and Lim et al. (Reference Lim, Chong, Ding and Xia2019), it is also desirable to investigate further turbulent flows with different stabilizing forces, to get an even more complete understanding of various stabilizing–destabilizing flow systems.