2. Governing equations

We consider a domain composed of a porous region and a surrounding pure fluid one. The flow domain of interest is illustrated in figure 1. The working fluid is assumed to be an isotropic and Newtonian liquid. The solid matrix of the porous medium is assumed to be rigid with constant material properties. Also, it is internally heated, so it transfers heat to the interstitial fluid. Moreover, as mentioned above, thermal non-equilibrium is assumed between the solid matrix and the interstitial fluid, inside the porous region. In the following, quantities related to the solid matrix are denoted with the subscript ‘s’, whereas quantities related to the fluid phase are not denoted with a subscript.

Figure 1.

Illustration of the computational domain with the immersed porous region. The dimensions of the computational domain are $\hat {h}_{f} \times (\hat {h}_{p} + \hat {h}_{f}) \times \hat {h}_{f}$ and those of the porous region are $\hat {l}_{p} \times \hat {h}_{p} \times \hat {l}_{p}$.

Our numerical study is based on the single-domain approach and, more specifically, on the thermo-mechanical model of Papalexandris & Antoniadis (Reference Papalexandris and Antoniadis2015). This model relies on a mixture-theoretic formalism according to which both the solid matrix and the fluid phase are treated as coexisting but immiscible thermodynamic continua. Then, the porosity (volume fraction), $\phi$, is introduced as a concentration variable that provides, at each point of the entire domain, the density of volume occupied by the fluid. In mathematical terms, the axiomatic definition of the porosity as a density function follows from the Radon–Nikodym theorem (Goodman & Cowin Reference Goodman and Cowin1972; Varsakelis & Papalexandris Reference Varsakelis and Papalexandris2010). Inside the porous region, the porosity is less than unity, $\phi ({\boldsymbol {x}})< 1$; whereas in the pure fluid one, $\phi ({\boldsymbol {x}}) = 1$. Further, each phase is endowed with its own set of kinematic and thermodynamic variables and its own set of balance laws. Herein, the solid matrix is assumed to consist of non-reacting and rigid fibres at rest; therefore, its mass and momentum balance laws are trivially satisfied. Thus, the mathematical model consists of the mass, momentum and energy balance laws for the fluid phase and the energy balance law for the solid matrix.

Further, in the flow under study, the fluid density can vary but such variations are due to temperature gradients only and are not induced by mechanical means. Accordingly, the fluid velocity is very small with respect to the speed of sound and, therefore, compressibility effects are negligible. For this reason, there is no need to consider the full, compressible set of equations for mixed domains but its low-Mach-number approximation, which is also provided by Papalexandris & Antoniadis (Reference Papalexandris and Antoniadis2015). This approximation is derived by singularly perturbing the full compressible equations at the zero-Mach-number limit (Paolucci Reference Paolucci1982; Lessani & Papalexandris Reference Lessani and Papalexandris2006). According to it, the thermodynamic pressure is uniform while the momentum and energy equations of the fluid decouple. Thus the mass, momentum and energy balance laws for the fluid phase and the energy balance law for the solid matrix read, in dimensionless form,

(2.1)$$\begin{gather} \phi \frac{\partial \rho }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi \rho \boldsymbol{u}) = 0 , \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial \phi \rho \boldsymbol{u} }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\phi \rho \boldsymbol{uu}) + \phi \boldsymbol{\nabla} p = \frac{1}{Re} \boldsymbol{\nabla}\boldsymbol{\cdot}(\phi \mu \boldsymbol{V}) - \phi \rho Ri \tilde{\boldsymbol{y}} - \boldsymbol{\delta} \boldsymbol{u} , \end{gather}$$

(2.3)$$\begin{gather}\phi c_p \frac{\partial \rho T }{\partial t} + c_p \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\phi \rho T \boldsymbol{u}\right) = \frac{1}{Pr\,Re} \boldsymbol{\nabla}\boldsymbol{\cdot} (\phi k \boldsymbol{\nabla} T) + \tilde{h} (T_{s} - T) , \end{gather}$$

(2.4)$$\begin{gather}c_{p_{s}} (1-\phi) \rho_{s} \frac{\partial T_{s}}{\partial t} = \frac{1}{Pr\,Re} \boldsymbol{\nabla}\boldsymbol{\cdot} \left((1-\phi) \boldsymbol{k}_{s} \boldsymbol{\nabla} T_{s} \right) - \tilde{h} (T_{s} - T) + r . \end{gather}$$

In the above equations, $\rho$, $\boldsymbol {u} = (u, v, w)$, $p$ and $h$ denote respectively the density, velocity vector, dynamic pressure and specific enthalpy of the fluid, whereas $T_{s}$ denotes the temperature of the solid matrix. It is noted that $p$ represents the second-order term in the low-Mach-number expansion of the fluid pressure. The first-order term of this expansion, interpreted as the thermodynamic pressure of the fluid, is constant since the domain under study is open. Also, $\mu$, $k$ and $c_{p}$ are respectively the shear viscosity, thermal conductivity and specific heat capacity of the fluid. Further, $\tilde {{\boldsymbol {y}}}$ stands for the unit vector in the vertical direction (parallel to the gravity vector).

In (2.2), the deviatoric part of the viscous stress tensor $\boldsymbol {V}$ is given by

(2.5)\begin{equation} \boldsymbol{V} = (\boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^{\top}) - \tfrac{2}{3} \left( \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} \right) {\boldsymbol{I}}, \end{equation}

where ${\boldsymbol {I}}$ is the identity matrix. In the same equation, $\boldsymbol {\delta } \boldsymbol {u}$ represents the interphasial drag force. The drag parameter $\boldsymbol {\delta }$ is a second-order tensor since the solid matrix is, in general, an anisotropic medium. For the same reason, the heat conductivity of the solid matrix, $\boldsymbol {k}_{s}$, is a second-order tensor too. The term $\tilde {h} (T_{s} - T)$ in the energy equations (2.3) and (2.4) represents the heat transfer between the two phases, with $\tilde {h}$ being the heat transfer coefficient. Further, $r$ stands for the heat-source term of the solid matrix. It is noted that bulk viscous stresses have not been included in (2.2) since they are assumed to be negligible. It is also worth adding that for steady creeping flow inside a porous medium, the momentum equation (2.2) reduces to the Darcy–Forchheimer model wherein the permeability tensor is expressed in terms of $\phi$ and ${\boldsymbol {\delta }}$.

The dimensionless groups appearing in the above system are the Reynolds, $Re = {\hat {u}_{r} \hat {h}_{r}}/{\hat {\nu }_{r}}$, Richardson, $Ri = {\hat {g} \hat {h}_{r}}/{\hat {u}_{r}^2}$, and Prandtl, $Pr = {\hat {\nu }_{r}}/{\hat {\alpha }_{r}}$, numbers, where $\hat {u}_{r}$ is the reference velocity, $\hat {h}_{r}$ the reference length, $\hat {\nu }_{r}$ the reference kinematic viscosity, $\hat {\alpha }_{r}$ the reference thermal diffusivity and $\hat {g}$ the gravitational acceleration. Dimensional variables are denoted with the hat symbol $\, \hat {}\, $ throughout this work.

The governing system (2.1)–(2.4) is presented in dimensionless form for purposes of compactness and clarity. Also, the above dimensionless form of the governing system is precisely that which we have solved numerically in our simulations. However, it is well known that in natural convection, the relevant dimensionless group is the Rayleigh number and not the Reynolds or Richardson numbers. However, in the problem under study, the flow domain consists of two distinct parts: above and below the horizontal porous–pure fluid interface. Accordingly, there are two top–bottom temperature differences (one for each part of the domain) and two heights. The Rayleigh number, $Ra_{f}$, in the pure fluid region, i.e. $y > h_{p}$, is then given by

(2.6)\begin{equation} Ra_{f} = \frac{\hat{g} \hat{\beta}_{r} \Delta \hat{T}_{f} \hat{h}_{f}^3}{\hat{\alpha}_{r} \hat{\nu}_{r}} , \end{equation}

where $\hat {\beta }_{r}$ is the thermal expansion coefficient at the reference state and $\Delta \hat {T}_{f}$ is the mean temperature difference across the pure fluid region. However, $Ra_{f}$ is not really relevant to the lower domain, i.e. that encompassing the porous medium. For example, it would not be convenient to non-dimensionalise the energy equation of the solid matrix (2.4) with respect to $Ra_{f}$. However, the introduction of the Reynolds number facilitates the use of available correlations for the components of the drag parameter $\boldsymbol {\delta }$. Consequently, we resort to non-dimensionalising the governing system with respect to $Re$ and $Ri$. Nonetheless, the analysis of the results in the pure fluid domain is presented in terms of $Ra_{f}$.

The system (2.1)–(2.4) is closed with an equation of state for the fluid. In our study, the fluid is liquid water; we, therefore, employ an ‘isobaric equation of state’ for the density in the form of a $\rho -T$ relation. This is a fourth-order polynomial fit of the tabulated data (Wagner & Kretzschmar Reference Wagner and Kretzschmar2007) for water density at one atmosphere and over the temperature range of interest. The other parameters of water that appear in the governing system ($k$, $\mu$ and $c_p$) are also computed via a fourth-order polynomial fit of data taken from Wagner & Kretzschmar (Reference Wagner and Kretzschmar2007).

As mentioned in the introduction, for modelling purposes, we employ the one-domain approach. Accordingly, the governing equations of the fluid phase (2.1)–(2.3) are valid in both the porous and pure fluid regions, which eliminates the need for matching conditions at the interface between the porous and pure fluid regions. However, the energy balance law for the solid matrix (2.4) is valid only inside the porous region. An additional condition for $T_{s}$ is, therefore, necessary at the interface. Herein, we prescribe the following condition:

