2.1. Theoretical framework

Based on the theorem of local volumetric average, the macroscopic equations governing flows in porous media were often derived by applying both volume and time averaging to the Navier–Stokes equations (de Lemos Reference de Lemos2012a ,Reference de Lemos b ). In order to enable the capability of real-time prediction, our approach employs only volume averaging, as demonstrated in (2.1) below. We assume that the fluid is incompressible, and include a heat source term, leaving the density constant despite variations in temperature. The solid phase is static and has a temperature field ${T}_{s}$ whose evolution usually depends on extra solid properties. However, it can be found that the complexities of heat transfer within porous media primarily reside in the fluid phase, either explicitly or implicitly. Therefore, this study focuses on thermal convection within the fluid phase, and only displays the solid-phase energy equation formally. Nevertheless, the modelling framework in this study is applicable to a wide range of geometries and coarse-grained thermodynamic properties of the solid phase. This most crucially requires the incorporation of heat exchange between solid and fluid phases, which involves the local solid–fluid temperature difference. The equations are

(2.1a) \begin{align}& \partial _i\big(\phi {{\langle {{u}_{i}} \rangle }^{f}}\big)=0, \end{align}

(2.1b) \begin{align}& \partial _t\big(\phi {{\langle {{u}_{i}} \rangle }^{f}}\big) +\partial _j\big(\phi {{\langle {{u}_{i}} \rangle }^{f}}{{\langle {{u}_{\!j}} \rangle }^{f}}\big) =-\partial _i\big(\phi {{\langle p \rangle }^{f}}\big)/\rho +\phi {{g}_{i}}-\phi {{\hat {R}}_{i}} +\partial _j{{{\tilde {\tau }}}_{\textit{ij}}}, \end{align}

(2.1c) \begin{align}&\qquad\qquad\quad \partial _t\big(\phi {{\langle {{T}_{f}} \rangle }^{f}}\big) +\partial _i\big(\phi {{\langle {{u}_{i}} \rangle }^{f}}{{\langle {{T}_{f}} \rangle }^{f}}\big) =\phi \dot {q}-\phi \hat {H} +\partial _i {\tilde {J}_{i}}, \end{align}

(2.1d) \begin{align}&\qquad\qquad\qquad\qquad \partial _t\big(\phi \eta _c{{\langle {{T}_{s}} \rangle }^{s}}\big) =\phi \eta _c\dot {q}^s+\phi \hat {H} +\partial _i {{{\tilde {J}}}^s_{i}}, \end{align}

where porosity $\phi$ represents the fraction of void space to total volume in porous media. Also, ${\left \langle \varphi \right \rangle }^{f}$ and ${\left \langle \varphi \right \rangle }^{s}$ denote the intrinsic average values of any fluid and solid field within each REV. The spatial deviation ${}^{f}\varphi =\varphi -{{\left \langle \varphi \right \rangle }^{f}}$ is defined as the difference between the microscopic value and the intrinsic volume-averaged value. In (2.1b ), ${g}_{i}$ is an applied constant body force, ${\hat {R}}_{i}$ is the Darcy–Forchheimer drag force, ${\tilde {\tau }}_{\textit{ij}}$ is a second-order pseudo-stress tensor, $T_f$ is the temperature for the fluid phase, $\dot {q}$ is a heat source term, $\hat {H}$ is the interphase energy resistance, and ${\tilde {J}}_{i}$ is pseudo-flux vector. Terms in the formal solid energy equation (2.1d ), including the REV solid–fluid heat capacity ratio $\eta _c=(1-\phi )\rho ^s c_{\!p}^{s}/(\phi \rho c_p)$ (volume ratio times the volumetric isobaric heat capacity ratio), the solid heat source term $\dot {q}^s$ , and the solid-phase pseudo-flux vector ${\tilde {J}}^s_{i}$ , can be arbitrary and are thus not specified in the present work. The dimensions of the variables include length, time, mass and temperature, which are not explicitly defined in this context.

Figure 1.

A schematic of a transversely isotropic porous matrix made up of a cuboid, with ratio 2.0, and $\phi =0.81$ .

