1. Introduction
Thermally driven convective fluid motions represent a fundamental process widely observed in both geophysical and astrophysical systems, and they are essential for numerous industrial applications. Classical models, such as Rayleigh–Bénard convection (RBC) (Bodenschatz, Pesch & Ahlers Reference Bodenschatz, Pesch and Ahlers2000; Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Zhang et al. Reference Zhang, Xia, Zhou and Chen2020; Ecke & Shishkina Reference Ecke and Shishkina2023; Lohse & Shishkina Reference Lohse and Shishkina2024) and vertical natural convection (Ng et al. Reference Ng, Ooi, Lohse and Chung2015; Shishkina Reference Shishkina2016), have been extensively studied to uncover the fundamental mechanisms governing heat and momentum transport. In RBC, the fluid is heated from below and cooled from above, while in vertical natural convection, the fluid is confined between two vertical walls with differing temperatures. These models provide critical insights into the dynamics of buoyancy-driven flows and their roles in large-scale environmental and engineering systems.
Vertical natural convection holds significant practical importance, with applications ranging from thermal insulation in double-pane windows and horizontal heat transport in water reservoirs to crystal growth processes and the cooling of electronic devices. Furthermore, it serves as an idealised model for studying thermally driven atmospheric and oceanic circulations, providing a framework for understanding large-scale environmental flow phenomena (Lappa Reference Lappa2009). Despite extensive research being devoted to vertical natural convection, questions remain regarding the evolution of flow and heat transfer as the Rayleigh number
$Ra$
increases across various flow regimes.
A key research question in vertical natural convection is the scaling exponents of Nusselt number
$\textit{Nu}\sim R{{a}^{\gamma _1 }}$
and Reynolds number
$R{{e}}\sim R{{a}^{\gamma _2 }}$
. Research studies conducted by Trias et al. (Reference Trias, Soria, Oliva and Pérez-Segarra2007, Reference Trias, Gorobets, Soria and Oliva2010), Ng et al. (Reference Ng, Ooi, Lohse and Chung2015, Reference Ng, Spandan, Verzicco and Lohse2020), Shishkina (Reference Shishkina2016) and Wang et al. (Reference Wang, Xia, Yan, Sun and Wan2019, Reference Wang, Liu, Verzicco, Shishkina and Lohse2021) report that the scaling exponent
$\gamma _1$
varies between 1/4 and 1/3, while
$\gamma _2$
varies between 1/2 and 4/9, depending on the specific ranges of
$Ra$
and Prandtl number
$\textit{Pr}$
. For
$\textit{Pr}\gg 1$
, Shishkina (Reference Shishkina2016) applied laminar boundary layer theory to theoretically derive the scaling exponents
$\gamma _1 ={1}/{4}$
and
$\gamma _2 ={1}/{2}$
. These findings were corroborated by Ng et al. (Reference Ng, Spandan, Verzicco and Lohse2020) through direct numerical simulations (DNS) of three-dimensional vertical convection with periodic boundary conditions. Consequently, within the small range of
$Ra$
,
$\gamma _1$
tends to approximate to 1/4 (Trias et al. Reference Trias, Soria, Oliva and Pérez-Segarra2007, Reference Trias, Gorobets, Soria and Oliva2010; Wang et al. Reference Wang, Xia, Yan, Sun and Wan2019). However, as
$Ra$
increases,
$\gamma _1$
approaches approximately 0.31 (Ng et al. Reference Ng, Ooi, Lohse and Chung2015). In the larger range of
$Ra$
,
$\gamma _1$
is approximately 1/3, while
$\gamma _2$
is approximately 4/9 (Wang et al. Reference Wang, Liu, Verzicco, Shishkina and Lohse2021), which aligns with the results of RBC and is consistent with Grossman–Lohse (GL) theory (Grossmann & Lohse Reference Grossmann and Lohse2000; Lohse & Xia Reference Lohse and Xia2010).
Compared with the pure thermal convection in the above studies, porous-medium thermal convection is receiving increased attention. Investigating thermal convection in porous media is essential due to its wide-ranging applications in engineering and environmental systems (Nield & Bejan Reference Nield and Bejan2006; Rebhi et al. Reference Rebhi, Menni, Lorenzini and Ahmad2022), including the design of heat exchangers, geothermal energy extraction and cooling systems for electronic devices (Odenthal et al. Reference Odenthal, Steinmann, Eck and Laing2014). Pore-scale turbulent convection enhances heat transfer by actively driving fluid through the porous structure, thereby optimising thermal performance in systems where heat transfer alone is insufficient. As industries increasingly demand higher energy efficiency and precision, the investigation of porous-medium thermal convection becomes imperative for effective heat management.
The scaling laws of porous-medium thermal convection represent a fundamental area of research. Currently, two primary methods are employed to investigate this question: DNS (Liu et al. Reference Liu, Jiang, Chong, Zhu, Wan, Verzicco, Stevens, Lohse and Sun2020, Reference Liu, Jiang, Wang and Sun2021; Li, Yang & Sun Reference Li, Yang and Sun2024) and macroscopic models (Pirozzoli et al. Reference Pirozzoli, De Paoli, Zonta and Soldati2021; De Paoli Reference De Paoli2023; Hu & Yang Reference Hu and Yang2024; Zhu, Fu & De Paoli Reference Zhu, Fu and De Paoli2024). However, the existing studies predominantly focus on porous-medium RBC. The mechanisms underlying porous-medium vertical natural convection, especially with pore-scale turbulence, remain inadequately understood. Consequently, this study aims to address the following primary questions: (i) What are the major regimes of porous-medium vertical natural convection? (ii) What are the scaling laws of
$\textit{Nu}$
and
$\textit{Re}$
? (iii) Are there unified scaling models for
$\textit{Nu}(Ra_D, \phi , Da)$
and
$\textit{Re}(Ra_D, \phi , Da)$
? Here,
$Ra_D, \phi$
and
$Da$
refer to the Rayleigh–Darcy number, porosity and Darcy number, respectively.
This study employs theoretical analysis and DNS to investigate porous-medium vertical natural convection across a wide range of
$R{{a}_{D}}$
,
$\phi$
and
$Da$
. We identify three distinct regimes based on the characteristics of the flow states. Subsequently, we theoretically derive scaling laws for
$\textit{Nu}$
and
$\textit{Re}$
within each regime. Furthermore, we propose unified scaling models at
$\textit{Pr}=1$
. The models are validated by DNS results across a broad spectrum of
$(R{{a}_{D}},\phi , Da)$
. Through this study, we aim to enhance the understanding of regime transitions and the development of unified scaling laws in porous-medium vertical natural convection.
This paper is organised as follows: § 2 discusses the governing equations, the geometry set-ups and the numerical methods employed. In § 3, we present the corresponding simulation results and engage in a discussion of these findings. The posed questions are comprehensively addressed. Finally, § 4 presents the conclusions drawn from this study.
2. Governing equations, basic set-ups and numerical methods
2.1. Governing equations
The governing equations for DNS consist of the incompressible Navier–Stokes equations and the thermal transport equation, as presented in (2.1a – c )
\begin{equation} {\partial {{u}_{i}}}/{\partial {{x}_{i}}} =0, \\[-12pt]\nonumber \end{equation}
\begin{equation} {\partial {{u}_{i}}}/{\partial t} +{\partial \big ( {{u}_{i}}{{u}_{\!j}} \big )}/{\partial {{x}_{\!j}}} =-{{\rho }^{-1}}{\partial p}/{\partial {{x}_{i}}} +\nu {{{\partial }^{2}}{{u}_{i}}}/{\partial x_{\!j}^{2}} -{{g}_{i}}\beta (T-{{T}_{0}}), \\[-12pt]\nonumber \end{equation}
\begin{equation} {\partial T}/{\partial t} +\partial ({{u}_{i}}{ T})/{\partial {{x}_{i}}} ={\alpha {{\partial }^{2}}T}/{\partial x_{i}^{2}}. \end{equation}
Here,
${u}_{i}$
is the velocity,
$p$
is the pressure,
$T$
is the temperature with reference value
$T_0$
,
$\beta$
is the thermal expansion coefficient,
${g}_{i}$
is the applied gravitational acceleration,
$\nu$
is the kinematic viscosity and
$\alpha$
is the thermal diffusivity. The Oberbeck–Boussinesq approximation, assuming variations in density
$\rho \approx {{\rho }_{0}} [ 1-\beta ( T-{{T}_{0}} ) ]$
, are negligible except in the buoyancy term, is adopted.
2.2. The set-ups of geometry and important parameters
As illustrated in figure 1, the two-dimensional square cavity with size
$H$
is filled with a porous medium. The porous medium consists of a
$N_r\times N_r$
array of representative elementary volumes (REVs) with size
$S=H/N_r$
, which accurately captures the pore-scale properties of the entire domain. Each REV features a solid square element with size
$d$
at its centre, which adheres to no-slip and adiabatic (i.e. impermeable in the context of solutal convection) boundary conditions, resulting in porosity
$\phi =1-{{ ({d}/{S} )}^{2}}$
. The gravitational acceleration
$\boldsymbol{g}=(0,-g,0)$
acts downward, directing the buoyancy forces of the hot fluid upwards. The thermal boundary conditions for the top and bottom walls are adiabatic, while the left and right walls are maintained at constant temperatures of
$T_0+\Delta /2$
and
$T_0-\Delta /2$
, respectively.
In relation to the porous medium within the cavity, the Darcy number is defined as
\begin{equation} Da={K}/{{{H}^{2}}}, \end{equation}
where
$K$
represents permeability, quantifying the ability of the porous medium to facilitate fluid flow. Permeability is determined using Darcy’s law
$K={u}_{Di}\nu /g_i$
, where
$u_{Di}$
is the volume-averaged velocity in a REV for Stokes flow (Rao & Zhang Reference Rao and Zhang2025c
). It is important to note that
$K$
is calculated directly through DNS rather than relying on empirical equations, thus enhancing its accuracy. Regarding the presence of porous medium, the Rayleigh–Darcy number
$R{{a}_{D}}$
is employed instead of the conventional Rayleigh number
$Ra={g\beta \Delta {{H}^{3}}}/{ ( \nu \alpha )}$
and is defined as
\begin{equation} R{{a}_{D}}=Ra\boldsymbol{\cdot }Da. \end{equation}
The Nusselt number is defined as the ratio of total heat transfer to static conductive heat transfer at the conducting boundary in the fluid, and is expressed as follows:
\begin{equation} \textit{Nu}=\left ({{ {\langle {{\partial }_{1}}T \rangle _A}/{{{\langle {{\partial }_{1}}T_{Ra=0} \rangle _A }}} }}\right )|_{x_1=\pm H/2}, \end{equation}
where
$\langle \boldsymbol{\cdot }\rangle _A|_{x_1=\pm H/2}$
represents the time- and area-averaged value at position
$x_1=\pm H/2$
. The global Reynolds number is defined as
\begin{equation} \textit{Re}={ U_{\!D} H}/{\nu }, \end{equation}
where the root-mean-square velocity
$U_{\!D}=\phi U$
and
$U=\sqrt {\langle {u_iu_i}\rangle _{\!f}}$
. Here,
${\left \langle \boldsymbol{\cdot }\right \rangle }_{\!f}$
denotes the time- and intrinsic volume-averaged value for the fluid phase. The Prandtl number
$\textit{Pr}={\nu }/{\alpha }$
is defined as the ratio of momentum diffusivity (viscosity) to thermal diffusivity. In this study, we mainly focus on the scaling laws with respect to
$R{{a}_{D}}$
,
$\phi$
and
$Da$
. Therefore, we fix the Prandtl number at 1 for most simulation cases.
The schematic geometry and set-ups employed in DNS.
2.3. Numerical methods
In this study, the finite volume method (FVM) and lattice Boltzmann method (LBM) are employed for simulations. The solver pimpleFoam, based on OpenFOAM, has been modified to effectively couple the flow and temperature dynamics. A second-order central difference scheme is utilised for the spatial discretisation of both velocity and temperature, achieved through the interpolation of variables at the interfaces of grid cells. The Poisson equation is solved at each time step to update the pressure, while the pressure-implicit with splitting of operators scheme ensures accurate coupling between pressure and velocity. The LBM is employed to indirectly solve microscopic equations by determining particle distributions. The fundamental equation in LBM is a discretised version of the Boltzmann equation, with the collision operator addressed by the Bhatnagar–Gross–Krook model. Both the FVM and LBM solvers have undergone rigorous validation and verification in prior research (Rao, Kuznetsov & Jin Reference Rao, Kuznetsov and Jin2020; Rao Reference Rao2022; Rao & Jin Reference Rao and Jin2022; Rao & Zhang Reference Rao and Zhang2025a
,
Reference Rao and Zhangb
,
Reference Rao and Zhangc
); therefore, they are directly utilised in this study. All cases are calculated using two distinct mesh resolutions, employing
${({60S}/{d}+1)}^{3}$
grid points for each REV for laminar flow and
${({90S}/{d}+1)}^{3}$
grid points for each REV for turbulent flow. This methodology aims to achieve a balance between mesh independence and computational cost. The computational duration is established at a minimum of
$100$
time periods
$t={H}/{{{u}_{m}}}$
, where
${{u}_{m}}={g\beta \Delta K}/{\nu }$
is the characteristic velocity. Furthermore, the calculated Nusselt number is consistent for the hot and cold plates.
3. Results and discussion
3.1. Flow states
To answer the questions raised in § 1, we conducted a series of DNS simulations across a range of
$R{{a}_{D}}$
from 0 to
$10^8$
. As illustrated in figure 2(a–c), which shows time-averaged temperature fields, our observations reveal three distinct regimes, characterised by the angle
$\theta$
between the isothermal layer
$\overline {T}=T_0$
and the
$x_1$
axis. It is important to note that these three regimes pertain specifically to the bulk region. For clarity, in figure 2(d–f) we depict the velocity vectors and amplitude contours by mapping the time- and REV-averaged velocity
$\overline {u}_{Di}$
onto a pure cavity of identical dimensions. Here,
$\overline {\,\boldsymbol{\cdot }\,}$
denotes the time average.
(a,b,c) Time-averaged temperature
$\overline {T}$
contour with isothermal layer
$\overline {T}=T_0$
and (d,e, f) filtered time- and REV-averaged velocity vectors and magnitude contour in three regimes at
$\phi =0.69$
,
$Da=2.1\times {{10}^{-5}}$
. Panels show (a,d)
$R{{a}_{D}}=2.1\times {{10}^{-2}}$
; (b,e)
$R{{a}_{D}}=4.2\times {{10}^{1}}$
; (c, f)
$R{{a}_{D}}=2.1\times {{10}^{4}}$
.
The zoom-in snapshots of instantaneous turbulent structures near the bottom-right corner in the turbulent regime. Panels show (a)
$R{{a}_{D}}=2.1\times {{10}^{4}}$
,
$\phi =0.69$
; (b)
$R{{a}_{D}}=1.0\times {{10}^{8}}$
,
$\phi =0.69$
; (c)
$R{{a}_{D}}=1.0\times {{10}^{8}}$
,
$\phi =0.31$
.
Regime I (steady regime, distinct vertical isothermal layer): heat transfer is predominantly governed by conduction. The flow is characterised by a steady state with a small velocity amplitude, which allows for the applicability of the Stokes approximation in the momentum equation. The viscous term nearly balances the buoyancy term, resulting in a weak circulation around the cavity, as illustrated in figure 2(d). In this regime, the temperature distribution is almost linear across the REVs. Subsequently, a large-scale convection cell is established, which resembles porous-medium RBC. Although both velocity and heat transfer are expected to increase with
$R{{a}_{D}}$
, their impact on the distribution of temperature is minimal, leading to distinct vertical isothermal layer where
$T\approx T_{Ra_D=0}$
,
$\textit{Nu}\approx 1$
and
$\theta \approx \pi /2$
.
Regime II (transitional regime, distinct mixed isothermal layer): heat transfer is primarily governed by both thermal conduction and convection, allowing
$\textit{Nu}-1$
to reach the scale of
$O(1)$
. Figure 2(e) suggests that
$|\overline {\boldsymbol{u}}_D|$
begins to concentrate in elongated vertical regions near the sidewalls, indicating the formation of kinetic and thermal boundary layers. Similar to porous-medium RBC, the large-scale structure in the flow field begins to exhibit disorder. Consequently, thermal convective instability gradually develops near the hot and cold walls, impacting heat transfer. Furthermore, as illustrated in figure 2(b), the isothermal layer deviates from the
$x_2$
axis towards the
$x_1$
axis, signifying a distinct mixed isothermal layer with
$0\lt \theta \lt \pi /2$
.
Regime III (turbulent regime, distinct horizontal isothermal layer): heat transfer is predominantly governed by a weak and pore-scale turbulent flow (Rao et al. Reference Rao, Kuznetsov and Jin2020; Rao & Jin Reference Rao and Jin2022; Rao & Zhang Reference Rao and Zhang2025a
,
Reference Rao and Zhangb
,Reference Rao and Zhang
c
). The onset of turbulence in this regime induces spatio-temporal chaos: the flow originating from the bottom near the left wall rises directly to the top, while the flow originating from the top near the right wall descends to the bottom, as illustrated in figure 2(c, f). The vorticity, defined as
$\omega ={\partial {{u}_{2}}}/{\partial {{x}_{1}}}-{\partial {{u}_{1}}}/{\partial {{x}_{2}}}$
, is utilised to qualitatively identify turbulent structures within porous media, as it effectively captures vortex structures of varying strengths, as illustrated in figure 3. The results indicate that an increase in
$Ra_D$
correlates with an intensified vorticity strength, while a decrease in
$\phi$
is associated with a reduction in the high-vorticity region. Turbulence is primarily concentrated near the sidewalls, resulting in a significant increase in the Nusselt number. The temperature distribution leads to the formation of a distinct horizontal isothermal layer with
$\theta \approx 0$
.
(a) Filtered time- and REV-averaged temperature contour with isothermal layer
$\overline {T}=T_0$
. (b) The profile of
$\theta \text{-}R{{a}_{D}}$
at
$\phi =0.69$
,
$Da=2.1\times {{10}^{-5}}$
.
Furthermore, we apply two statistical methods to identify the three regimes. Method 1 involves calculating the volume-averaged temperature within each REV and subsequently mapping all REV-averaged values to a pure cavity of identical size. The angle
$\theta =\arctan ( {\Delta {{x}_{2}}}/{\Delta {{x}_{1}}} )$
between the isothermal layer and the
$x_1$
axis is then calculated, as illustrated in figure 4(a). Method 2 employs
$\theta =\pi /2\ -\arctan ( {{{\left \langle T\boldsymbol{\cdot }{{x}_{2}} \right \rangle }_{f}}}/{{{\left \langle T\boldsymbol{\cdot }{{x}_{1}} \right \rangle }_{f}}})$
to determine the corresponding angle. Figure 4(b) shows the profile of
$\theta -R{{a}_{D}}$
across three distinct regimes. In both regimes I and III,
$\theta$
increases slowly with values close to
$\pi /2$
and
$0$
, respectively. In the transitional regime II,
$\theta$
increases dramatically over a small range of
$R{{a}_{D}}$
, showing that convection begins to significantly influence the temperature distribution. In addition, the profile of
$\theta -R{{a}_{D}}$
is almost rotationally symmetric about the point at the centre of regime II and
$\theta =\pi /4$
. The angle obtained by method 2 is larger than that obtained by method 1, particularly in regime II, as method 2 incorporates the thermal distribution across the entire area, while method 1 only accounts for the thermal distribution in the central region.
3.2. Scaling laws of
$\textit{Nu}$
and
$\textit{Re}$
This section aims to theoretically derive unified scaling laws for
$\textit{Nu}$
and
$\textit{Re}$
in porous-medium vertical natural convection. This analysis encompasses a broad spectrum of input parameters
$(R{{a}_{D}}, \phi , Da, \textit{Pr})$
, as these parameters are crucial for facilitating rapid predictions of heat transfer efficiency and flow intensity. To ensure a smooth derivation, two coarse-grained physical coefficients, the permeability
$K$
and static effective diffusivity
${\tilde {\alpha }}_{0}$
, must be determined in advance. Here,
${\tilde {\alpha }}_{0}$
is defined by the cavity-scale extension of Fourier’s law. In this study, we focus solely on the channels within the porous matrix, as they primarily determine the values of permeability and static effective diffusivity. Based on the approximate analysis of Stokes flow and static conduction, these two terms can be expressed as follows:
\begin{align} K\approx {{H}^{2}}{( 1-\sqrt {1-\phi } )^3}/\big( 12N_{r}^{2} \big), \\[-28pt] \nonumber \end{align}
\begin{align} {{\tilde {\alpha }}_{0}}=-(\alpha {{\left \langle {{\partial }_{1}}{{T}_{Ra=0}} \right \rangle }_{A}}{{|}_{{{x}_{1}}=\pm H/2}}){H}/{\Delta }\approx \alpha \big( 1-\sqrt {1-\phi } \big). \\[-2pt] \nonumber \end{align}
Regarding the presence of the porous matrix, the expressions for the global kinetic and thermal dissipation rates are adjusted, as shown below (see Appendix A)
\begin{align} {{\varepsilon }_{u}}=\nu {{\langle {{\partial }_{i}}{{u}_{\!j}}{{\partial }_{i}}{{u}_{\!j}} \rangle }_{V}} =\beta g{{\left \langle {{u}_{2}}T \right \rangle }_{V}}, \\[-28pt] \nonumber \end{align}
\begin{align} {{\varepsilon }_{T}}=\alpha {{\langle {{\partial }_{i}}T{{\partial }_{i}}T \rangle }_{V}} ={{\tilde {\alpha }}_{0}}{{\Delta }^{2}}{{H}^{-2}}\textit{Nu} ={{\tilde {\alpha }}_{0}}{{\Delta }^{2}}{{H}^{-2}}+\alpha {{\left \langle {{\partial }_{i}}\delta T{{\partial }_{i}}\delta T \right \rangle }_{V}}, \\[-2pt] \nonumber \end{align}
where
${\left \langle \boldsymbol{\cdot }\right \rangle }_{V}$
denotes the time- and whole volume-averaged value including the fluid and solid phases. Here,
$\delta T=T-{{T}_{Ra=0}}$
and
${T}_{Ra=0}$
represents the temperature distribution when
$Ra=0$
.
3.2.1. Steady regime
In this regime, the total drag caused by the porous matrix, as approximated by Darcy’s law, equals the buoyancy term induced by the temperature difference and
$T\approx T_{Ra=0}$
. Consequently, it is straightforward to derive the scaling law of
$\textit{Re}$
\begin{equation} \phi \nu {{K}^{-1}}{{u}_{Di}}\sim -\phi {{g}_{i}}\beta \left ( {{T}_{Ra=0}}-{{T}_{0}} \right ) \,\Rightarrow \, {{u}_{Di}}\sim g_i\beta \Delta K{{\nu }^{-1}} \,\Rightarrow \, \textit{Re}\sim R{{a}_{D}}{\textit{Pr}^{-1}}. \end{equation}
This also indicates that the transitional regime, representing the onset of deviation from Stokes approximation, should be characterised by
$R{{a}_{D}}$
and
$\textit{Pr}$
.
Notably,
$\partial _{i}^{2}{{T}_{Ra=0}}=0$
and
${{\partial }_{1}}{{T}_{Ra=0}}\sim -\Delta {{H}^{-1}} ( {{{\tilde {\alpha }}}_{0}}/\alpha )$
in this regime. Consequently, the scaling of
$\delta T$
can be derived based on steady macroscopic energy equation (Rao & Zhang Reference Rao and Zhang2025a
,
Reference Rao and Zhangb
)
\begin{equation} {{\phi }^{-1}}{{u}_{Di}}{{\partial }_{i}}{{T}_{D,Ra=0}}={{\tilde {\alpha }}_{0}}{{{\nabla} }^{2}}\delta {{T}_{D}} \,\Rightarrow \, \delta T\sim {{\phi }^{-1}}\Delta H{{\alpha }^{-1}}{{u}_{Di}}\sim {{\phi }^{-1}}R{{a}_{D}}\Delta . \end{equation}
The scaling law for
$\textit{Nu}$
is derived by combining (3.2b
) and (3.4), as shown below
\begin{equation} \textit{Nu}-1={\big( {{{\tilde {\alpha }}}_{0}}{{\Delta }^{2}}{{H}^{-2}} \big)^{-1}}\alpha {{\langle {{\partial }_{i}}\delta T{{\partial }_{i}}\delta T \rangle }_{V}}\sim Ra_{D}^{2}{\big( \phi \big( 1-\sqrt {1-\phi } \big) \big)^{-1}}. \end{equation}
3.2.2. Turbulent regime
In this regime, the GL theory and the process of its derivation (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009) are applied with adjustments concerning the porous matrix. The total heat flux can be further elaborated based on the definition of
$\textit{Nu}$
\begin{equation} {{\tilde {\alpha }}_{0}}\Delta {{H}^{-1}}\textit{Nu}=-\alpha {{\langle {{\partial }_{1}}T \rangle }|_{x=\pm H/2}}={{\langle {{u}_{1}}T-\alpha {{\partial }_{1}}T \rangle }_{A}}\approx {{\langle {{u}_{1}}T \rangle }_{V}}. \end{equation}
It can be inferred from figure 2(c, f) that
$\langle u_1T\rangle _V$
is dominated by the convective heat flux along the stratified regions near the horizontal walls, while
$\langle u_2T\rangle _V$
is dominated by the convective heat flux along the boundary layers near vertical walls. With the heat input/output from vertical walls, the vertical heat fluxes caused by
$u_2$
gradually increases along their flow paths, and change directions at the top-left/bottom-right corner, transitioning to the heat flux primarily driven by
$u_1 T$
. Consequently, it is reasonable to estimate that
\begin{equation} (u_1T\delta x_2)_{\rightleftarrows }\approx 2(u_2T\delta x_1)_{\uparrow \downarrow } \,\Rightarrow \, \langle u_1T\rangle _V\sim \langle u_2T\rangle _V, \end{equation}
where
$\rightleftarrows$
refers to the rightward/leftward flow near the upper/lower plate, and
$\uparrow \downarrow$
refers to the upward/downward flow near the left/right plate. Therefore, by combining (3.2a
), (3.6) and (3.7) we can get
\begin{equation} {{\varepsilon }_{u}}\sim \beta g{{\tilde {\alpha }}_{0}}\Delta {{H}^{-1}}\textit{Nu}. \end{equation}
Considering the presence of porous matrix, the no-slip and adiabatic boundary conditions imposed on the porous matrix,
${{\left \langle \boldsymbol{\cdot }\right \rangle }_{V}}=\phi {{\left \langle \boldsymbol{\cdot }\right \rangle }_{f}}$
,
${{\partial }_{i}}{{u}_{\!j}}\sim {U}/{ ( S\times ( 1-\sqrt {1-\phi } ) )}={ ( U\times {{N}_{r}} )}/$
$ ( H\times ( 1-\sqrt {1-\phi } ) )$
,
${{\partial }_{i}}T\sim {\Delta }/{ ( H\times ( 1-\sqrt {1-\phi } ) )}$
. By integrating these with the original definitions of
${\varepsilon }_{u}$
and
${\varepsilon }_{T}$
in (3.2), the global kinetic and thermal dissipation rates are adjusted
\begin{align} {{\varepsilon }_{u}}\sim {{U}^{3}}{{H}^{-1}}\times \phi \times {( H/( S-d ) )^{2}}\sim {{U}^{3}}{{H}^{-1}}N_{r}^{2}\phi {\big( 1-\sqrt {1-\phi } \big)^{-2}}, \\[-28pt] \nonumber \end{align}
\begin{align} {{\varepsilon }_{T}}\sim {{\Delta }^{2}}{{H}^{-1}}U{{\lambda }_{T}}\lambda _u^{-1}\times \phi \times {( S/( S-d ) )^{2}}\sim {{\Delta }^{2}}{{H}^{-1}}U{{\lambda }_{T}}\lambda _{u}^{-1}\phi {\big( 1-\sqrt {1-\phi } \big)^{2}}, \\[-2pt] \nonumber \end{align}
where
$\lambda _u=\textit{HRe}^{-1/2}$
and
$\lambda _T=H(2\tilde {\alpha }_0 \textit{Nu}/\alpha )^{-1}$
. Combining (3.2b
), (3.8) and (3.9), we can finally get
\begin{align} \textit{Nu}\sim Ra_{D}^{1/3}{\big( \phi \big( 1-\sqrt {1-\phi } \big) \big)^{-1}}{\big( \big( 1-\sqrt {1-\phi } \big)/\phi \big)^{-5/3}}, \\[-28pt] \nonumber \end{align}
\begin{align} \textit{Re}\sim Ra_{D}^{4/9}{\textit{Pr}^{-2/3}}{\big( \big( 1-\sqrt {1-\phi } \big)/\phi \big)^{-8/9}}. \\[-2pt] \nonumber \end{align}
Validation of the scaling law exponents with respect to
$\textit{Pr}$
. Here,
$\phi =0.69$
,
$Da=2.1\times {{10}^{-5}}$
. Panels show (a)
$\textit{Nu}-1$
; (b)
$\textit{Re}$
. Panels (a) and (b) share the same legend.
3.3. Validation of scaling laws with respect to
$\textit{Pr}$
In this subsection, we validate the scaling laws of
$\textit{Nu}$
and
$\textit{Re}$
with respect to
$\textit{Pr}$
. Since the derived equations provide the scaling law exponents exclusively for the steady and turbulent regimes, as discussed in § 3.2, we concentrate our validation efforts on these two regimes. Consequently, we maintain fixed values of
$Ra_D$
and
$\phi$
, which would be validated in § 3.4. The Prandtl number is defined as
$\textit{Pr}=2^n$
, where the integer
$n$
varies from
$-3$
to
$3$
. The DNS results are illustrated in figure 5. Our findings indicate that the Nusselt number exhibits approximate independence from the Prandtl number. In the steady regime, the calculated Nusselt number is distinctly identified. Conversely, in the turbulent regime, we observe slight fluctuations, which may be attributed to the size matching between pore-scale and turbulent structures. Furthermore, the results imply that the scaling law exponent for
$\textit{Re}$
is
$-1$
in the steady regime and
$-2/3$
in the turbulent regime. The DNS results effectively validate the scaling laws of
$\textit{Nu}-1$
and
$\textit{Re}$
with respect to
$\textit{Pr}$
derived from theoretical frameworks. It is important to note that this study exclusively considers
$\textit{Pr}$
values within the range of [0.125, 8].
3.4. Unified scaling models and validation for
$\textit{Pr}=1$
This subsection aims to establish unified scaling models for
$\textit{Nu}$
and
$\textit{Re}$
across the three regimes of porous-medium vertical natural convection. The development of these unified scaling models is crucial, as they enable the rapid determination of corresponding
$\textit{Nu}$
and
$\textit{Re}$
values based on prior
$Ra_D$
and
$\phi$
, thereby benefiting industrial applications. In this study,
$\textit{Pr}$
is set to 1, and for the range of porosity considered,
$( 1-\sqrt {1-\phi } )/\phi$
can be approximated as a constant value. Consequently, bridge functions are proposed as illustrated below (Rao & Zhang Reference Rao and Zhang2025a
,
Reference Rao and Zhangb
)
\begin{align} \textit{Nu}-1={{a}_{1}}Ra_D^{1/3}{\big( 1+{\big( R{{a}_{D}}/Ra_{D}^{c1} \big)^{-5/3}} \big)^{-1}}{\big( \phi \big( 1-\sqrt {1-\phi } \big) \big)^{-1}}, \\[-28pt] \nonumber \end{align}
\begin{align} \textit{Re}={{a}_{2}}Ra_{D}^{4/9}{\big( 1+{\big( R{{a}_{D}}/Ra_{D}^{c2} \big)^{-5/9}} \big)^{-1}}, \\[-2pt] \nonumber \end{align}
where the prefactors
$a_1,a_2$
and transitional values
$Ra_{D}^{c1},Ra_{D}^{c2}$
are determined by fitting DNS results for
$\phi =0.69$
and equal
$(0.44,1.57,78.9,258.8)$
respectively. Subsequently, the unified scaling models are validated against additional DNS results across various porosities.
Dependence of (a)
$\textit{Nu}$
and (b)
$\textit{Re}$
on
$Da$
at
$\phi =0.69$
.
$N_r=4,10,20,40$
, corresponding to
$Da={K}/{{{H}^{2}}}=5.3\times {{10}^{-4}}, 8.5\times {{10}^{-5}}, 2.1\times {{10}^{-5}}, 5.3\times {{10}^{-6}}$
, respectively. Panels (a) and (b) share the same legend.
The relationships across a broad range of porosities. Panels show (a)
$\textit{Nu}-1\sim R{{a}_{D}}$
; (b)
$\textit{Re}\sim R{{a}_{D}}$
. Panels (a) and (b) share the same legend.
Our analysis reveals that
$\textit{Nu}$
and
$\textit{Re}$
are independent of
$Da$
by conducting DNS cases, as well as being demonstrated in the unified scaling models. Various
$Da$
values are achieved by altering the number
$N_r$
of porous matrices, while maintaining a constant
$\phi$
, as shown in (2.2) and figure 6. The results indicate that both
$\textit{Nu}$
and
$\textit{Re}$
remain nearly constant at each
$Ra_D$
, thereby confirming the independence of
$\textit{Nu}$
and
$\textit{Re}$
from
$Da$
. However, it is crucial to acknowledge that substantial errors may occur if
$Da$
continues to increase.
The unified scaling models are validated against DNS results, as illustrated in figure 7. The results demonstrate a strong alignment between the model results and the DNS results in both the steady and turbulent regimes across a range of porosities. Notably, the scaling exponents
$(\gamma _1=1/3, \gamma _2=4/9)$
in regime III conform to GL theory (Grossmann & Lohse Reference Grossmann and Lohse2000; Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Liu et al. Reference Liu, Jiang, Chong, Zhu, Wan, Verzicco, Stevens, Lohse and Sun2020), indicating the similarity between the present system and classical/porous-medium RBC. The critical
$R{{a}_{D}}$
value for
$\textit{Nu}$
is approximately 78.9, which closely corresponds with the theoretical critical
$R{{a}_{D}}$
of 98.3 and the experimental value of 66.2 reported in the research (Kwok & Chen Reference Kwok and Chen1987). In contrast, the critical value for
$\textit{Re}$
is approximately 258.8, which exceeds the critical
$R{{a}_{D}}$
value for
$\textit{Nu}$
. It is likely that
$\textit{Nu}$
is significantly influenced by the factor
$R{{a}_{D}}$
that causes it to transition to regime II earlier. Despite the minor discrepancies observed in regime II, the unified scaling models demonstrate universal applicability.
It is worth noting that the unified models presented in (3.11) do not imply that
$Da$
and
$\phi$
have no influence on the scaling laws of
$\textit{Nu}$
and
$\textit{Re}$
. The effect of
$Da$
is partially represented in the term
$R{{a}_{D}}=Ra\boldsymbol{\cdot }Da$
. The influences of
$\phi$
are also reflected in the permeability as described in (3.1a
), and this dependence is embedded in the definition of
$R{{a}_{D}}={g\beta \Delta KH}/{ ( \nu \alpha )}$
. Furthermore, (3.11a
) retains the term
${ ( \phi ( 1-\sqrt {1-\phi } ) )}^{-1}$
, which also affects the determination of
$\textit{Nu}$
. When
$\phi$
approaches 1, the scaling law exponents of
$\textit{Nu}-1\sim Ra_{D}^{\gamma _1 }\textit{Pr}^{\eta _1}$
and
$\textit{Re}\sim Ra_{D}^{\gamma _2 }\textit{Pr}^{\eta _2}$
align with the GL theory and pure vertical natural convection with
$\phi =1$
.
4. Conclusion
We have investigated vertical natural convection in a two-dimensional cavity filled with porous media over a broad range of Rayleigh–Darcy number
$R{{a}_{D}}$
, porosity
$\phi$
and Darcy number
$Da$
. Our findings reveal three major regimes characterised by the angle
$\theta$
of the isothermal layer and flow states. Consequently, we theoretically derive distinct scaling laws for
$\textit{Nu}-1\sim Ra_{D}^{\gamma _1 }\textit{Pr}^{\eta _1}$
and
$\textit{Re}\sim Ra_{D}^{\gamma _2 }\textit{Pr}^{\eta _2}$
. Regime I is dominated by Stokes flow and thermal conduction. We observe
$\theta \approx \pi /2$
and derive
$(\gamma _1,\gamma _2,\eta _1,\eta _2)=(2,1,0,-1)$
based on perturbation analysis. Regime III is characterised by turbulent flow and thermal convection with clear kinetic and thermal boundary layers. We observe
$\theta \approx 0$
and derive
$(\gamma _1,\gamma _2,\eta _1,\eta _2)=(1/3,4/9,0,-2/3)$
, which aligns with GL theory. Regime II represents a transitional flow state, with corresponding parameters that fall between those of regimes I and III. The scaling laws are well validated by DNS results. Furthermore, we have established unified scaling models for
$\textit{Nu}-1\sim Ra_{D}$
and
$\textit{Re}\sim Ra_{D}$
at
$\textit{Pr}=1$
, effectively integrating three major regimes, and these remain accurate across a broad spectrum of
$Ra_D$
and
$\phi$
. Our findings indicate that both
$\textit{Nu}$
and
$\textit{Re}$
are independent of
$Da$
when
$Da$
is sufficiently small. The unified scaling models are instrumental in predicting
$\textit{Nu}$
and
$\textit{Re}$
values based on the governing parameter
$R{{a}_{D}}$
, thereby offering significant benefits to thermal convection and solutal convection.
In this study, we primarily focus on vertical natural convection in regular porous media. Future research will aim to extend this study to irregular porous media, particularly investigating the effects of the random distance between the left/right walls and the nearest porous matrix on the kinetic and thermal boundary layers. Notably, in regime III, the minimum pore scale may influence pore-scale turbulence, potentially altering the corresponding scaling exponents compared with those observed in regular porous media.
Funding
This research was supported by Zhejiang Provincial Natural Science Foundation of China (Grant No. LQN25E060006) and National Natural Science Foundation of China (Grant No. 12302284).
Declaration of interests
The authors report no conflict of interest.
Appendix A. The global kinetic and thermal dissipation rates in porous-medium vertical natural convection
The global kinetic dissipation rate is equal to the global buoyancy work rate. The thermal dissipation rate can be expressed using
$\textit{Nu}$
in a manner similar to classical RBC, but with
$\tilde {\alpha }_0$
instead of
$\alpha$
. Therefore,
${\varepsilon }_{u}$
and
${\varepsilon }_{T}$
can be expressed as
\begin{align} {{\varepsilon }_{u}} & = \nu {{\langle \boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{\nabla }\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u} )-\boldsymbol{u}\boldsymbol{\cdot }{{{\nabla} }^{2}}\boldsymbol{u} \rangle }_{V}} =-{{\langle \boldsymbol{u}\boldsymbol{\cdot }( \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}+\boldsymbol{g}\beta ( T-{{T}_{0}} )+{{\rho }^{-1}}\boldsymbol{\nabla }p ) \rangle }_{V}} \nonumber \\ & = -{{\langle \boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{u}{{|\boldsymbol{u}|}^{2}}/2 )-{{u}_{2}}g\beta T+{{\rho }^{-1}}\boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{u}p ) \rangle }_{V}} =g\beta {{\langle {{u}_{2}}T \rangle }_{V}}, \\[-9pt] \nonumber \end{align}
\begin{align} {{\varepsilon }_{T}} & = \alpha {{\langle \boldsymbol{\nabla }\boldsymbol{\cdot }( T\boldsymbol{\nabla }T )-T{{{\nabla} }^{2}}T \rangle }_{V}} =\alpha {{\langle \boldsymbol{\nabla }\boldsymbol{\cdot }( T\boldsymbol{\nabla }T )-{{( 2\alpha )}^{-1}}\boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{u}{{T}^{2}} ) \rangle }_{V}} \nonumber \\ & = \alpha {{H}^{-2}}\textstyle {\oint _{\partial V_{\!f}}}T\boldsymbol{\nabla }T\boldsymbol{\cdot }\boldsymbol{n}{\rm d}S ={{{\tilde {\alpha }}}_{0}}{{\Delta }^{2}}{{H}^{-2}}\textit{Nu}, \\[9pt] \nonumber \end{align}
where
${V}_{f}$
denotes fluid phase and
$\boldsymbol{n}$
denotes unit vector pointing from the fluid to solid at their interfaces. Regarding
${{{\nabla} }^{2}}{{T}_{Ra=0}}=0$
and
$\delta T{{|}_{{{x}_{1}}=\pm H/2}}=0$
,
${\varepsilon }_{T}$
can be further derived
\begin{equation} \begin{aligned} {{\varepsilon }_{T}} & = \alpha {{\langle {{( \boldsymbol{\nabla }{{T}_{Ra=0}} )}^{2}}+{{( \boldsymbol{\nabla }\delta T )}^{2}}+2\boldsymbol{\nabla }{{T}_{Ra=0}}\boldsymbol{\cdot }\boldsymbol{\nabla }\delta T \rangle }_{V}} \\& ={{{\tilde {\alpha }}}_{0}}{{\Delta }^{2}}{{H}^{-2}}N{{u}_{Ra=0}}+\alpha {{\langle {{( \boldsymbol{\nabla }\delta T )}^{2}} \rangle }_{V}}+2\alpha {{\langle \boldsymbol{\nabla }\boldsymbol{\cdot }( \delta T\boldsymbol{\nabla }{{T}_{Ra=0}} )-\delta T{{{\nabla} }^{2}}{{T}_{Ra=0}} \rangle }_{V}} \\ & ={{{\tilde {\alpha }}}_{0}}{{\Delta }^{2}}{{H}^{-2}}+\alpha {{\langle {{( \boldsymbol{\nabla }\delta T )}^{2}} \rangle }_{V}}+2\alpha {{H}^{-2}}\textstyle {\oint _{\partial V_{\!f}}}{\delta T\boldsymbol{\nabla }{{T}_{Ra=0}}\boldsymbol{\cdot }{\rm d}\boldsymbol{S}} \\ & ={{{\tilde {\alpha }}}_{0}}{{\Delta }^{2}}{{H}^{-2}}+\alpha {{\langle {{( \boldsymbol{\nabla }\delta T )}^{2}} \rangle }_{V}}. \end{aligned} \end{equation}