(2.7)\begin{equation} (1- \phi) ({\boldsymbol{k}}_{s} \boldsymbol{\nabla} T_{s}) \boldsymbol{\cdot} \tilde{\boldsymbol{n}} = 0, \quad \text{at the boundaries of the porous medium,} \end{equation}

where $\tilde {{\boldsymbol {n}}}$ is the unit vector normal to the interface. For a discussion on the appropriateness and applicability of this condition, the reader is referred to the papers by Papalexandris & Antoniadis (Reference Papalexandris and Antoniadis2015) and Antoniadis & Papalexandris (Reference Antoniadis and Papalexandris2016).

In our study, the porous medium is of uniform porosity $\phi$. With regards to the solid matrix, we assume that it consists of an ensemble of identical and uniformly distributed circular cylinders that are all aligned with the vertical $y$ direction. Accordingly, the solid matrix is an orthotropic medium. As such, $\boldsymbol {\delta }$ and $\boldsymbol {k}_{s}$, the interphasial drag and the solid thermal conductivity, are diagonal tensors and, moreover, their diagonal components in the $x$ and $z$ directions are equal, $\delta _{11} = \delta _{33}$ and $k_{s11} = k_{s33}$. Since the vertical cylinders composing the porous medium are separated from each other, we assume that there is no conductive heat transfer in the horizontal directions; therefore, $k_{s11} = k_{s33} = 0$. Then, for the remaining non-zero element of the conductivity matrix, we set $k_{s22} =k_{s}$, where $k_{s}$ is the thermal conductivity of the material from which the solid matrix is made. It is noted that these coefficients do not account for dispersive heat transport. In principle, incorporating a heat-dispersion tensor in the model would improve the accuracy of the modelling of heat conduction. However, available expressions of such a tensor are quite cumbersome (Whitaker Reference Whitaker1999; Lesseux & Valdés-Parada Reference Lesseux and Valdés-Parada2017) and, to our knowledge, there are no available numerical data that we could employ for the problem in hand. Herein, we consider a generic porous medium consisting of aligned cylinders and we thus expect dispersive effects to be limited. Accordingly, for simplicity purposes, we have opted not to include such effects in the parametrisation of the conductivity tensor.

For the parametrisation of $\delta _{11}$ and $\delta _{22}$, we proceed as follows. We consider the expression proposed by Koch & Ladd (Reference Koch and Ladd1997) for the average drag force $\hat {f}$ exerted on a single cylinder of diameter $\hat {d}_{c}$ at low Reynolds numbers:

(2.8)\begin{equation} \hat{f}=\hat{d}_{c} \hat{\mu}\hat{u} \left(\kappa_0+\kappa_2 Re_{d_{c}}^2\right), \end{equation}

where $Re_{d_{c}}$ corresponds to the particle Reynolds number based on the diameter of a single cylinder, $\hat {d}_{c}$, and the local horizontal velocity. Also, $\kappa _0$ and $\kappa _2$ are empirical coefficients, the values of which depend on the porosity $\phi$; their expressions are provided by Koch & Ladd (Reference Koch and Ladd1997) and are not reproduced herein for purposes of economy of space. The above expression for the horizontal force acting on one cylinder is then multiplied by the number density of cylinders in the porous medium, $\hat {N}_{c} = {4 (1 - \phi )}/{{\rm \pi} \hat {d}_{c}^2 \hat {h}_{p}}$. Upon non-dimensionalisation, the following relation for $\delta _{11}$ is derived:

(2.9)\begin{equation} \delta_{11} = \delta_{33} = \frac{1}{Re} \frac{4 \mu }{\rm \pi} \frac{(1-\phi)}{d_{c}^2} \left(\kappa_0 + \kappa_2 * Re_{d_{c}}^2 \right) . \end{equation}

Similarly, for the component $\delta _{22}$ of the drag coefficient tensor, we consider the wall stress in 2-D flow over a flat plate $\tau _{w}=\frac {1}{2} c_{f} \hat {\rho }\hat {v}^2$. Herein, $c_{f}$ is the average skin friction coefficient over the height of the cylinder $\hat {h}_{p}$, $c_{f}=1.328\sqrt {{\hat {\mu }}/{\hat {\rho }\hat {v}\hat {h}_{p}}}$ (Incropera et al. Reference Incropera, DeWitt, Bergman and Lavine2007). This expression is then multiplied by the surface area of a single cylinder and the number density of the cylinders $\hat {N}_{c}$ (Antoniadis & Papalexandris Reference Antoniadis and Papalexandris2015). Upon non-dimensionalisation, the following approximation for $\delta _{22}$ is derived:

(2.10)\begin{equation} \delta_{22} = \frac{2.656}{\sqrt{Re}} \frac{(1- \phi)}{d_{c}} \sqrt{\frac{\mu \rho |v|}{h_{p}}} , \end{equation}

where $Re$ is calculated with respect to the macroscopic reference length and velocity.

With respect to the above parametrisation for $\delta _{11}$ and $\delta _{22}$, it is noted that they are based on the premise that the cylinder Reynolds number is sufficiently small and that the flow is steady. In our simulations, the particle Reynolds number is indeed small (of the order of or less than unity) but the flow under study is unsteady. However, to our knowledge, reliable correlations for the drag coefficient over cylinder arrays for unsteady conditions are currently unavailable. For this reason, we resorted to employing (2.9) and (2.10).

With regards to the interphasial heat exchange coefficient $\tilde {h}$ that appears in the energy equations (2.3) and (2.4), correlations or experimental data are not currently available either. Herein, the heat transfer in such a configuration is assumed to be the sum of two contributions: one due to the cross-flow of the interstitial fluid and the other due to the flow parallel to the cylinders comprising the solid matrix. Accordingly, an approximation of $\tilde {h}$ is obtained by adding these two contributions. For the first contribution, $\tilde {h}_{\perp }$, we employ the heat transfer coefficient for flow around a cylinder and multiply it by the number density of cylinders $\hat {N}_{c}$. The final result reads

(2.11)\begin{equation} \tilde{h}_{{\perp}} = \frac{1}{Re \,Pr} \frac{4 (1 - \phi ) }{d_{c}^2} k \overline{Nu}_{d_{c}}, \end{equation}

where the Reynolds and Prandtl numbers are the macroscopic ones while $\overline {Nu}_{d_{c}}$ corresponds to the Nusselt number based on the cylinder diameter $\hat {d}_{c}$. The expression for the second contribution is obtained by multiplying the heat-transfer coefficient for flow parallel to a flat plate with the perimeter of a single cylinder and the number density of the cylinders. The final result reads

(2.12)\begin{equation} \tilde{h}_{{\parallel} } = \frac{1}{Re \,Pr} \frac{4 (1 - \phi ) }{d_{c} h_{p}} k \overline{Nu}_{h_{p}}. \end{equation}

For $\overline {Nu}_{d_{c}}$ and $\overline {Nu}_{h_{p}}$, we employ the following empirical correlations (Incropera et al. Reference Incropera, DeWitt, Bergman and Lavine2007):

(2.13)\begin{equation} \overline{Nu}_{d_{c}} = C Re_{d_{c}}^m Pr^{1/3} , \end{equation}

where $C$ and $m$ are functions of $Re_{d_{c}}$, and

(2.14)\begin{equation} \overline{Nu}_{h_{p}} = 0.453 Re_{h_{p}}^n Pr^{1/3} , \end{equation}

where $Re_{h_{p}}$ is the Reynolds number based on the height of the porous medium $\hat {h}_{p}$. The final expression for $\tilde {h}$ reads, in dimensionless form,

(2.15)\begin{equation} \tilde{h}=\tilde{h}_{{\perp}} + \tilde{h}_{{\parallel}} = \frac{1}{Re \,Pr} \frac{4 (1 - \phi ) }{d_{c}} k \left( \frac{\overline{Nu}_{d_{ c}}}{d_{c}} + \frac{\overline{Nu}_{h_{p}}}{h_{p}} \right) . \end{equation}

The algorithm employed for the numerical treatment of the governing system (2.1)–(2.4) is an extension of the single-phase low-Mach-number solver (Lessani & Papalexandris Reference Lessani and Papalexandris2006) to mixed domains. It is based on a predictor-corrector time-marching scheme and employs a projection method that yields a linear elliptic equation for the dynamic pressure. The discretisation of the elliptic equation gives a linear system that is solved via a preconditioned conjugate gradient iterative method, combined with an algebraic multigrid preconditioner implemented in the library PETSc (Balay et al. Reference Balay, Gropp, Curfman McInnes and Smith1997, Reference Balay2019). Spatial discretisation of the governing system is performed via a second-order finite-volume scheme on a collocated grid system. This is combined with a flux-interpolation technique to avoid the pressure odd–even decoupling that is typically encountered when using such grids (Lessani & Papalexandris Reference Lessani and Papalexandris2006, Reference Lessani and Papalexandris2008). For more details on the numerical algorithm, the reader is referred to Lessani & Papalexandris (Reference Lessani and Papalexandris2006) and Papalexandris & Antoniadis (Reference Papalexandris and Antoniadis2015). The algorithm has been parallelised using the message passage interface (MPI) protocol. Coordination of the communication between processors is performed via functions from the PETSc library.

3. Numerical set-up