Far away from wall boundaries, the total drag ${R}_{i}$ is decomposed into two components: the surface integrals of the intrinsic average pressure and viscous forces ${R}_{i}^{\gamma }$ , and the surface integrals of the spatial deviations of pressure and viscous forces $\tilde {R}_{i}$ . Here, ${R}_{i}^{\gamma }$ can be approximated to zeroth order of corresponding values multiplied by the negative porosity gradient. Given that the porosity is considered a constant value for macroscopic simulation, ${R}_{i}^{\gamma }$ equals zero. Consequently, the total drag ${{R}_{i}}={{\tilde {R}}_{i}}$ , and is approximated by the Darcy–Forchheimer equation ${\hat {R}}_{i}$ , as shown below (Whitaker Reference Whitaker1986; Breugem, Boersma & Uittenbogaard Reference Breugem, Boersma and Uittenbogaard2006):

(2.2) \begin{equation} {{R}_{i}}=(\phi\, \Delta V)^{-1} {\int _{A}}({{{\hat {n}}_{i}}p/\rho }-{{\nu \hat {n}_j\,\partial _j}{{u}_{i}}})\,{\rm d}S\approx {{\hat {R}}_{i}}=\nu K^{-1}\left ( \sum _{n=0}^{\infty }{{{c}_{Fn}}\,Re_{K}^{n}} \right ){{u}_{Di}}, \end{equation}

which adopts scalar permeability $K$ and Taylor expansion coefficients ${c}_{Fn}$ to account for the total drag exclusively caused by isotropic porous media, with ${c}_{F0}$ being 1. Here, $\Delta V$ denotes the volume of one REV, $A$ denotes the fluid–solid interfacial area, and $\hat {n}_i$ (in surface integrals) denotes the normal direction from the fluid to porous matrix at the solid–fluid interface. To extend this model for anisotropic porous media, we should generalise $K$ and ${c}_{Fn}$ into tensors that reflect anisotropic dependence of $\hat {R}_i$ on flow direction and local Reynolds number ${\textit{Re}_{K}}={| {{u}_{D}} |\sqrt {| {{K}_{\textit{ij}}} |}}/{\nu }$ , where ${{u}_{Di}}=\phi {{\left \langle {{u}_{i}} \right \rangle }^{f}}$ :

(2.3) \begin{equation} {{\hat {R}}_{i}}=\nu K_{ik}^{-1}\left ( \sum_{n=0}^{\infty }{{{c}_{Fnkj}}\,\textit{Re}_{K}^{n}} \right ){{u}_{Dj}}. \end{equation}

Here, ${c}_{Fnkj}$ represent the coefficients of the Taylor series, with ${c}_{F0kj}$ being a unit matrix. For simplicity, in the present work we consider the porous matrices that are isotropic in the planes perpendicular to a specific direction $\tilde {\boldsymbol{n}}$ (transversely isotropic, as sketched in figure 1), which represents a wide range of anisotropic porous media. With additional symmetry required by energy conservation, ${K}_{\textit{ij}}$ and ${c}_{Fnkj}$ can be decomposed as follows:

(2.4) \begin{equation} {{K}_{\textit{ij}}}={{K}_{L}}{{\tilde {n}}_{i}}{{\tilde {n}}_{\!j}}+{{K}_{T}}( {{\delta }_{\textit{ij}}}-{{\tilde {n}}_{i}}{{\tilde {n}}_{\!j}}),\quad{{c}_{Fnij}}={{c}_{FnL}}{{\tilde {n}}_{i}}{{\tilde {n}}_{\!j}}+{{c}_{FnT}}( {{\delta }_{\textit{ij}}}-{{\tilde {n}}_{i}}{{\tilde {n}}_{\!j}}). \end{equation}

The characteristic direction $\tilde {{n}}_{i}$ is associated with the detailed geometry of the porous matrix, and ${\delta }_{\textit{ij}}$ denotes unit matrix. The permeability ${K}_{\textit{ij}}$ indicates the capacity of the porous medium to facilitate fluid flow, and $| {{K}_{\textit{ij}}} |=[3/(K_{L}^{-2}+K_{T}^{-2}+K_{T}^{-2})]^{1/2}$ . Dimensionless coefficients ${c}_{Fnij}$ account for the nonlinear increase in form drag with velocity. In this study, $L$ denotes the longitudinal direction parallel to the characteristic direction $\tilde {n}_i$ , while $T$ represents any transverse direction perpendicular to $\tilde {n}_i$ . It is crucial to clarify that $L$ and $T$ are used to indicate the orientations of the anisotropic porous matrix, as illustrated on the right-hand side of figure 1. Furthermore, coefficient values corresponding to all transverse directions are identical, as required by the transversely isotropic geometry of the porous matrix.