The computational domain is a cuboidal cavity of equal length and width $\hat {h}_{f}$, and height $\hat {h}_{tot}$; see figure 1. In other words, its dimensions in the $x$, $y$ and $z$ directions are $\hat {h}_{f} \times \hat {h}_{tot} \times \hat {h}_{f}$. It encloses a porous region of equal length and width, $\hat {l}_{p}$, and height $\hat {h}_{p}$; i.e. the dimensions of the porous medium are $\hat {l}_{p} \times \hat {h}_{p} \times \hat {l}_{p}$. The porous medium is placed at the bottom of the domain, as illustrated in figure 1. The overlaying pure fluid domain is a cube of height $\hat {h}_{f}$. Accordingly, the total height of the domain, $\hat {h}_{tot}$, is given by $\hat {h}_{tot} = \hat {h}_{f} + \hat {h}_{p}$. In our study, the height of the porous region represents half of the height of the pure fluid domain above, $\hat {h}_{p} = 0.5 \hat {h}_{f}$, which brings the total height of the computational domain to $\hat {h}_{tot} = 1.5 \hat {h}_{f}$. In our study, $\hat {h}_{f}$ is set to 0.038 m and, therefore, the total height $\hat {h}_{tot}$ is 0.057 m.

For non-dimensionalisation purposes, $\hat {h}_{f}$ is used as the reference length. Using dimensionless coordinates, the lower boundary (bottom wall) is located at $y = 0$, the upper boundary (free surface) at $y=h_{tot} = 1.5$, and the horizontal interface between the porous and pure fluid regions at $y=h_{p} = 0.5$. The side walls of the computational domain are located at $x = 0$ and 1, and at $z = 0$ and 1. For simplicity purposes, the top boundary of the porous medium at $y=h_{p}$ will be referred to as the ‘horizontal interface’. The cubical subdomain above, $y>h_{p}$, will be referred to as the ‘pure fluid region’, while the domain below $y\leqslant h_{p}$ will be simply called the ‘lower part’ of the flow domain. As such, the lower part comprises the porous medium and the surrounding fluid between the porous medium and the side walls.

Concerning the physical parameters of the solid structure composing the porous medium, we assume that its porosity is equal to $\phi =0.6$ and that the (dimensionless) diameter of each cylinder of the solid matrix, $d_{c}$, is equal to 0.001. The material composing the cylindrical elements is assumed to be stainless steel; its dimensional thermophysical properties are listed in table 1. The variations of these properties over the temperature range considered herein are negligible; accordingly, we treat these properties as constant. These values are then non-dimensionalised by the corresponding reference variables of liquid water.

Table 1.

Physical properties of the solid matrix.

The ratio between the horizontal surface area of the porous region and the overall horizontal area of the domain is referred to herein as the ‘area coverage’ and is denoted by $a$; accordingly, $a=l_{p}^2 / h_{f}^2$. In our study, we consider three cases of different area coverage: Case 1 with $a=0.64$, Case 2 with $a=0.884$ and Case 3 with $a=1$. In these three cases, the heating of the solid matrix is uniform in all directions. In various applications of interest, such as SFPs, the details of the heat-load distribution in the solid matrix can play an important role in the emerging convective patterns. To study this role, we consider herein an additional case in which some parts of the solid matrix are uniformly heated while the remaining parts are not heated at all; this is Case 4 and the area coverage is the same as in Case 1, $a=0.64$. The area coverage of the porous medium and the heat-source distribution for all four cases are illustrated in figure 2.

Figure 2.

Porosity and heat-load distribution throughout a horizontal plane for four studied cases: (a) Case 1; (b) Case 2; (c) Case 3; and (d) Case 4. The inner rectangle represents the porous region and the line-filled areas indicate the parts of the porous medium that are heated. The heat-source term $r$ is constant along the vertical direction.

However, the total heat-supply rate, $\hat {Q}$, is kept the same in all cases. The relation between $\hat {Q}$ and the heat-source term $\hat {r}$ in (2.4) is $\hat {Q}=\int _{V_{p}} (1 - \phi ) \, \hat {r} \, \mathrm {d} V$, with $V_{p}$ being the volume of the porous domain. Accordingly, the value of $\hat {r}$ varies from one case to another. For uniform heating, Cases 1–3, $\hat {r}$ decreases as the area coverage $a$ increases. Similarly, in Case 4, the value of $\hat {r}$ in the heated parts is higher than in all other cases. Table 2 summarises the details of the heat-load distribution for all cases.

Table 2.

Area coverage and heat-supply details for the four cases considered in our study. In Case 4, the heat-source term is positive in the heated parts of the solid matrix and zero in the non-heated ones.

With regards to the boundary conditions, free slip is prescribed at the upper boundary, $y = h_{tot}$, so as to approximate the free surface of a liquid. The no-slip kinematic condition is prescribed at the bottom and side walls. Concerning the thermal boundary conditions, the temperature of the top wall, $T_{top}$, remains fixed and uniform, while the bottom and side walls are assumed to be adiabatically isolated. As for initial conditions, we assume that the fluid is initially at rest and of uniform temperature $T= T_{top}$ in the entire computational domain. The initial temperature of the solid matrix is also uniform and equal to $T_{top}$.

In our simulations, we start taking statistics once the flow becomes fully developed and statistically steady, based on the area-averaged temperature at the plane of the horizontal interface and the bottom wall. In terms of notation, the time-average of a generic quantity $\varphi$ is denoted by $\langle \varphi \rangle _{t}$ and the area-average by $\langle \varphi \rangle _{xz}$. The time and area average is denoted by $\langle \varphi \rangle$. In particular, $\langle \hat {T}_{int} \rangle$, denotes the time-and-area-averaged temperature over the entire plane of the horizontal interface, $y = h_{p}$.

The reference temperature, $\hat {T}_{r}$, is chosen to be the arithmetic mean between $\hat {T}_{top}$ and $\langle \hat {T}_{int} \rangle$, i.e. $\hat {T}_{r} = ( \hat {T}_{top} + \langle \hat {T}_{int} \rangle ) / 2$, so as to facilitate comparisons with classical RBC. It is noted that the reference values of the material properties are taken at $\hat {T}_{r}$ and at an ambient pressure of 1 bar. The values of the reference, top and bottom temperatures are listed in table 3, together with the values of $\langle \hat {T}_{int} \rangle$. In our analysis, we will also make use of the temperature difference across the pure fluid region, $\Delta T_{f} = \langle \hat {T}_{int} \rangle - \hat {T}_{top}$, the values of which are also provided in this table.

Table 3.

Time-and-area-averaged temperatures at three planes, namely the top ($y = h_{tot}$), interface ($y = h_{p}$) and bottom planes ($y = 0$). Here, $\hat {T}_{top}$ is fixed for all cases and $\langle \hat {T}_{int} \rangle$ and $\langle \hat {T}_{bot} \rangle$ are the time-and-area-averaged temperatures at the plane of the horizontal interface and at the bottom wall. Additionally, $\Delta T_{f}$ stands for the temperature difference across the pure fluid region, $\Delta T_{f} = \langle \hat {T}_{int} \rangle - \hat {T}_{top}$.

Also, the free-fall velocity, $\hat {u}_{r} = \sqrt {\hat {\beta }_{r}\, \hat {g}\, \hat {h}_{f}\, \Delta T_{f}}$, and free-fall time, $\hat {t}_{r} = \hat {h}_{f} / \hat {u}_{r}$, are chosen as the reference velocity and time unit respectively, as is typical in studies of natural convection. The reference values of the various thermophysical properties (density $\hat {\rho }_{r}$, thermal expansion coefficient $\hat {\beta }_{r}$, kinematic viscosity $\hat {\nu }_{r}$, thermal conductivity $\hat {k}_{r}$ and thermal diffusivity $\hat {k}_{r}$) are those of the fluid (liquid water) at the reference temperature $\hat {T}_{r}$. These values, along with the free-fall velocity and the Prandtl number, are listed in table 4 for all cases of our study. For completion purposes, it is worth adding that for Cases 1–4, $Re=1333$, 1435, 1762 and 1224, respectively; also $Ri=388$, 342, 245 and 443, respectively.

Table 4.

Reference values of the various thermophysical parameters and Prandtl number for the cases considered herein.

The values of $Ra_{f}$ for all cases of our study, as determined by the numerical simulations, are provided in table 5. In the same table, we also provide the values of the Nusselt number in the pure fluid region $Nu_{f}$.

Table 5.

Time-averaged results. Here, $Ra_{f}$ and $Nu_{f}$ denote respectively the Rayleigh and Nusselt numbers in the pure fluid region.

3.1. Definition of the Nusselt number

In the problem at hand, convection is driven by the internal heating of the solid matrix which is modelled as a volumetric effect. To provide a definition of the Nusselt number across the entire domain that can also be used for comparisons with RBC in a meaningful manner, we introduce the area-averaged heat supply $J_{y}$:

(3.1)\begin{equation} J_{y} = \int_0^y \langle r \rangle_{xz} \, {{\rm d} y}^\prime, \end{equation}

which is known beforehand.

In our case, since the heat source $\langle r \rangle _{xz}$ is a Heaviside step function of the $y$ coordinate, $J_{y}$ is simply a ramp function; i.e. it increases linearly in the lower part of the flow domain and remains constant in the pure fluid region. Crucially, in fully developed flow, $J_{y}$ is equal to the total heat flux across a horizontal plane. As such, it can be obtained by time and area averaging the sum of the energy equations (2.3) and (2.4) and integrating the result between the bottom wall and a horizontal plane:

(3.2)\begin{equation} J_{y} = \int_0^y \langle r \rangle_{xz} \, {{\rm d} y}' = \underbrace{\langle c_{p} \phi \rho v T \rangle}_{J_{conv}} \ \underbrace{- \frac{1}{Re \,Pr} \langle \phi k \boldsymbol{\nabla} T \rangle}_{J_{ diff}} \ \underbrace{- \frac{1}{Re\, Pr} \langle (1 - \phi) \boldsymbol{k}_{s} \boldsymbol{\nabla} T_{s} \rangle}_{J_{diff,s}}. \end{equation}

Accordingly, the Nusselt number can be defined on the basis of $J_{y}$ and the diffusive contribution, $\overline {J_{diff}}$, averaged over the pure fluid region:

(3.3)\begin{equation} \overline{J_{diff}} = \frac{1}{h_{f}} \, \int_{h_{p}}^{h_{tot}} J_{diff} \, {{\rm d} y} =- \frac{1}{Re \,Pr} \frac{k \Delta T_{f}}{h_{f}}. \end{equation}

More specifically, the local Nusselt number $Nu_{y}$ is given as the ratio between $J_{y}$ and $\overline {J_{diff}}$:

(3.4)\begin{equation} Nu_{y} = \frac{J_{y}}{\overline{J_{diff}}} = \underbrace{\frac{J_{ conv}}{\overline{J_{diff}}}}_{Nu_{conv}} + \underbrace{\frac{J_{ diff}}{\overline{J_{diff}}}}_{Nu_{diff}} + \underbrace{\frac{J_{ diff,s}}{\overline{J_{diff}}}}_{Nu_{diff,s}}. \end{equation}

According to the above relation, $Nu_{y}$ is the sum of the convective and diffusive components related to the fluid and the diffusive component related to the solid. It is important to note that, in the above equation, the numerator of $Nu_{y}$ is known a priori, but the denominator is not because it involves the temperature difference $\langle T_{int} \rangle - T_{top}$. This is opposite to RBC wherein the temperature difference across the domain is known beforehand but not the heat flux $J_{y}$.

Finally, the global Nusselt number in the pure fluid region $Nu_{f}$ is the depth average of $Nu_{y}$ for $y>h_{p}$. In this region, $Nu_{y}$ is constant along $y$ and equal to $Nu_{f}$. The values of $Nu_{f}$ for all cases considered herein are provided in table 5.

3.2. Grid-resolution criteria

The present study consists of DNS and, as such, the computational grid has to be sufficiently refined. Typically, in classical RBC, two resolution criteria are employed to ensure the adequacy of the grid resolution (Hay & Papalexandris Reference Hay and Papalexandris2020). The first concerns the resolution of the hydrodynamic and thermal boundary layers (Shishkina et al. Reference Shishkina, Stevens, Grossmann and Lohse2010), and the second one the cell size in the bulk (Grötzbach Reference Grötzbach1983).

With regards to the near-boundary regions, we employ the criterion of Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010) for the minimum number of cells in the hydrodynamic and thermal boundary layers. These criteria require the knowledge of $Nu_{f}$. In our study, $Nu_{f}$ is estimated via a preliminary simulation for each case. These preliminary simulations are run until the flow becomes fully developed and statistically steady, based on the area-averaged temperature at the plane of the horizontal interface $y = h_{p}$. Subsequently, $Nu_{f}$ is computed by taking statistics over 200 free-fall times. Instantaneous distributions of the flow variables from the preliminary simulations are then used for the initialisation of the fully resolved DNS via grid interpolation.

The criterion of Grötzbach (Reference Grötzbach1983) for classical RBC is based on the principle that the cell size in the bulk must be at most equal to the smallest of the Kolmogorov and Corrsin scales of the flow turbulence, $\eta = (\nu ^3 / \epsilon )^{{1}/{4}}$ and $\eta _{c} = \eta /Pr^{{3}/{4}}$, respectively, with $\epsilon$ being the kinetic-energy dissipation. Then the computational mesh is stretched according to a hyperbolic-tangent function from the boundaries to the bulk of the pure fluid region so as to satisfy this criterion. Concerning the lower part of the flow domain, $y \leqslant h_{p}$, the resolution at the bottom wall also meets the requirements of Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010). Also, for purposes of numerical stability, the porosity jump across the porous–pure fluid interfaces has been smeared across two computational cells. Also, as the macroscopic interface corresponds to a discontinuity of the porosity, the number of cells in that area is increased consequently. A schematic representation of the computational mesh is shown in figure 3, showing the regions of refinement and the relative size of the cells throughout the computational domain.

Figure 3.

Schematic representation of the computational mesh, showing the regions of refinement and the relative size of the cells throughout the computational domain. For visualisation purposes, only a fraction of the computational cells is shown.

The number of cells in the three spatial directions for the fully resolved DNS, denoted $N_{x}$, $N_{y}$ and $N_{z}$, are provided in table 6. In the same table, we also provide the number of cells inside the hydrodynamic and thermal boundary layers, $N_{u}$ and $N_{th}$, respectively. In those columns, the numbers in parentheses are the minimum requirements according to the criterion of Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010). According to our numerical simulations, in Case 1, the area-averaged Kolmogorov scale in the pure fluid domain is estimated at $\eta =1.1 \times 10^{-2}$ near the porous–pure fluid interface and $1.8\times 10^{-2}$ in the middle plane ($y=1$). While the Corrsin scale $\eta _{c}$ is equal to $3\times 10^{-3}$ near the porous–pure fluid interface and $4.7\times 10^{-3}$ in the middle plane. Based on the values of the cell sizes listed in table 6, we deduce that our simulations sufficiently resolve all relevant hydrodynamic and thermal scales. The smallest thermal fluctuations in the bulk are not resolved; however, the energy carried by these is exceedingly small.

Table 6.

Grid-resolution criteria. Here, $N_{x}$, $N_{y}$ and $N_{z}$ respectively are the number of cells in the $x$, $y$ and $z$ directions. Additionally, $N_{u}$ and $N_{th}$ are the number of cells in the hydrodynamic and thermal boundary layers, with numbers in parentheses providing the minimum requirement according to Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010). The terms $\Delta y_{min}$ and $\Delta y_{max}$ respectively correspond to the smallest and largest cell size, with the numbers in parentheses providing the maximum cell size according to Grötzbach (Reference Grötzbach1983).

Finally, it is noted that in Cases 1, 2 and 4 (which, as it turns out, are the cases where the horizontal flow inside the porous medium is not negligible), the minimum grid resolution inside the porous domain is $\Delta x = \Delta z = 3.11\times 10^{-5}$, $7.57\times 10^{-6}$ and $2.71\times 10^{-5}$, respectively. On the basis of these values and those for the porosity $\phi =0.6$ and cylinder diameter $d_{c}=0.001$, it is readily found that for Cases 1, 2 and 4, the minimum coverage is 15.8, 3.8 and 14.8 cylindrical elements per cell, respectively. In turn, this implies that the flow scales resolved in our simulations are considerably larger than the characteristic scale of the microstructure of the porous medium.

4. Numerical results and discussion

The results section is divided into two parts. In the first part, we present the results for Case 1 and analyse the flow properties and convective structures in both the pure fluid region and the lower part of the domain. In the second part, we present the results for the other cases and make comparisons between them. All statistical results presented below are time-averaged over a period of 300 free-fall times after the flow has become statistically steady. The simulation time for the averaging period of 300 $t_{r}$ is approximately 25 000 h per core in a cluster of Intel-Xeon-Gold-6252${\circledR}$ processors.

4.1. Results of Case 1

4.1.1. Global flow properties

We first examine the distribution of the normalised temperature $\theta$, defined as

(4.1)\begin{equation} \theta = \frac{T - T_{top}}{\langle T_{bot} \rangle - T_{top}}. \end{equation}

According to this definition, $\langle \theta \rangle$ at the bottom plane is equal to unity. Plots of the time-averaged normalised temperature $\langle \theta \rangle _{t}$ at the top, horizontal-interface and bottom planes are shown in figure 4(a). At the bottom plane, the fluid inside the porous medium is heated substantially, as expected, while the temperature of the fluid outside remains quite low, thereby leading to large thermal gradients across this plane. In particular, $\langle \theta \rangle _{t}=2.3$ at the centre of the plane, while $\langle \theta \rangle _{t}=0.22$ near the side walls. Further, from the plots at the bottom plane and horizontal interface, we readily infer that the fluid temperature inside the porous medium is not uniform either, even though the solid matrix is uniformly heated. Evidently, this is due to fluid motion inside the porous medium.

Figure 4.

Case 1: (a) time-averaged normalised temperature, $\langle \theta \rangle _{t}$, at the top ($y = h_{tot}$), horizontal-interface ($y = h_{p}$) and bottom planes ($y = 0$); (b) superimposed vector plots of the time-averaged velocity at the top plane (free surface), showing the impingement point and the horizontal motion of the fluid. The red dashed line denotes the vertical plane parallel to the side walls that contains the impingement point; this is the point from which the velocity arrows start to diverge.

In figure 4(b), we present vector plots of the velocity at the free surface. We observe that the impingement point of the upward fluid motion is not aligned with a diagonal plane, as is often the case in turbulent RBC at this range of Rayleigh numbers. Instead, it is located closer to the centre of the free surface. Then, the fluid reaches the free surface at the impingement point and moves horizontally (‘spreads’) in all directions. Subsequently, it starts descending once it approaches the side walls.

In figure 5, we present different views of time-averaged streamlines coloured by the vertical velocity component. Therein, we can discern the impingement point and the downward motion of the fluid close to the side walls once it leaves the free surface. We also note that the fluid originating from the corners of the free surface continues to descend below the plane of the horizontal interface and reaches the bottom wall. In the lower part of the domain, $y< h_{p}$, the descending fluid starts to circumvent the porous medium and, therefore, its vertical velocity decreases. Also, from figure 5(c), we infer that the circumventing currents eventually reach the side boundaries of the porous medium. Subsequently, such a current splits into two parts. One part actually enters the porous medium, while the other moves upwards along its side faces. In the same figure, we can readily discern the perimeter of the porous medium because of the change in the direction of the streamlines: inside the porous medium, they point towards the centre. In fact, near the bottom and in the periphery of the porous medium, the fluid moves mostly horizontally and towards the centre. However, in the central region, the streamlines look intertwined, which is indicative of vertical motion (plumes) due to thermal convection.

Figure 5.

Case 1: time-averaged streamlines coloured by the time-averaged vertical velocity component, $\langle v \rangle _{t}$. In these plots, $\langle v \rangle _{t}$ has been non-dimensionalised by $\bar {\alpha }/\hat {h}_{tot}$, where $\bar {\alpha }$ is the thermal diffusivity of water at $(\hat {T}_{top} + \langle \hat {T}_{bot} \rangle )/ 2$. (a,b) Side views. (c) Bottom view.

Further, the fluid moving upwards along the side boundaries of the porous medium eventually reaches the plane of the horizontal interface, $y = h_{p}$, and gets entrained in the convective structures of the pure fluid region $y>h_{p}$. For example, the vortex shown in figure 5(b) and its counterpart in figure 5(b) appear to have been formed by the interaction of this ascending fluid and the downward stream close to the side walls.

In figure 6, we show colour plots of the time-averaged vertical velocity component $\langle v \rangle _{t}$ on two vertical planes, superimposed with the time-averaged in-plane velocity vectors. The first plane, $p_1$, is parallel to the side walls and includes the impingement point at the free surface, while the second one is the diagonal plane $d_2$, as depicted in figure 4(b). According to these figures, the flow in the pure fluid domain organises itself in a large-scale circulation (LSC), which is one of the principal coherent structures in turbulent convection. In this case, the impingement point is not located at a diagonal plane and therefore, the resulting LSC does not have the single-roll structure that is typically encountered in classical RBC. However, the impingement point is not located at the centre of the free surface either and, therefore, a dual-roll LSC such as that reported by Hay et al. (Reference Hay, Martin, Migot and Papalexandris2021) is not formed either. Instead, the emerging LSC has an intermediate structure between those of the one roll and dual roll; at $p_1$, it comprises a single roll whereas at $d_2$, it comprises two rolls. It is noted that reorientation events of the LSC occur over time. During the simulation, we monitored the evolution of the LSC by placing probes as was done by Foroozani et al. (Reference Foroozani, Niemala, Armenio and Sreenivasan2017) and Marichal & Papalexandris (Reference Marichal and Papalexandris2022), and did not record any reorientation events. Nonetheless, we expect that such events will occur at later times, in which case, the location of the impingement point will change.

Figure 6.

Case 1: contour of time-averaged vertical velocity $\langle v \rangle _{t}$ superimposed with time-averaged in-plane velocity vectors: (a) plane $p_{1}$; (b) diagonal plane $d_{2}$. Here, $\langle v \rangle _{t}$ is non-dimensionalised by $\bar {\alpha }/\hat {h}_{tot}$. The black line represents the boundary of the porous medium.

It is interesting to observe that, based on the velocity amplitudes, the upward motion at $p_1$ in the pure fluid region is almost twice as strong as at $d_2$. This is explained by the fact that the upward-moving fluid at $p_1$ is right below the impingement point, i.e. at the core of the ascending plume. However, the motion of the descending fluid close to the side walls is stronger in $d_2$ than in $p_1$. In particular, at $d_2$, the downward plume close to the walls splits at the plane of the horizontal interface; one part continues to descend while the other part turns and follows the LSC as it moves along the horizontal interface. Whereas at $p_1$, most of the descending fluid joins the LSC so that the downward motion in the lower part, $y< h_{p}$, is very weak. Part of this is due to the fact that at the corners, the distance between the side walls and the vertical boundaries of the porous medium is larger so that boundary-layer effects are comparatively less pronounced. Further, from the velocity amplitude, we observe that, in this case, the ascending motion in the pure fluid region is stronger than the descending one.

With regards to the bottom part, in figure 6, we discern the upward motion of cold fluid along the vertical (side) boundaries of the porous medium that was mentioned above. Accordingly, this motion is mostly due to inertial effects rather than buoyancy. We can also identify some weak convective structures inside the porous medium. These are alternating ascending and descending plumes. Finally, from the in-plane velocity plots, we infer that the horizontal velocity inside the porous medium is very low, which in turn implies strong shear at the horizontal interface.

4.1.2. Analysis of the flow properties in the pure fluid region ($y > h_{p}$)

We proceed to describe in more detail the flow properties in the pure fluid region $y> h_{p}$, i.e. in the subdomain above the plane of the horizontal interface. To this end, we first introduce the normalised temperature $\theta _{f}$, defined by

(4.2)\begin{equation} \theta_{f}= \frac{T - T_{top}}{\langle T_{int} \rangle - T_{top}}. \end{equation}

It is noted that $\theta _{f}$ is defined with respect to the mean temperature at the horizontal interface, whereas $\theta$ is defined with respect to the temperature at the bottom wall, so as to allow for comparisons with the corresponding RBC case.

In figure 7, we present plots of the time-averaged normalised temperature $\langle \theta _{f} \rangle _{t}$ at the two diagonal planes $d_1$ and $d_2$, superimposed with time-averaged in-plane velocity vectors. At these two planes, the LSC consists of two large rolls. As mentioned above, the descending plumes are located close to the side walls, whereas the two rolls merge into a strong ascending plume in the region underneath the impingement point. As expected, the temperature peaks occur right above the horizontal interface. We also observe that, in both planes, there are patches of relatively cold fluid inside the rolls of the LSC, located below the mid-plane of the pure fluid region. They consist of fluid originating from the descending plumes and of fluid that has moved upwards along the vertical boundaries of the porous medium. The fluid above these patches is at a higher temperature. By contrast, such a pattern is not observed when the porous medium extends throughout the horizontal plane (so that its area coverage is unity).

Figure 7.

Case 1: colour plots of the time-averaged normalised temperature $\langle \theta _{f} \rangle _{t}$ in the pure fluid region, superimposed with time-averaged in-plane velocity vectors: (a) diagonal plane $d_1$; (b) diagonal plane $d_2$. The colour scale is the same for both panels.

Next, we present flow statistics in the pure fluid region. Also, we make comparisons with the equivalent RBC case, i.e. convection between a uniformly heated bottom wall and a cooler free surface (with fixed temperature) at the same Rayleigh and Prandtl numbers. Herein, this equivalent case is denoted by RBC-FS.

In figure 8(a), we have plotted the profiles of the mean (area and time-averaged) normalised temperature $\langle \theta _{f} \rangle$. The profile for the problem in hand shows large thermal gradients right above the horizontal interface, which implies a significant diffusive heat flux from the porous medium to the pure fluid region above. We further observe that in the RBC-FS, the profile is ‘shifted’ to the left (Hay & Papalexandris Reference Hay and Papalexandris2019), i.e. towards the temperature of the free surface. Whereas in classical RBC (between two horizontal walls), the mean temperature profile is centred. This is due to the absence of a hydrodynamic boundary layer at the free surface, which promotes convective motions of the fluid therein. However, in our case, the mean temperature profile in the pure fluid region is slightly ‘shifted’ to the right, i.e. towards the mean temperature at the plane of the horizontal interface, despite the presence of the free surface. As it will be shown below, this shift of $\langle \theta _{f} \rangle$ to the right is not observed in Case 3 in which the area coverage of the porous medium is unity. Therefore, it cannot be attributed solely to the presence of the porous medium and, instead, is explained as follows. The porous medium inhibits significantly the horizontal motion of fluid due to drag. Therefore, at the horizontal interface, at $y=h_{p}$, there is a shear layer whose effect on the convective motions is analogous to that of boundary layers at a rigid wall. However, when the area coverage of the porous medium decreases, the shear effects at the horizontal interface decrease too, thereby enhancing convective motions of the nearby fluid. Additionally, convective motions are enhanced by the flow in the side gaps of the lower part of the domain. In turn, these effects result in a shift of the temperature profile in the pure fluid domain towards the temperature at the plane of the interface $\langle T_{int} \rangle$. From the profile in figure 8(a), we observe that the effect of the horizontal interface has counterbalanced the effect of the free surface.

Figure 8.

Profiles of the normalised temperature in the pure fluid region, $y \geqslant h_{p}$: (a) mean profile, $\langle \theta _{f} \rangle$, where the thick vertical line represents the arithmetic mean, $\langle \theta _{f} \rangle =0.5$, between the values at the top and at $y=h_{p}$; (b) r.m.s. profile, $\langle \theta _{rms,f} \rangle$, where the legend is Case 1 (——) and RBC-FS ($\cdots \cdots$).

Further, in the same figure, we observe that $\langle \theta _{f} \rangle$ is not monotonous with $y$ but, instead, has a lower local minimum and an upper local peak. Such overshoots have been observed previously in the soft-turbulence regime; see, for example, Horn, Shishkina & Wagner (Reference Horn, Shishkina and Wagner2013). In our case, however, the value of $Ra_{f}$ indicates that we are already in the hard-turbulence regime. Therefore, we attribute these overshoots to the presence of patches of colder fluid as mentioned above.