In regions characterised by non-negligible macroscopic velocity gradient, the drag deviation $\phi ({{R}_{i}}-{{\hat {R}}_{i}})$ , the turbulent momentum dispersion $\phi {{\langle {}^{f}{{u}_{i}}{}^{f}{{u}_{\!j}}\rangle }^{f}}$ and molecular diffusion $2\nu {{S}_{Dij}}$ ( ${{S}_{Dij}}= ( {{\partial }_{i}}{{u}_{Dj}}+{{\partial }_{\!j}}{{u}_{Di}} )/2$ ) are significantly influenced by the gradient of the macroscopic velocity. To account for anisotropic effects, a second-order pseudo-stress tensor ${\tilde {\tau }}_{\textit{ij}}$ should be introduced instead as follows:

(2.5) \begin{equation} \phi \big({{R}_{i}}-{{\hat {R}}_{i}}\big) +\nu\, \partial _{jj} \big(\phi {{\langle {{u}_{i}} \rangle }^{f}}\big) -\partial _j \big(\phi {{\langle {}^{f}{{u}_{i}}{}^{f}{{u}_{\!j}} \rangle }^{f}}\big) =\partial _j {{{\tilde {\tau }}}_{\textit{ij}}}, \end{equation}

where ${\tilde {\tau }}_{\textit{ij}}$ must be symmetric due to angular momentum conservation, and can be made traceless due to incompressibility. The total drag ${R}_{i}$ is modelled solely based on the assumption of homogeneity. This assumption may introduce errors, which can be conveniently compensated by adjusting the pseudo-stress tensor. Consequently, ${{R}_{i}}-{{\hat {R}}_{i}}$ is included in the definition of ${\tilde {\tau }}_{\textit{ij}}$ in (2.5). Linear mappings from ${S}_{Dij}$ to ${\tilde {\tau }}_{\textit{ij}}$ must be written as

(2.6) \begin{equation} {{\tilde {\tau }}_{\textit{ij}}}={{\nu }_{ijkl}}{{S}_{Dkl}}, \end{equation}

with the effective viscosity ${\nu }_{ijkl}$ being a fourth-order tensor with 81 components. The degrees of freedom of ${\nu }_{ijkl}$ can be reduced to $25$ using the symmetric and traceless properties of ${S}_{Dij}$ and ${\tilde {\tau }}_{\textit{ij}}$ , and further reduced to $3$ using the transverse isotropy of porous media. In addition, for the stability of simulation and calibration, the allowed ${{\tilde {n}}_{i}}{{\tilde {n}}_{\!j}}{{\tilde {n}}_{k}}{{\tilde {n}}_{l}}$ term with an excessive degree of anisotropy is temporarily neglected, which will be justified by modelling accuracy. With further arrangements, the formal expression of ${\tilde {\tau }}_{\textit{ij}}$ is

(2.7a) \begin{align}&\qquad\qquad\qquad\qquad\qquad {{\tilde {\tau }}_{\textit{ij}}}=2{{\nu }_{L}}S_{Dij}^{||}+2{{\nu }_{T}}S_{Dij}^{\bot }, \end{align}

(2.7b) \begin{align}& S_{Dij}^{\parallel }=( {{\tilde {n}}_{i}}{{S}_{Djk}}{{\tilde {n}}_{k}}+{{\tilde {n}}_{\!j}}{{S}_{Dik}}{{\tilde {n}}_{k}} )/2-\left ( {{\tilde {n}}_{k}}{{S}_{Dkl}}{{\tilde {n}}_{l}} \right ){{\delta }_{\textit{ij}}}/3,\text{ }S_{Dij}^{\bot }={{S}_{Dij}}-S_{Dij}^{\parallel }. \end{align}

Here, ${\nu }_{L}$ and ${\nu }_{T}$ are respectively the effective longitudinal and transversal viscosities. According to the pore-scale prevalence hypothesis (Jin et al. Reference Jin, Uth, Kuznetsov and Herwig2015), the characteristic length scale of flow in porous media is defined by the pore size $s$ . In this study, we utilise $\sqrt {| {{K}_{\textit{ij}}} |}$ to characterise the flow, as it exhibits a linear relationship with $s$ for a fixed shape of the porous matrix. A Taylor expansion can be expressed with respect to ${\textit{Re}_{d}}={| {{S}_{Dik}}{{K}_{kj}} |}/{\nu }$ , as

(2.8) \begin{equation} {{{\nu }_{L}}}/{\nu }=\sum_{n=0}^{\infty }{{{c}_{LBn}}\,\textit{Re}_{d}^{n}},\quad {{{\nu }_{T}}}/{\nu }=\sum_{n=0}^{\infty }{{{c}_{\textit{TBn}}}\,\textit{Re}_{d}^{n}}. \end{equation}

The interphase energy resistance $H$ should be modelled using the effective interfacial convective heat transfer coefficient $h$ , which primarily accounts for the heat transfer ${\int _{A}{\alpha {{\hat {n}}_{i}}{{\partial }_{i}}{{T}_{f}}\,{\rm d}S}}/{(\phi\, \Delta V})$ between the fluid and solid phases, expressed as

(2.9) \begin{equation} H=-(\phi \Delta V)^{-1}{\int _{A}}{\alpha {{\hat {n}}_{i}}\,{{\partial }_{i}}{{T}_{f}}\,{\rm d}S}\approx \hat {H}=ha\big( {{\langle {{T}_{f}} \rangle }^{f}}-{{\langle {{T}_{s}} \rangle }^{s}} \big), \end{equation}

where $\alpha$ denotes molecular thermal diffusivity, and $a={A}/{\Delta V}$ represents the ratio of the porous matrix surface area to the volume within a REV. From DNS results, two distinct $h\sim \textit{Re}_K^\beta$ scalings are newly identified: $\beta =2$ in the Darcy regime ( $\textit{Re}_K\lesssim 1$ ), and $\beta =0.5$ in the Forchheimer regime ( $\textit{Re}_K\gtrsim 1$ ). Since the larger/smaller scaling occurs in the low/high ${\textit{Re}_{K}}$ regime correspondingly, simple addition like the Taylor expansion in (2.3) would incorrectly exchange scalings in these regimes, thus a bridge function should be introduced instead:

(2.10) \begin{equation} {h\,\sqrt {| {{K}_{\textit{ij}}} |}}/{\alpha }={{c}_{h0}}+{{c}_{h1}\, \textit{Re}_{K}}^{{1}/{2}}{{\big( 1+{{\left ( {\textit{Re}_{K}}/{{{c}_{h2}}} \right )}^{-{3}/{2}}} \big)}^{-1}}. \end{equation}

The angle between fluid velocity and $\tilde {\boldsymbol{n}}$ will affect heat transfer and must be considered:

(2.11a) \begin{align}&\,\, {{c}_{h0}}={{{h}_{0}}\,\sqrt {| {{K}_{\textit{ij}}} |}}/{\alpha },\quad {{c}_{h1}}={{c}_{h1ij}}{{{u}_{Di}}{{u}_{Dj}}}/{{{\left | {{u}_{D}} \right |}^{2}}},\quad {{c}_{h2}}={{c}_{h2ij}}{{{u}_{Di}}{{u}_{Dj}}}/{{{\left | {{u}_{D}} \right |}^{2}}}, \end{align}

(2.11b) \begin{align}& {{c}_{h1ij}}={{c}_{h1L}}{{\tilde {n}}_{i}}{{\tilde {n}}_{\!j}}+{{c}_{h1T}}\left ( {{\delta }_{\textit{ij}}}-{{\tilde {n}}_{i}}{{\tilde {n}}_{\!j}} \right )\!,\quad {{c}_{h2ij}}={{c}_{h2L}}{{\tilde {n}}_{i}}{{\tilde {n}}_{\!j}}+{{c}_{h2T}}\left ( {{\delta }_{\textit{ij}}}-{{\tilde {n}}_{i}}{{\tilde {n}}_{\!j}} \right )\!, \end{align}

where ${h}_{0}$ denotes the heat transfer coefficient in the absence of thermal convection. In this study, the Prandtl number $\textit{Pr}={\nu }/{\alpha }$ is set to 1 for simplicity and is therefore omitted.

Similar to the modelling of momentum dispersion, interphase energy resistance $\phi ( H-\hat {H} )$ , molecular thermal diffusion ${\alpha\, \partial \phi {{\langle {{T}_{f}} \rangle }^{f}}}/{\partial {{x}_{i}}}$ and thermal dispersion $\phi {{\langle {}^{f}{{u}_{i}}{}^{f}{{T}_{f}} \rangle }^{f}}$ are accounted for by the divergence of a pseudo-flux vector ${\tilde {J}}_{i}$ of heat, rather than the commonly employed Laplacian of temperature, to provide a more precise depiction of heterogeneity:

(2.12) \begin{equation} \phi \big( H-\hat {H} \big) +\alpha\, \partial _{ii} \big(\phi {{\langle {{T}_{f}} \rangle }^{f}}\big) -\partial _i \big(\phi {{\langle {}^{f}{{u}_{i}}{}^{f}{{T}_{f}} \rangle }^{f}}\big) =\partial _i {{{\tilde {J}}}_{i}}. \end{equation}

To account for anisotropy, we introduce a second-order anisotropic diffusivity tensor ${\tilde {\alpha }}_{\textit{ij}}$ to map the temperature gradient to ${\tilde {J}}_{i}$ :

(2.13a) \begin{align}& {{\tilde {J}}_{i}}={{\tilde {\alpha }}_{\textit{ij}}}\,\partial _j \big(\phi {{\langle {{T}_{f}} \rangle }^{f}}\big),\quad {{\tilde {\alpha }}_{\textit{ij}}}={{\alpha }_{L}}{{\tilde {n}}_{i}}{{\tilde {n}}_{\!j}}+{{\alpha }_{T}}\left ( {{\delta }_{\textit{ij}}}-{{\tilde {n}}_{i}}{{\tilde {n}}_{\!j}} \right )\!, \end{align}

(2.13b) \begin{align}&\qquad\,\,\, {{{\alpha }_{L}}}/{\alpha }=\sum_{n=0}^{\infty }{{{c}_{\textit{Ln}}}\,\textit{Re}_{d}^{n}},\quad {{{\alpha }_{T}}}/{\alpha }=\sum_{n=0}^{\infty }{{{c}_{\textit{Tn}}}\,{\textit{Re}}_{d}^{n}}, \end{align}

where ${\alpha }_{L}$ and ${\alpha }_{T}$ represent the longitudinal and transversal thermal diffusivities.

As shown above, in the closure of macroscopic equations for anisotropic porous media, the tensor approach is used intensively to account for heterogeneity and anisotropy, and the sufficiency and necessity will be shown in the following subsections. The corresponding Taylor expansions are determined based on the specific geometry of porous matrix, as elaborated in § 2.2. For Reynolds number ranges not too large, two leading terms (i.e. the constant and linear terms) of each Taylor expansion can be sufficient.

2.2. Calibration strategy

Without loss of key mechanisms, in all calibration and validation cases we fix the temperature ${T}_{s}$ of the porous matrix at zero for simplicity, which can still sufficiently capture the laws of non-equilibrium heat transfer between fluid and solid. This allows the solid phase to function as either a heat sink or a heat source when the surrounding fluid temperature deviates from zero. This approach is analogous to the commonly employed no-slip condition where the solid exhibits velocity zero, thereby acting as a momentum sink or source.

The model coefficients ${\hat {R}}_{i}$ , ${\tilde {\tau }}_{\textit{ij}}$ , $\hat {H}$ and ${\tilde {J}}_{i}$ are geometric parameters that are independent of the flow conditions and must be determined prior to solving the macroscopic equations. Only the two leading terms of the corresponding Taylor expansions are adopted, which were previously sufficient for modelling purposes (Rao et al. Reference Rao, Kuznetsov and Jin2020; Rao & Zhang Reference Rao and Zhang2025b ). To determine these coefficients, a schematic geometry of the porous matrix with ratio ${{{d}_{1}}}/{{{d}_{2}}}=2.0$ and $\phi =0.81$ is utilised, as illustrated in figure 1. The characteristic direction related to the detailed geometry of the porous matrix is $\tilde {n}=(1\ 0\ 0)$ . This porous matrix consists of 256 REVs arranged in an $8s\times 8s\times 4s$ configuration, oriented both parallel and perpendicular to the flow direction. To derive a complete set of model coefficients characterising anisotropy, different types of channel flows are necessary, with a unit vector ${n}_{w}$ denoting the wall-normal direction, and a driving force $g$ perpendicular to ${n}_{w}$ . Periodic boundary conditions are applied in the other two directions parallel to the channel walls. The heat source is fixed at $\dot {q}=1$ , and the temperature of the channel walls and the porous matrix is fixed at ${{T}_{s}}=0$ . The Reynolds number is varied by adjusting ${g}_{i}$ . The configuration is adequate to facilitate the development of turbulent flow and ensure the accurate calculation of the model coefficients.