In figure 8(b), we present profiles of the root mean square (r.m.s.) of temperature fluctuations. In both Case 1 and RBC-FS, the profile has a local peak below the free surface due to the thermal boundary layer. Subsequently, it drops, as thermal mixing becomes more efficient. Then, in the bulk of the domain, the gradients of the temperature r.m.s. are much smaller, as is typically the case in thermal convection. However, in Case 1, close to the plane of the horizontal interface, the fluctuations increase substantially and attain a maximum at $y=h_{p}$. However, in the corresponding RBC-FS case, the fluctuations peak at the edge of the bottom thermal boundary layer and vanish at the wall. The high peak in Case 1 is attributed primarily to the limited area coverage of the porous medium which gives rise to the development of the secondary flow between the solid walls and the porous medium, thereby resulting in large horizontal temperature gradients, as depicted in figure 4(a).

Profiles of the r.m.s. fluctuations of the vertical and horizontal velocity components are shown in figures 9(a) and 9(b), respectively. In the bulk, the $v_{rms}$ profiles have the same shape for both Case 1 and RBC-FS. In particular, $v_{rms}$ increases with $y$ and peaks between the mid-plane and the top free surface. Subsequently, it decreases until it vanishes at the top. The fact that the peaks are located above the mid-plane is due to the difference in the strengths of the ascending and descending motions. This has been documented before (Zikanov et al. Reference Zikanov, Slinn and Dhanak2002; Hay & Papalexandris Reference Hay and Papalexandris2019) and is attributed to the presence of the top free surface. Further, we observe that while in RBC-FS, $v_{rms}$ vanishes at the bottom wall, in Case 1, its value at the plane of the horizontal interface is substantial. As a result, the values of $v_{rms}$ are everywhere higher in Case 1 than in RBC-FS. This, in turn, implies stronger plumes (stronger vertical motion) in Case 1. Once again, this is attributed to the limited area coverage of the porous medium that enhances convection.

Figure 9.

Profiles of the r.m.s. of velocity fluctuations: (a) $v_{rms}$ and (b) $u_{rms}$. The velocities have been non-dimensionalised by $\hat {\alpha }_{r} / \hat {h}_{f}$. The legend is Case 1 (——) and RBC-FS ($\cdots \cdots$).

Plots of the r.m.s. of the horizontal velocity, $u_{rms} = \sqrt {\langle u^2\rangle + \langle w^2\rangle - \langle u\rangle ^2 - \langle w\rangle ^2}$ are provided in figure 9(b). As expected, similarly to the r.m.s. of the vertical velocity, in Case 1 $u_{rms}$ does not vanish at the plane of the horizontal interface, whereas it does vanish at the bottom wall in RBC-FS. Interestingly, in Case 1, it increases rapidly and peaks at a short distance above the horizontal interface. This is due to the shear layer that is developed across this interface. Further, in both Case 1 and RBC-FS, $u_{rms}$ is approximately constant in the bulk of the domain but becomes larger close to the free surface. Actually, in both cases, $u_{rms}$ attains its maximum value right at the top. This peak is a standard feature of turbulent convection with a free surface on top (Zikanov et al. Reference Zikanov, Slinn and Dhanak2002; Hay et al. Reference Hay, Martin, Migot and Papalexandris2021; Marichal & Papalexandris Reference Marichal and Papalexandris2022) and is due to the strong vertical motion of fluid beneath the free surface which, in the absence of a hydrodynamic boundary layer, leads to strong surface currents. In other words, the absence of a hydrodynamic boundary layer greatly reduces viscous dissipation in the vicinity of the free surface. The fluid coming from below has been accelerated by the buoyancy force and, once it reaches the free surface, it turns while maintaining most of its momentum and ‘spreads’ to all horizontal directions; this, in turn, leads to the aforementioned surface currents.

Lastly, in figures 10(a) and 10(b), we respectively plot the convective and diffusive contributions of the local Nusselt number. As is typically the case, $Nu_{conv}$ is constant in the bulk and decreases near the boundaries. However, in Case 1, the lower limit of the pure fluid region, $y=h_{p}$, is not a boundary but simply encompasses the horizontal interface. Accordingly, at that plane, the vertical velocity is not zero and the temperature is not uniform. This implies that $Nu_{conv}$ does not vanish therein. As a result, $Nu_{conv}$ is higher in Case 1 than in RBC-FS throughout the pure fluid domain. Nonetheless, in both cases, $Nu_{conv}$ vanishes at the free surface. As a result, and in contrast to RBC-FS, the profile of $Nu_{conv}$ in Case 1 is not symmetric with respect to the mid-plane.

Figure 10.

Profiles of the time-and-area-averaged components of the local Nusselt number $Nu_{y}$ along the $y$ direction. (a) Convective component, $Nu_{conv}$. (b) Diffusive component, $Nu_{diff}$. The legend is Case 1 (——) and RBC-FS ($\cdots \cdots$).

With regards to the diffusive component $Nu_{diff}$ shown in figure 10(b), it follows the opposite trend, i.e. it decreases in the bulk and increases near the boundaries, as expected in turbulent convection. It is worth noting that $Nu_{diff}$ is slightly negative in the bulk due to the overshoots of the mean temperature profile. In this manner, the local Nusselt number, $Nu_y = Nu_{conv}+Nu_{diff}$, remains constant along $y$ and equal to the global Nusselt number $Nu_{f}$. Accordingly, not only $Nu_{conv}$ but also $Nu_y$ is higher in Case 1 than in RBC-FS. Also, since similarly to $Nu_{conv}$, the profile of $Nu_{diff}$ is symmetric only in RBC-FS but not in Case 1.

4.1.3. Analysis of the flow properties in the lower part ($y \leqslant h_{p}$)

Herein, we examine in more detail the flow properties in the lower part of the domain, which includes the porous medium. To this end, we first define the normalised temperatures in this part for the two phases via

(4.3a,b)\begin{equation} \theta_{p} = \frac{T - \langle T_{int} \rangle}{\langle T_{bot} \rangle - \langle T_{int} \rangle}\quad \textrm{and} \quad \theta_{s,p} = \frac{T_{s} - \langle T_{int} \rangle}{\langle T_{bot} \rangle - \langle T_{int} \rangle}. \end{equation}

In figure 11, we present colour plots of the time-averaged normalised temperatures, $\langle \theta _{p} \rangle _{t}$ and $\langle \theta _{s,p} \rangle _{t}$, respectively at the vertical plane $p_1$ and the diagonal plane $d_2$. The most noticeable feature is that, inside the porous medium, there are strong thermal gradients in the horizontal direction even though the solid matrix is uniformly heated. As mentioned above, this is due to cold fluid entering from the vertical faces of the porous medium but also due to convective plumes in its interior. Convective motion inside the porous medium is much weaker than in the pure fluid region above because of interphasial drag. Also, the convective rolls are long and thin (i.e. finger-like), which implies that the fluid motion therein is almost vertical, i.e. with a very small horizontal velocity. This is explained by the fact that with this particular microstructure of the solid matrix (vertically aligned cylinders), the interphasial drag is much higher in the horizontal direction than in the vertical one, which favours vertical motion. However, since fluid enters the porous medium sideways, this means that mixing between the fluid in the central region of the porous medium (where the plumes are located) and the fluid outside the central region is reduced. Further, from figure 11, we infer that the plumes reach up to the horizontal interface and interact with the shear layer formed right above it.

Figure 11.

Case 1: time-averaged normalised temperatures of the fluid $\langle \theta _{p} \rangle _{t}$ (top) and the solid matrix $\langle \theta _{s,p} \rangle _{t}$ (bottom), superimposed with time-averaged in-plane velocity vectors. (a) Vertical plane $p_1$. (b) Diagonal plane $d_2$. The arrows mark the fluid velocity, but are enhanced in the bottom panels for better visibility.

Another interesting aspect is that the temperatures of the two phases are practically equal, which implies efficient heat transfer from the solid matrix to the interstitial fluid. For example, as cold fluid enters the porous medium from the sides, it absorbs heat from the solid matrix. Consequently, $\theta _{s}$ is much lower near the sides of the porous medium than in its interior. It is noted, however, that these observations concern a fully developed, time-averaged statistically steady flow. In particular, in the early stages of the evolution of the flow, the difference between the temperatures of the two phases is significant.

The velocity vector plots in figure 11 reaffirm the aforementioned observation that cold fluid descends mostly from the corners of the domain and then circumvents the porous medium. Accordingly, a part of it enters the porous medium while the rest starts ascending along its vertical sides. Further, the fluid entering the porous medium changes direction rapidly and moves upwards even though its temperature is still relatively low. In other words, the upward motion of the fluid near the side boundaries of the porous medium is mostly driven by inertial effects rather than by buoyancy. It is noted that the velocity vectors shown in figure 11 are used mainly for visualisation purposes of the flow direction. More specifically, the size of the velocity vectors in the bottom subfigures has been made larger than in the top ones for better visibility.

In figure 12, we have plotted the profile of the mean (time and area-averaged) normalised temperature of the fluid, $\langle \theta _{p} \rangle$. We can observe that $\langle \theta _{p} \rangle$ remains almost constant up to $y\approx 0.13$, which roughly corresponds to 1/4 of the height of the porous medium. The absence of a lower thermal boundary layer is due to the adiabatically isolated bottom wall. Then, as $y$ increases, $\langle \theta _{p} \rangle$ drops. Also, the rate of decrease of $\langle \theta _{p} \rangle$ gets higher with $y$. Interestingly, similar temperature profiles have been reported by Thirlby (Reference Thirlby1970) who studied convection in an internally heated layer with an isolated bottom wall.

Figure 12.

Plots of the normalised temperature in the lower part of the domain: time- and area-averaged temperature of Case 1 (——) and computed static temperature $T_{st}$ normalised accordingly (– – –).