The finite volume method (FVM) is employed to address the Navier–Stokes equations, while the FVM and the spectral method (SM) are specifically utilised for solving the macroscopic equations. The FVM and SM solvers have already been extensively validated and verified in prior studies (Rao et al. Reference Rao, Kuznetsov and Jin2020; Rao Reference Rao2022; Rao & Jin Reference Rao and Jin2022; Rao & Zhang Reference Rao and Zhang2025a , Reference Rao and Zhangb ), demonstrating the reliability and accuracy of the solvers.

In the core region of the fully developed channel flow, the time-averaged velocity and temperature are strongly homogeneous, and (2.1b , c ) can be simplified to

(2.14) \begin{equation} {{\hat {R}}_{i}}={{g}_{i}},\quad\hat {H}=\dot {q}. \end{equation}

To determine the coefficients ${\hat {R}}_{i}$ and $\hat {H}$ , we first specify ${g}_{i}$ and subsequently calculate the macroscopic $\phi {{\left \langle {{u}_{i}} \right \rangle }^{f}}$ and $\phi {{\langle {{T}_{f}} \rangle }^{f}}$ using DNS results. We then utilise (2.14) to obtain the values of ${\hat {R}}_{i}$ and $\hat {H}$ . A series of DNS cases is conducted, allowing the model coefficients in (2.3) and (2.10) to be obtained directly by fitting DNS results, as illustrated in figure 2. It can be implied from figure 2 that the total drag ${\hat {R}}_{i}$ and the interphase energy resistance $\hat {H}$ exhibit anisotropic properties, with the overall drag and energy resistance coefficients in the transverse direction being larger than those in the longitudinal direction. Regardless of the direction of ${g}_{i}$ , $h$ approaches ${h}_{0}$ as ${\textit{Re}_{K}}\to 0$ , as depicted in figure 2(b).

Figure 2.

The determination of coefficients by fitting DNS results, ${{g}_{2}}-{{n}_{w}}=(1\ 0\ 0)$ (black, squares, left-hand axis) and ${{g}_{1}}-{{n}_{w}}=(0\ 1\ 0)$ (red, circles, right-hand axis), for (a) ${\hat {R}}_{i}$ , (b) $\hat {H}$ , ${{c}_{h0}}=1.6$ .

In the first REV close to the wall, which is the only region where momentum and thermal dispersions are non-negligible, the coupled macroscopic equations can be simplified as

(2.15) \begin{equation} \partial _j {{{\tilde {\tau }}}_{\textit{ij}}}=\phi \big( {{{\hat {R}}}_{i}}-{{g}_{i}} \big),\quad \partial _i {{{\tilde {J}}}_{i}}=\phi \big( \hat {H}-\dot {q} \big). \end{equation}

The errors associated with ${\hat {R}}_{i}$ and $\hat {H}$ in the central region may affect the accuracy of the calculated ${\tilde {\tau }}_{\textit{ij}}$ and ${\tilde {J}}_{i}$ . To reduce this uncertainty, the coefficient values of ${\hat {R}}_{i}$ and $\hat {H}$ are fine-tuned to ensure that the computed velocity $\phi {{\left \langle {{u}_{i}} \right \rangle }^{f}}$ and temperature $\phi {{\langle {{T}_{f}} \rangle }^{f}}$ profiles at the centreline align with the DNS results. In line with our previous research (Rao et al. Reference Rao, Kuznetsov and Jin2020; Rao & Zhang Reference Rao and Zhang2025b ), the corresponding effective viscosity and thermal diffusivity are determined by fitting DNS results for both low and high ${\textit{Re}_{K}}$ , as shown in figure 3(a,d). And the transition between the two scaling regimes occurs at approximately ${\textit{Re}_{K}}=1$ . The negative effective viscosity shown in figure 3(a) may infer a positive contribution to macroscopic kinetic energy, and is already explained in flows without porous media (Bayly & Yakhot Reference Bayly and Yakhot1986; Bhaduri Reference Bhaduri2021). The negative ${\nu }_{L}$ in the present work may result from an overestimation of total drag, serving to compensate for this error.

3.1. Flow in a channel occupied by porous media with one length scale