To assess the effect of convection in the lower part of the flow domain, we consider an equivalent one-dimensional purely diffusive process, i.e. in absence of fluid motion and hence convection. More specifically, we examine the static temperature distribution $T_{st}$ in an internally heated quiescent medium which has the same height as the porous medium and whose thermal conductivity equals the average conductivity of the lower part, i.e. $\langle \phi k + (1 - \phi ) k_{s} \rangle _{xz}$. Further, this medium is subject to exactly the same internal heating as the solid matrix, $\langle r \rangle _{xz}$. In this problem, we prescribe $T_{st}=\langle T_{int} \rangle$ at $y=h_{p}$ and zero-Neumann condition at $y=0$. The solution is a parabolic profile for $T_{st}$. In figure 12, we have included the plot of the normalised static temperature $\theta _{st,p}$, where the normalisation is performed according to (4.3a,b). We observe that the total variation of the static temperature in the purely diffusive process is four times larger than in Case 1; equivalently, the static temperature at the bottom is four times higher than in Case 1. We, therefore, deduce that convective heat transfer has a strong impact and reduces considerably the thermal gradients in the vertical direction, even though convective motions are significantly constrained inside the porous medium.

The relative strength of the various modes of heat transfer can be also assessed by comparing the different components of the Nusselt number, as provided in (3.4). Their plots are shown in figure 13. We observe that the local Nusselt number $Nu_y$ increases linearly with $y$ due to internal heating. Since the bottom wall is adiabatically isolated, the major component of $Nu_y$ in the near-wall region is the convective one. Also, both $Nu_{conv}$ and the diffusive component of the solid matrix $Nu_{diff,s}$ increase considerably away from the bottom wall. Actually, $Nu_{diff,s}$ has a peak close to the horizontal interface, and then it decreases rapidly and vanishes right at the horizontal interface. This is due to the prescribed thermal boundary condition (2.7) for the solid matrix which stipulates zero diffusive heat flux at points where the porosity becomes unity. However, the diffusive component of the fluid, $Nu_{diff}$, remains very small for most of the lower part of the domain and becomes noticeable only very close to the horizontal interface and at the point where $Nu_{diff,s}$ starts to drop. As a result, convection is the dominant mode throughout the lower part of the domain, in the sense that at any given $y$, $Nu_{conv}$ represents more than 50 % of $Nu_y$. The fact that $Nu_{diff}$ is much smaller than $Nu_{diff,s}$ in the bulk of the lower part of the domain is due to the difference in the thermal conductivities of the two phases.

Figure 13.

All contributions of the Nusselt number in the lower part of the domain normalised by the global Nusselt number in the pure fluid region $Nu_{f}$. The legend is total Nusselt number $Nu_{y}$ (——), convective contribution $Nu_{conv}$ ($\cdots \cdots$), diffusive contribution of the fluid $Nu_{ diff}$ ($-{\cdot }-{\cdot }-{\cdot }$) and of the solid $Nu_{diff,s}$ (– – –).

4.2. Comparisons of Cases 1 to 4

In this section, we make comparisons among Cases 1 to 4. As mentioned above, Cases 1–3 are in ascending order of area coverage $a$ of the porous medium and with uniform heating. However, in Case 4, the area coverage is the same as in Case 1 but the heat-load distribution is not uniform in the horizontal direction. Details of each case are summarised in table 2. It is also worth recalling that the total amount of heat added to the system is kept the same. By considering these cases, we examine the effects of the area coverage and the type of heat loading on the emerging flow structures and turbulence statistics.

4.2.1. Global flow properties

In figures 14, 15 and 16, we provide plots of time-averaged streamlines coloured by the time-averaged vertical velocity for Cases 2, 3 and 4, respectively. According to them, in Case 2 (whose area coverage is 88.4 %), the flow in the upper pure fluid region organises itself in a single-roll LSC aligned to a diagonal plane, similar to the LSC in classical RBC. However, a noticeable difference is that whereas in classical RBC there are two pairs of large counter-rotating vortices in the second diagonal plane, in Case 2, there is only one pair. Also, the observations made above for Case 1 regarding the descending fluid motion close to the corners and the circumvention of the porous medium are still valid for Case 2. Nonetheless, since the area coverage is much larger in Case 2 (88.4 % versus 64 %), these patterns are significantly weaker.

Figure 14.

Case 2: time-averaged streamlines coloured by non-dimensional time-averaged vertical velocity, $\langle v \rangle _{t}$. (a,b) Side view. (c) Bottom view.

Figure 15.

Case 3: time-averaged streamlines coloured by non-dimensional time-averaged vertical velocity, $\langle v \rangle _{t}$. (a,b) Side view. (c) Bottom view.

Figure 16.

Case 4: time-averaged streamlines coloured by non-dimensional time-averaged vertical velocity, $\langle v \rangle _{t}$. (a,b) Side view. (c) Bottom view.

According to figure 15, in Case 3, the flow in the pure fluid region also organises itself in an LSC. However, the structure of this LSC is different from Cases 1 and 2, as there is significant upward motion close to the side walls. Evidently, this is due to the area coverage (100 % in Case 3) which significantly reduces the descent of cold fluid towards the lower part of the domain. Also, the plots suggest that the vertical motion in the pure fluid region has a smaller magnitude compared to what we observed in Case 1. It is also worth mentioning that the streamlines and flow patterns in the pure fluid region, $y \geqslant h_{p}$, are very similar to those of the RBC-FS (not shown herein for the sake of economy of space) that was discussed in the previous section.

With regards to Case 4, from figure 16, we infer that the average streamlines in the pure fluid domain resemble those of Case 1. Also, as in Case 1, fluid descends from the corners and then circumvents the porous medium. Eventually, a portion of the fluid enters the porous medium from the side while the rest ascends along its side boundaries. However, a noticeable difference is that, in Case 4, there is no clear formation of convective rolls inside the porous medium; instead there is upward motion in the heated parts and weak downward motion in the non-heated parts.

In figure 17, we present colour plots of $\langle T_{int} \rangle _{t} / \langle T_{int} \rangle$ at $y = h_{p}$, i.e. of the time-averaged temperature at the plane of the horizontal interface, normalised by its mean (time-and-area-averaged) value. In this manner, we can assess the amplitude of the temperature variations at this plane. We readily infer that these variations are the smallest in Case 3. This is mainly due to the fact that the heating is uniform in the horizontal directions and the area coverage is 100 %. As the area coverage decreases, then the downward motion close to the corners becomes stronger, the heat-source term $\hat {r}$ gets higher and, therefore, the temperature variations increase. Accordingly, these variations are larger in Case 2 and even more so in Case 1. Also, as expected, the thermal gradients are even larger in Case 4 due to the non-uniformity of the heat loading.

Figure 17.

Colour plots of $\langle T_{int} \rangle _{t} / \langle T_{int} \rangle$ at the plane of the horizontal interface $y = h_{p}$. (a) Case 1. (b) Case 2. (c) Case 3. (d) Case 4.

With respect to figure 17, it is also interesting to observe that the high-temperature areas constitute in a sense the ‘footprints’ of the convective patterns that develop in the pure fluid region and inside the porous medium. For example, in figure 17(b), we can infer the alignment of the LSC with one diagonal plane. Also, the five ‘hot spots’ seen in this figure correspond to ascending plumes inside the porous medium. The fact that these hot spots are distributed symmetrically with respect to the diagonal of the LSC suggests an interaction between the plumes in the porous medium, the shear layer above the interface and the LSC. In addition, fluid coming out from these plumes eventually merges with the LSC in the pure fluid domain. Similar observations hold for the hot spots observed in figures 17(a) and 17(c). Concerning Case 4, the hot spots shown in figure 17(d) correspond to the heated parts of the porous medium.

For purposes of visualisation of the thermal field, in figure 18, we present contour plots of the time-averaged fluid temperature for the cases examined in this study. Therein, we can observe the principal plumes and infer the orientation of the LSC in the various cases.

Figure 18.

Iso-surfaces of the time-averaged fluid temperature throughout the flow domain. (a) Case 1. (b) Case 2. (c) Case 3. (d) Case 4.

Concerning the temperature variations in the vertical direction, our simulations predict that the global temperature difference, $\langle T_{bot} \rangle - T_{top}$, increases with the area coverage of the porous medium, even though the heat supply remains the same. In particular, it equals 15.9, 19.7 and 29 K for Cases 1, 2 and 3, respectively. This is attributed to the fact that large area coverage inhibits fluid motion inside and around the porous medium which, in turn, reduces thermal mixing since the heat supply comes from the solid matrix. Also, in Case 4, $\langle T_{bot} \rangle - T_{top}=12.4$ K, which is mildly lower than in Case 1. This is a manifestation of the fact that a non-uniform heat-load distribution promotes convective motions. The same trend is observed when comparing the temperature difference across the pure fluid region $\Delta \hat {T}_{f}$ and the Rayleigh number $Ra_{f}$. In other words, these quantities increase with the area coverage but decrease for non-uniform heat-loading. The values of the mean temperatures at the bottom, horizontal interface and top planes have been summarised in table 3 and those of the Rayleigh number in table 5.

Next, we compare our numerical predictions for the local Nusselt number. Profiles of $Nu_y$ across the entire domain for the four cases are presented in figure 19. As mentioned above, in all cases, $Nu_y$ increases linearly in the lower part of the domain due to the internal heating, while it remains constant in the pure fluid region and equal to $Nu_{f}$. Also, according to these plots, $Nu_y$ decreases substantially with the area coverage of the porous domain. This corroborates the fact that lower area coverage enhances convection. Moreover, the rate of decrease of $Nu_y$ becomes higher as the area coverage approaches 100 %. However, $Nu_y$ increases for non-uniform heat-load distributions, which is also due to the fact that such non-uniformities facilitate fluid recirculation in the pure fluid domain and hence convection.