The macroscopic momentum equation is integrated with the macroscopic energy equation to simulate turbulent forced thermal convection within a channel filled with either anisotropic or isotropic porous media. The schematic arrangement of a porous matrix in figure 1 is utilised. Simulations are conducted for both laminar and turbulent flow regimes. Here, $\phi$ is fixed at $0.81$ , with the ratio fixed at 2 or 0.5 representing an anisotropic porous matrix, or ratio being 1 representing the isotropic case. The coefficients ${\hat {R}}_{i}$ , ${\tilde {\tau }}_{\textit{ij}}$ , $\hat {H}$ and ${\tilde {J}}_{i}$ are determined by fitting DNS results using the methods outlined in § 2.2. Subsequently, the corresponding expressions for ${\hat {R}}_{i}$ , ${\tilde {\tau }}_{\textit{ij}}$ , $\hat {H}$ and ${\tilde {J}}_{i}$ are directly applied in (2.1) to validate the macroscopic models against additional DNS cases. The pseudo-stress and pseudo-flux significantly influence the velocity and temperature distribution only within the first REV adjacent to the wall boundary. Consequently, the model results are compared with DNS results for the volume-averaged velocity and temperature in the first REV next to the wall, as shown in figure 3.

Figure 3.

The volume-averaged velocity and temperature in the first REV adjacent to the wall between DNS results (hollow) and model results (solid), for $\phi =0.81$ and $\textit{Pr}=1.0$ . Here, ${\langle {{u}_{D2}} \rangle }^{V}$ is from ${{g}_{2}}-{{n}_{w}}=(1\ 0\ 0)$ cases (black, squares), and ${\langle {{u}_{D3}} \rangle }^{V}$ is from ${{g}_{3}}-{{n}_{w}}=(0\ 1\ 0)$ cases (red, circles). Plots of (a–c) velocity, and (d–f) temperature, with (a,d) ratio 2.0, (b,e) ratio 1.0, (c, f) ratio 0.5.

The model results remain consistent with those from DNS, verifying the sufficiency of the present modelling framework and the truncated Taylor expansions for channel flow in the present range of ${\textit{Re}_{s}}$ . As ${\textit{Re}_{s}}$ increases, the velocity increases, leading to the decrease of temperature. More importantly, when the ratio is larger, equal to and smaller than 1, the velocity driven by the longitudinal force is larger, equal to and smaller than that driven by the transverse force of the same amplitude, leading to the opposite trend of temperature.

Since the anisotropic properties of momentum and thermal dispersions are intricately linked to the specific characteristics of the porous matrix, it is essential to directly compare the fitted dispersion coefficients across various ratio values. In isotropic porous media with ratio 1 (figure 3 b), the transverse components of both momentum and thermal dispersion coefficients ( ${\nu }_{T}$ and ${\alpha }_{T}$ ) are equal to the corresponding longitudinal components ( ${\nu }_{L}$ and ${\alpha }_{L}$ ), and both values are positive, indicating balanced dispersion capability in all directions. In anisotropic porous media with ratio 2 (figure 3 a), transverse dispersion coefficients exceed longitudinal values and the values in isotropic porous media of the same porosity. Such anisotropy reflects the enhanced capacity to transport momentum and heat in the transverse direction according to (2.7), (2.8) and (2.13). In contrast, in anisotropic porous media with ratio 0.5 (figure 3 c), transverse dispersion coefficients are less than the longitudinal values. Most importantly, the differences indicate the necessity of anisotropic parameters governing momentum and thermal dispersions.

3.2. Flow in a channel occupied by porous media with two length scales

This subsection aims to validate the coupled macroscopic equations in porous media characterised by two distinct length scales. In the DNS, the porous matrices consist of cuboids and bars featuring a large element size ${{d}_{l}}=8s$ and infinite length, as illustrated in figure 4(a). In contrast, the macroscopic simulation utilises porous matrices composed solely of square bars as shown in figure 4(b). The entire simulation domain is configured as $32s\times 32s\times 16s$ for the cuboids with ratio 2, ensuring that turbulent flow can develop fully without being constrained by the computational domain. Periodic boundary conditions are applied in all directions, with the body force ${g}_{i}$ imposed in the streamwise ( ${x}_{1}$ ) direction.

Figure 4.

(a,c, e) The DNS, and (b,d, f) the macroscopic simulation, for ${\textit{Re}_{s}}=793$ and $\textit{Pr}=1.0$ : (a,b) the schematic of the porous matrix with two length scales; (c,d) the instantaneous velocity magnitude; (e, f) the instantaneous temperature.

Figure 5.

Distributions of ${\langle {{u}_{D1}} \rangle }^{V,{{x}_{3}}}$ , ${\langle {{u}_{D2}} \rangle }^{V,{{x}_{3}}}$ and ${\langle T \rangle }^{V,{{x}_{3}}}$ along the streamwise ( ${x}_{1}$ ) lines, for ${\textit{Re}_{s}}=793$ and $\textit{Pr}=1.0$ . Hollow symbols with lines represent DNS results, while solid symbols indicate model results. The three lines are represented by circles (red), squares (black), and triangles (blue) in the sequence, ${{{x}_{2}}}/{s}=4.5$ , 16.5 and 27.5: (a) ${\langle {{u}_{D1}} \rangle }^{V,{{x}_{3}}}$ , (b) ${\langle {{u}_{D2}} \rangle }^{V,{{x}_{3}}}$ , (c) ${\langle T \rangle }^{V,{{x}_{3}}}$ .

The instantaneous velocity and temperature results for both DNS and the macroscopic model are presented in figure 4(c–f). The DNS results show no evidence of large-scale turbulence. Similarly, the model results do not exhibit any turbulence as the turbulence confined to the pore scale has already been filtered through the macroscopic models, which aligns with prior research (Rao & Jin Reference Rao and Jin2022), emphasising that macroscopic models based on the pore-scale prevalence hypothesis effectively account for pore-scale turbulence effects without displaying large-scale turbulence. Figure 5 illustrates the time, REV and spanwise ( ${x}_{3}$ ) averaged components ${\left \langle {{u}_{D1}} \right \rangle }^{V,{{x}_{3}}}$ , ${\left \langle {{u}_{D2}} \right \rangle }^{V,{{x}_{3}}}$ and ${\left \langle T \right \rangle }^{V,{{x}_{3}}}$ along the streamwise ( ${x}_{1}$ ) direction. The three selected streamwise lines are positioned farthest from the bars and intersect the centres of the bars, providing a robust basis for comparing DNS and model results. These quantitative findings indicate a satisfactory accuracy of our macroscopic models. The minor differences observed in REV-averaged velocity between the DNS and model results could potentially be further minimised by incorporating additional terms from the Darcy–Forchheimer equation (2.3).

3.3. Further application: prediction of overall time evolution

This subsection highlights the application of coupled macroscopic equations for predicting real-time changes in velocity and temperature. This functionality is critical for effective thermal management, as it necessitates continuous prediction to prevent excessive and uncontrolled temperature increases, thereby ensuring the safety and longevity of devices. The schematic geometry of porous media with two length scales in figure 4(a,b) is utilised. The geometry is considered a complete system, and no macroscopic turbulence emerges, thus the macroscopic equations can be simplified into ordinary differential equations by eliminating the volume $V$ , as

(3.1) \begin{equation} \partial _t{ {{{\bar {u}}}_{i}}}=\phi {{g}_{i}}-\phi {{\hat {R}}_{is}},\quad \partial _t{ \bar {T}}=\phi \dot {q}-\phi {{\hat {H}}_{s}}. \end{equation}

Here, $\bar {\varphi }$ represents the overall volume average. The effective overall drag ${\hat {R}}_{is}$ and interphase energy resistance ${\hat {H}}_{s}$ can be modelled using the same framework introduced in our previous works as shown in (2.3) and (2.10). The key distinction in this implementation is that we now have fully closed macroscopic equations, which permit direct computation of ${\hat {R}}_{is}$ and ${\hat {H}}_{s}$ , leading to a significant reduction in computational costs. This subsection primarily utilises the macroscopic equations (2.1) to track real-time temperature variations. The evolution of volume-averaged velocity and temperature is computed using the Dormand–Prince algorithm implemented as ode45 in MATLAB.

Figure 6.

The changes of coupled temperature and velocity versus time, for ${\textit{Re}_{s}}=793$ and $\textit{Pr}=1.0$ : (a) velocity, (b) temperature.

Figure 6 compares the DNS and model results, demonstrating that the coupling model effectively captures the growth of overall velocity and temperature with high consistency. Although a minor discrepancy exists, it may be attributed to the errors in ${\hat {R}}_{is}$ and ${\hat {H}}_{s}$ . Notably, the model coefficients are derived only from data obtained under steady or statistically steady flow states. The strong agreement between the model and DNS results suggests that the relaxation time scale for microscopic turbulence to reach an approximate steady state is shorter than the transient growth time scale of mean velocity and temperature. Additionally, it is important to note that temperature responds more slowly than velocity, as it requires time to adjust to variations in velocity. Overall, these coupled macroscopic equations provide a reliable framework for predicting real-time changes in both velocity and temperature.