Figure 19.

Profiles of the local Nusselt number $Nu_{y}$. The legend is Case 1 (——), Case 2 ($\cdots \cdots$), Case 3 (– – –) and Case 4 ($-{\cdot }-{\cdot }-{\cdot }$). The (red) horizontal line shows the location of the horizontal interface.

In summary, as the area coverage increases, $Nu_{f}$ and $Ra_{f}$ follow different trends; the former decreases while the latter increases. At a first glance, this appears to be counterintuitive because the correlation between $Nu_{f}$ and $Ra_{f}$ in classical RBC is the opposite. This difference is explained as follows. As the area coverage $a$ increases, there is less recirculation and convective motion, so the Nusselt number $Nu_{f}$ decreases. However, in the lower part of the domain, less recirculation reduces thermal mixing and leads to higher temperatures because the liquid stays longer in that part and, therefore, receives more heat from the solid matrix. Accordingly, as $a$ increases, then $\langle T_{int} \rangle$ gets higher and, therefore, the Rayleigh number $Ra_{f}$ gets higher too. Another way to interpret this is that, for the flows under consideration and contrary to classical RBC, it is not $\langle T_{int} \rangle$ but the heat supply from the solid matrix that is kept constant; this was also mentioned in § 3.1. Then, since $\langle T_{int} \rangle$ and hence $\Delta \hat {T}_{f}$ lie in the denominator of (3.4) for $Nu_{f}$ but in the numerator of (2.6) of $Ra_{f}$, it is natural for $Nu_{f}$ and $Ra_{f}$ to follow different trends as the area coverage increases. The same is also applied to explain the different trends of the Nusselt and Rayleigh numbers when the heat-load distribution becomes non-uniform.

The profiles of the various terms comprising $Nu_y$ in the lower part of the domain are presented in figure 20. Actually, in these plots, the terms have been normalised by the global Nusselt number of each case, $Nu_{f}$, to assess their relative contribution. Figure 20(a) shows plots of the convective term $Nu_{conv}/Nu_{f}$ for the various cases considered herein. According to them, in all cases, the convective contribution to the Nusselt number increases monotonically with $y$. Further, it decreases substantially with the area coverage but is increased with a non-uniform heat-load. For example, in Case 3 where $a=1$, the convective contribution to $Nu_{f}$ is less than 10 %, whereas in Case 1 ($a=0.64$), its contribution is approximately 50 %. Also, when the heat-load is non-uniform, the convective contribution rises even further, to 70 %. However, according to figure 20(b), the diffusive contributions for both solid and fluid follow the opposite trend, i.e. they increase with the area coverage but decrease with the non-uniform heat-load distribution.

Figure 20.

Profiles of the various terms comprising $Nu_{y}$ in the lower part of the domain, $y \leqslant h_{p}$, normalised by the global Nusselt number $Nu_{f}$: (a) $Nu_{conv} / Nu_{f}$; (b) $Nu_{diff} / Nu_{f}$ (grey lines) and $Nu_{diff,s} / Nu_{f}$ (black lines). The legend is Case 1 (——), Case 2 ($\cdots \cdots$), Case 3 (– – –) and Case 4 ($-{\cdot }-{\cdot }-{\cdot }$).

4.2.2. Flow statistics

In figure 21(a), we plot the mean (time-and-area-averaged) profiles of the normalised fluid temperature $\langle \theta \rangle$. One noticeable feature of these profiles is the strong temperature gradient in the vicinity of the horizontal interface which implies strong diffusive flux right above this interface; this is due to the fact that the fluid below the horizontal interface is heated while the fluid above is not. Also, in figure 21(a), the temperature profiles in the pure-fluid region shift to the right, i.e. towards $\langle \theta _{int} \rangle$, as the convective motions become stronger, either due to the reduction of the area coverage or non-uniform heat loading. We further observe that the overshoots of the temperature profile near the top and the horizontal interface become smaller as the convective motions decrease. In particular, these overshoots have disappeared in Case 3 wherein the temperature profile has become monotonically decreasing all across the pure-fluid region. This is in accordance with the previous observation that the overshoots are due to the recirculation of relatively cold fluid ascending from the side boundaries of the porous medium.

Figure 21.

(a) Time-and-area-averaged temperature $\langle \theta \rangle$ and (b) r.m.s. of the temperature fluctuations, $\theta _{rms}$. The inset is a zoom on the macroscopic interface. The legend is Case 1 (——), Case 2 ($\cdots \cdots$), Case 3 (– – –) and Case 4 ($-{\cdot }-{\cdot }-{\cdot }$). The (red) horizontal line denotes the plane of the horizontal interface.

With regards to the temperature distribution in the lower part, we observe that the temperature gradients and hence the diffusive flux increase as the area coverage becomes larger. As a result, the temperature at the interface $\langle \theta _{int} \rangle$ decreases with $a$.

The profiles of the r.m.s. of temperature fluctuations are plotted in figure 21(b). In the pure fluid region, and starting from the top, the profiles look similar to those of RBC-FS up to the vicinity of the horizontal interface. In particular, they attain a maximum at the edge of the thermal boundary layer at the free surface and decrease in the bulk. Overall, in the various cases, the values of $\langle \theta _{rms} \rangle$ in the bulk are close; nonetheless, they increase mildly as the area coverage becomes smaller and with the non-uniform heat loading. Then, in the vicinity of the horizontal interface, the r.m.s. profiles start to increase substantially. In Case 3, where $a=1$, the r.m.s. profile peaks a little above the interface and then decreases up to the bottom wall. However, in all other cases, the r.m.s. has a cusp right at the plane of the horizontal interface and continues to increase in the lower part of the domain. This difference is due to the fact that in Case 3, the area coverage is equal to unity and, therefore, fluid motion is sharply constrained in the lower part of the domain. In particular, in Case 3, fluid can enter or exit the lower part only by crossing the horizontal interface.

In the other cases, $\langle \theta _{rms} \rangle$ increases substantially with the depth in the lower part of the domain. This is due to the increase in temperature with the depth and the recirculation of cold fluid coming from the pure fluid region. Further, with regards to the temperature r.m.s. in the lower part, our simulations have predicted that they decrease as the area coverage decreases. This is explained by the fact that a larger porous region significantly limits fluid motion in the lower part of the domain as well as the descending motion of fluid coming from above. However, we observe that $\langle \theta _{rms} \rangle$ in the lower part of the domain is higher in Case 1 than in Case 4. We presume that this is attributed to the fact that in Case 1, the plumes in the central part of the porous medium do not remain fixed but slightly move from time to time, thereby increasing $\langle \theta _{rms} \rangle$; in contrast, such a motion is not observed in Case 4.

Finally, in figures 22(a) and 22(b), we have plotted the profiles of $v_{rms}$ and $u_{rms}$, respectively, for all cases considered in our study. Once again, in the pure fluid region, the $v_{rms}$ profiles are similar to those in RBC-FS up to the vicinity of the horizontal interface; in other words, they increase away from the free surface, attain a peak in the bulk and then decrease substantially as the horizontal interface is approached. In the lower part of the domain, the profiles also decrease with the depth of the domain until they vanish at the bottom wall. Overall, the $v_{rms}$ values become smaller as the size of the porous medium increases because vertical motion becomes more constrained. However, the $v_{rms}$ profiles for Cases 1 and 4 are quite similar throughout the flow domain, which means that the non-uniformity of the heat loading does not affect $v_{rms}$. The only noticeable difference is that in Case 4, the peak of $v_{rms}$ is located closer to the mid-plane of the pure fluid domain. This is because the difference between the strengths of the ascending and descending motions is less pronounced in Case 4 than in Case 1. Moreover, as mentioned above, the basic difference between the predicted $v_{rms}$ profiles and those in RBC-FS (Hay & Papalexandris Reference Hay and Papalexandris2019) is that now they do not vanish at the lower limit of the pure fluid domain due to the fluid motion below.

Figure 22.

Plots of the r.m.s. of velocity fluctuations in the whole domain: (a) $v_{rms}$ and (b) $u_{rms}$. The velocities have been made dimensionless using $\bar {\alpha }/\hat {h}_{tot}$. The legend is Case 3 (– – –), Case 1 (——), Case 2 ($\cdots \cdots$) and Case 4 ($-{\cdot }-{\cdot }-{\cdot }$). The (red) horizontal line denotes the plane of the horizontal interface.

According to figure 22(b), the $u_{rms}$ profiles are also similar to those in RBC-FS. More specifically, $u_{rms}$ is quite high at the top but decreases substantially with the depth. Subsequently, it starts to increase, peaks right above the horizontal interface and then decreases considerably as the interface is crossed. The large gradients observed in the vicinity of the horizontal interface are indicative of strong shear, the amplitude of which increases with the interface area. Then, in the lower part of the domain, $u_{rms}$ remains almost constant and drops to zero at the bottom wall. As with the other quantities, inside the porous medium, the amplitude of $u_{rms}$ drops significantly with the area coverage because flow recirculation becomes more restricted. Further, the $u_{rms}$ profiles between Cases 1 and 4 are overall similar. The main difference between these two profiles is the value at the free surface; in Case 1, the peak is higher than in Case 4 due to stronger ascending motion which in turn leads to stronger (horizontal) currents at the free surface.