Computational Modeling of Natural Convection in Nanofluid-Saturated Porous Media: An Investigation into Heat Transfer Phenomena

Abstract

A numerical study was carried out to analyze the phenomenon of natural convection in a porous medium saturated with nanofluid. In the study, the boundary element method was used for computational modeling. The fluid flow through a porous matrix is described using the Darcy–Brinkman–Forchheimer momentum equation. In addition, a mathematical model for nanofluids was used, which follows a single-phase approach and assumes that the nanoparticles within a fluid can be treated as an independent fluid with effective properties. A combination of single- and sub-domain boundary element methods was used to solve the relevant set of partial differential equations. The method was originally developed for pure flow scenarios, but also proves to be effective in the context of fluid flow through porous media. The results are calculated for the case of two- and three-dimensional square cavities. In addition to various values of dimensionless control parameters, including the porous Rayleigh number (Ra_p), Darcy number (Da), porosity (ϕ) and nanoparticle volume fractions (φ), the effects of the inclination angle of the cavity on the overall heat transfer (expressed by the Nusselt number (Nu)) and fluid flow characteristics were investigated. The results indicate a pronounced dependence of the overall heat transfer on the introduction of nanoparticles and inclination angle. The heat transfer in a two-dimensional cavity is increased for higher values of Darcy number in the conduction flow regime, while it is suppressed for lower values of Darcy number in the Darcy flow regime. In the case of a three-dimensional cavity, increasing the volume fraction of nanoparticles leads to a decrease in heat transfer, and furthermore, increasing the inclination angle of the cavity considerably weakens the buoyancy flow.

1. Introduction

A nanofluid is a colloidal suspension containing nanoparticles between 1 and 100 nm in size, including materials such as metals, oxides and carbides, dispersed in a base liquid such as water, oil or ethylene glycol. This breakthrough concept has been proven to improve the efficiency of cooling and heating processes in various industries and was originally introduced by Choi [1]. Recent efforts include both experimental and theoretical studies aimed at investigating the effects of improved heat transfer properties in various configurations and applications. A comprehensive overview of these studies can be found in the references [2,3,4]. Convective heat transfer plays a central role in various natural and engineering systems, including heating and cooling systems, geothermal systems and drying processes. The efficiency of heat transfer in these systems is considered crucial and can be significantly improved by the incorporation of nanoparticles, as described in the literature [5,6,7,8].

The thermal efficiency of several industrial applications, such as heat exchangers, can be improved by utilizing porous media (due to their large surface areas) and, furthermore, with the addition of nanoparticles. Most recent applications have been of microchannel heat sinks as innovative cooling devices, as introduced by Deng et al. [9] and Ghazvini and Shokouhmand [10]. Moreover, recent advances are comprehensively reviewed in the references [11,12,13,14,15], highlighting some recurring observations. In particular, the fusion of porous media and nanofluid leads to higher heat transfer rates, due to the presence of larger surface areas and more intense mixing. Nevertheless, there are still some unresolved issues, including the need for efficient numerical methods, exploring nanofluid flow in porous media under turbulent flow conditions, and conducting experimental studies involving different geometries and flow conditions.

Various mathematical models have been used in published studies to describe buoyancy-driven flow in porous media. Darcy’s law is the most commonly used mathematical model for the governing momentum equation, which is particularly valid in the laminar flow regime (when Reynolds number is $Re<10$), where viscous forces dominate over inertial forces at low velocities. In analogy to the Navier–Stokes equations, an extension of the governing momentum equations was established that includes the Brinkman term to account for viscous diffusion and the Forchheimer term to study the inertial effects in free convection [16].

There are two main approaches to the mathematical modelling of nanofluids: the single-phase and two-phase models. The single-phase approach assumes that nanoparticles behave similarly to water molecules and have similar local velocities. This assumption is particularly valid for low concentrations of nanoparticles ($2.5–5\%$) and for solid particles with a size between $1$ and $100 \mathrm{n}\mathrm{m}$ [17]. However, the two-phase model is better suited to describing mixtures of nanoparticles and base liquids in a physically correct way. This model includes mechanisms that consider the relative motion between the liquid and the nanoparticles, including Brownian diffusion and thermophoresis [18].

Various numerical methods have been used to investigate the complexity of heat transfer in porous media and to model the flow of fluids. These include the finite element method, the finite difference method, and the finite volume method, which are frequently used. Recently, the boundary element method (BEM) has emerged as an alternative, and is particularly popular for its efficiency in solving potential problems in fluid mechanics, such as inviscid flow and heat conduction. The main advantage of the BEM lies in its ability to solve partial differential equations by determining the boundary unknowns alone, bypassing the need to discretize the entire domain. However, this advantage is compromised if a suitable fundamental solution cannot be found, resulting in contributions of the domain remaining in the integral equation. In scenarios with inhomogeneous and non-linear problems such as diffusion–convection problems, the classical BEM is extended to treat domain integrals in addition to boundary integral equations. In particular, the evaluation of domain matrices becomes a central problem, since they are completely filled and asymmetric and require a lot of memory. Various strategies have been proposed to overcome this complication, including methods based on the expansion of the integral kernel [19], the dual reciprocity method [20] or the compression of the resulting complete matrices [21]. The numerical algorithm used in this study is divided into a single-domain and a sub-domain BEM, where the single-domain BEM deals with the kinematic aspect, while the sub-domain BEM is used to solve diffusion–advection type equations [22].

Previous studies have predominantly examined two-dimensional geometries in the context of convective flow within nanofluid-saturated porous media, with limited attention given to three-dimensional configurations. This research addresses this gap by presenting enhanced numerical simulations of convective flow within both two- and three-dimensional rectangular enclosures filled with a porous medium saturated with nanofluid. The study of three-dimensional geometry enables a more comprehensive understanding of fluid dynamics and heat transfer behavior in complex, real-world applications such as thermal management in electronic devices, enhanced geothermal systems and industrial heat exchangers. These simulations provide insights that can be used to design and optimize systems where precise control of heat transfer and fluid flow is critical. Furthermore, the study investigates the effect of the cavity’s inclination angle on flow behavior and heat transfer characteristics, offering a more comprehensive understanding of the dynamics involved in such systems. As a test case, a suspension of solid Cu nanoparticles in water is used as the base fluid. The calculations were performed using a BEM-based algorithm. The fluid flow in porous media is characterized by the Darcy–Brinkman–Forchheimer momentum equation. A single-phase nanofluid mathematical model is used since low concentrations of nanoparticles (2.5–5%) are considered. The effectiveness of the numerical code has been demonstrated in various applications, including pure fluid flow scenarios [22,23] and flow simulations in porous media [24,25]. In this study, the effects of inclination angle and nanoparticle volume fraction on convective heat transfer are systematically examined within both conduction-dominated and convection-dominated regimes. The combined influence of these parameters leads to complex interactions within the convective flow field. In a conduction-dominated regime, a low inclination angle can enhance the effects of the increased volume fraction of nanoparticles and maximize heat transfer. In contrast, in a convection-dominated regime at larger inclination angles, the effects of nanoparticles limit the enhancement of heat transfer. This nuanced understanding enables targeted optimization for specific applications where precise control of heat transfer is essential.

2. Mathematical Model

2.1. Equations Governing Flow in Porous Media

The mathematical model for the heat transfer of nanofluids in porous media is based on the conservation equations for mass, momentum and energy. These equations are derived from the Navier–Stokes equations, which are normally formulated at the microscopic level for pure fluid flows. However, due to the irregular and generally complex geometry of porous media, a microscopic description is impractical for modeling fluid flows. Therefore, all equations are averaged over the representative elementary volume (REV), taking into account only a part of the computational domain for the flow. The specific details of the averaging procedure described by Bear [26] are not discussed in this article.

The continuity equation expresses the conservation of mass for an incompressible fluid:

$$\stackrel{⃑}{\nabla }·\stackrel{⃑}{v}=0.$$

(1)

Here, the symbol $\stackrel{⃑}{v}$ represents the volume-averaged velocity vector.

The Brinkman–Forchheimer momentum equation is used in the present study, which reads as

$${\displaystyle \frac{1}{\varphi }}{\displaystyle \frac{\partial \stackrel{⃑}{v}}{\partial t}}+{\displaystyle \frac{1}{{\varphi }^2}}\left(\stackrel{⃑}{v}·\stackrel{⃑}{\nabla }\right)\stackrel{⃑}{v}={\displaystyle \frac{1}{{\rho }_{nf}}}\stackrel{⃑}{\nabla }p-{\beta }_{nf}\left(T-T_0\right)\stackrel{⃑}{g}+{\displaystyle \frac{1}{\varphi }}{\displaystyle \frac{{\mu }_{nf}}{{\rho }_{nf}}}{\nabla }^2\stackrel{⃑}{v}-{\displaystyle \frac{1}{K}}{\displaystyle \frac{{\mu }_{nf}}{{\rho }_{nf}}}\stackrel{⃑}{v}-{\displaystyle \frac{F\stackrel{⃑}{v}\left|\stackrel{⃑}{v}\right|}{K^{\frac{1}{2}}}}.$$

(2)

In the provided context, $\varphi$ represents porosity, $t$ signifies time, $p$ denotes pressure, $T$ stands for temperature, $\stackrel{⃑}{g}$ represents gravitational acceleration, $K$ denotes permeability, $F$ represents the Forchheimer coefficient, ${\rho }_{nf}$ represents the density of nanofluid, ${\beta }_{nf}$ represents the nanofluid thermal expansion coefficient, and ${\mu }_{nf}$ denotes the dynamic viscosity of nanofluid.

The Brinkman–Forchheimer momentum equation comprises two viscous and two inertia terms:

–

The Brinkman viscous term positioned third on the r.h.s., accounts for viscous forces and ensures compliance with non-slip boundary conditions along a boundary. It has similarities with the Laplacian term found in the Navier–Stokes equations formulated for pure fluid flow [16].

–

The Darcy term, positioned fourth on the r.h.s., is a linear term that establishes a connection between the velocity field and pressure difference. This relationship involves fluid viscosity and permeability ($K$), which is contingent upon the geometry of the porous medium and typically represents a second-order tensor. In the case of assuming an isotropic porous medium, the permeability becomes a scalar.

–

The Forchheimer term, positioned as the last term on the r.h.s., also referred to as the dimensionless form-drag constant, varies based on the characteristics of the porous medium and can be expressed using the Ergun model, as proposed in Nield and Bejan [16]:

$$K={\displaystyle \frac{{\varphi }^3{d}_p^2}{a{(1-\varphi )}^2}}, F={\displaystyle \frac{b}{\sqrt{a{\varphi }^3}}}.$$

(3)

Here, Ergun’s constants are denoted by $a$ and $b$, with specific values assigned as $a=150$ and $b=1.75$, as per [27]. Additionally, $d_p$ represents the average particle size of the bed.

Finally, the energy equation can be formulated as follows:

$$\sigma {\displaystyle \frac{\partial T}{\partial t}}+\left(\stackrel{⃑}{v}·\stackrel{⃑}{\nabla }\right)T={\displaystyle \frac{k_e}{{\left(\rho c_p\right)}_{nf}}}{\nabla }^{2}T,$$

(4)

here, $\sigma$ stands for the specific heat ratio, defined as $\sigma =\varphi +\left(1-\varphi \right){\left(\rho c_p\right)}_p/{\left(\rho c_p\right)}_{nf}$, where ${\left(\rho c_p\right)}_p$ and ${\left(\rho c_p\right)}_{nf}$ denote the heat capacities of the solid phase and the nanofluid phase, respectively. In addition, $k_e$ stands for the effective conductivity of the porous medium. Following the work of [28], it is assumed that the thermal properties of the solid matrix and the nanofluid are identical, which leads to $\sigma =1$ and $k_e=k_{nf}$.

2.2. Non-Dimensional Equations

To render Equations (1), (2) and (4) into non-dimensional form, the following dimensionless variables are utilized:

$$\stackrel{⃑}{v}\to {\displaystyle \frac{\stackrel{⃑}{v}}{v_0}}, \stackrel{⃑}{r}\to {\displaystyle \frac{\stackrel{⃑}{r}}{L}}, t\to {\displaystyle \frac{v_{0}t}{L}}, \stackrel{⃑}{g}\to {\displaystyle \frac{\stackrel{⃑}{g}}{g_0}}, p\to {\displaystyle \frac{p}{p_0}}, T\to {\displaystyle \frac{(T-T_0)}{∆T}}$$

(5)

In the given expressions, the parameters are defined as follows: the characteristic velocity ($v_0$) is expressed as $v_0=k_f/{(\rho c_p)}_{f}L$, where $k_f$ represents the fluid thermal conductivity, ${(\rho c_p)}_f$ denotes the heat capacity for the fluid phase, and $L$ is the characteristic length. Additionally, the characteristic temperature ($T_0$) is determined as $T_0=(T_2-T_1)/2$, and the characteristic temperature difference ($\Delta T$) is given by $\Delta T=(T_2-T_1)/2.$ The characteristic pressure ($p_0$) is established as $p_0=1 \mathrm{b}\mathrm{a}\mathrm{r}$, and the gravitational acceleration is denoted as $g_0=9.81 \mathrm{m}/{\mathrm{s}}^2$.

2.3. Velocity–Vorticity Formulation

The governing equations are reformulated through the introduction of the velocity–vorticity formulation, effectively dividing the computational scheme into kinematic and kinetic components. The vorticity vector is defined as the curl of the velocity, $\stackrel{⃑}{\omega }=\stackrel{⃑}{\nabla }\times \stackrel{⃑}{v}$, ensuring that both the velocity and vorticity fields adhere to the solenoidal condition, $\stackrel{⃑}{\nabla }·\stackrel{⃑}{v}=0$, $\stackrel{⃑}{\nabla }·\stackrel{⃑}{\omega }=0$. The kinematics equation, derived from the mass conservation law (Equation (1)), is a vector elliptic partial differential equation of a Poisson type and can be expressed as follows:

$${\nabla }^2\stackrel{⃑}{v}+\stackrel{⃑}{\nabla }\times \stackrel{⃑}{\omega }=0.$$

(6)

Applying the curl operator to the momentum equation (Equation (2)) allows for the derivation of the vorticity transport equation which represents the kinetic computational aspect along with the energy transport equation:

$$\begin{gathered} \left(\stackrel{⃑}{v}·\stackrel{⃑}{\nabla }\right)\stackrel{⃑}{\omega }=\left(\stackrel{⃑}{\omega }·\stackrel{⃑}{\nabla }\right)\stackrel{⃑}{v}-C_A\mathrm{Pr}{Ra}_T {\varphi }^2 \stackrel{⃑}{\nabla }\times T \stackrel{⃑}{g}+{\mathrm{C}}_B\mathrm{Pr}\varphi {\nabla }^2\stackrel{⃑}{\omega }- \\ {C}_B{\displaystyle \frac{Pr}{Da}}{\varphi }^2\stackrel{⃑}{\omega }-{\displaystyle \frac{F}{Da}}{\varphi }^2\left|\stackrel{⃑}{v}\right|\stackrel{⃑}{\omega }, \end{gathered}$$

(7)

$$\left(\stackrel{⃑}{v}·\stackrel{⃑}{\nabla }\right)T=C_C {\nabla }^{2}T.$$

(8)

The independent non-dimensional parameters featured in the momentum equation are as follows:

–

The fluid Rayleigh number:

$${Ra}_T={\displaystyle \frac{g {\beta }_T \Delta T L^3 {\rho }_f {\left(\rho c_p\right)}_f}{({\mu }_f k_f)}},$$

(9)

–

The Prandtl number:

$$Pr={\displaystyle \frac{{\mu }_f c_p}{k_f}},$$

(10)

–

The Darcy number:

$$Da={\displaystyle \frac{K}{L^2}}.$$

(11)

Furthermore, the porous Rayleigh number can be defined in dependance of the fluid Rayleigh number ${Ra}_T$ and Darcy number $Da$ as:

$${Ra}_P={Ra}_T·Da.$$

(12)

In the equations provided above, the parameters $C_A$, $C_B$, and $C_C$ represent the properties of the nanofluid and are defined by the following expressions:

$$C_A={\displaystyle \frac{{\mu }_{nf}}{{\mu }_f}}{\displaystyle \frac{{\rho }_f}{{\rho }_{nf}}}, C_B={\displaystyle \frac{{\beta }_{nf}}{{\beta }_f}}, C_C={\displaystyle \frac{{\alpha }_{nf}}{{\alpha }_f}},$$

(13)

here ${\alpha }_{nf}$ represents the thermal diffusivity of the nanofluid, defined as ${\alpha }_{nf}=k_{nf}/{\left(\rho c_p\right)}_{nf}$ where ${\alpha }_f$ is the thermal diffusivity of the pure fluid, given by ${\alpha }_f=k_f/{\left(\rho c_p\right)}_f.$ The nanofluid properties are determined using the expressions outlined in the next section. For simulating pure fluid flow, the parameters are set as $C_A=C_B=C_C=1$. As this paper focuses solely on steady flow simulations, the vorticity and energy transport equations are considered without time derivatives $\partial \stackrel{⃑}{\omega }/\partial t=\partial T/\partial t=0$.

The general momentum equation (Equation (2)) contains the pressure term in gradient form, which can lead to numerical instabilities. In the velocity–vorticity formulation, the pressure term is removed from the momentum equation as a primary variable, resulting from the application of the curl operator.

The determined partial differential Equations (6)–(8) form a non-linear system of equations that controls the unknown fields of velocity, vorticity and temperature. The characterization of the heat and mass transfer in the area of the porous media domain is clearly determined by the specification of the buoyancy coefficient, the Rayleigh, Prandtl, Lewis and Darcy numbers as well as the relevant initial and boundary conditions.

2.4. Nanofluid Properties

Nanofluid properties are expressed in terms of the relationships between the properties of the pure fluid and the pure solid. In all subsequent expressions, the indices $f$ and $s$ denote the fluid and solid phases, respectively.

Initially, the nanofluid’s solid volume fraction ($\phi$) is established as the proportion of the solid particles’ volume ($V_s$) to the combined volume of solid particles and fluid ($V_s+V_f$):

$$\phi ={\displaystyle \frac{V_s}{V_s+V_f}}.$$

(14)

Different models are used to describe the relationships between nanofluid and pure fluid properties, and a thorough investigation of different models can be found in [29]. For this study, the nanoparticles are assumed to have a spherical shape and all the assumed models apply to scenarios characterized by small temperature gradients.

The expressions for the density (${\rho }_{nf}$), the effective dynamic viscosity ${\mu }_{nf}$, the heat capacity of nanofluid ${\left(\rho c_p\right)}_{nf}$ and the thermal expansion coefficient ${\left(\rho \beta \right)}_{nf}$ are as follows:

$${\rho }_{nf}=\left(1-\phi \right){\rho }_f+\phi {\rho }_s,$$

(15)

$${\mu }_{nf}={\displaystyle \frac{{\mu }_f}{{\left(1-\phi \right)}^{2.5}}},$$

(16)

$${\left(\rho c_p\right)}_{nf}=\left(1-\phi \right){\left(\rho c_p\right)}_f+\phi {\left(\rho c_p\right)}_s,$$

(17)

$${\left(\rho \beta \right)}_{nf}=\left(1-\phi \right){\left(\rho \beta \right)}_f+\phi {\left(\rho \beta \right)}_s.$$

(18)

The Wasp model [30] provides the expression for effective thermal conductivity ($k_{nf}$) as follows:

$$k_{nf}=k_f{\displaystyle \frac{k_s+2k_f-2\phi (k_f-k_s)}{k_s+2k_f+\phi (k_f-k_s)}}.$$

(19)

Additional assumptions for the model used include that the nanoparticles are in thermal equilibrium with the base fluid and that the boundary condition of slip resistance is met. The fluid flow is assumed to be laminar, stationary, Newtonian and incompressible. The relationship between density and temperature can be described using the Boussinesq approximation as follows:

$${\rho }_{nf}={\rho }_0\left(1-{\beta }_{nf}\left(T-T_0\right)\right).$$

(20)

here, the subscript 0 denotes a reference state.

3. The Boundary Element Method

The BEM-based algorithm is used to solve the relevant non-linear partial differential equations (Equations (6)–(8)). Although the classical boundary element method (BEM) has been extended to incorporate domain integrals, several key advantages persist, justifying its use in the present study. The extended BEM can significantly reduce computational demands compared to fully volume-based methods, as it requires discretization primarily on boundaries rather than throughout the domain. This boundary-centric approach delivers high accuracy in scenarios where boundary interactions are essential. Furthermore, even though domain integrals in non-linear problems produce fully populated matrices, the BEM’s boundary-focused framework generally results in lower memory usage, especially in sparse or boundary-dominant cases. Overall, the extended BEM provides a versatile and efficient computational framework, particularly advantageous in boundary-dominated or unbounded domain problems, as it reduces computational overheads while enhancing precision at critical boundary regions.

In order to determine exact values for the boundary vortices, the algorithm is divided into components for a single domain and for sub-domains. The kinematics equation (Equation (6)) is solved with the single-domain BEM, which provides the values for the boundary vorticity. The sub-domain BEM then solves the vorticity and energy transport equations (Equations (7) and (8)) to determine the unknown values for vorticity and temperature. The algorithm was originally developed for pure fluid flow simulations [22] and subsequently adapted for nanofluids [23] and flow simulations in porous media [24,25].

The computational approach leads to a completely filled system of equations, which limits the maximum grid size due to memory constraints. To mitigate this drawback, the fast BEM is used, which utilizes sparse approximations of the full matrices [31]. The main advantage of using a single-domain BEM for boundary vortex values is the ability of the algorithm to preserve the mass in complicated geometries—a property that is not present when using velocity derivatives to compute boundary vortex values.

The numerical algorithm is structured as follows: At the beginning, the boundary conditions, which can be of Dirichlet or/and Neumann type, for velocity and temperature are specified. These conditions are used to solve the kinematic equation (Equation (6)) for the velocity values of the domain and the energy equation (Equation (8)) for the temperature values of the domain. In addition, temperature and temperature flux conditions are defined on the solid walls, together with the imposition of no-slip boundary conditions. The initial boundary conditions for vorticity values are unknown and are determined as part of the algorithm using the single-domain BEM from the kinematics equation (Equation (6)) [32]. The vorticity values within the domain are calculated by employing a sub-domain boundary element method (BEM) applied to the vorticity transport equation (Equation (7)). The overall structure of the numerical algorithm is as follows:

–

Determination of fluid and porous media properties.

–

Calculation of vorticity values on the boundary using single-domain BEM from the kinematics equation (Equation (6)).

–

Computation of velocity values within the domain using the sub-domain BEM from the kinematics equation (Equation (6)).

–

Determination of temperature values within the domain using the sub-domain BEM from the energy equation (Equation (8)).

–

Determination of vorticity values within the domain using the sub-domain BEM from the vorticity equation (Equation (7)).

Checking of convergence—repeat steps 2 to 5 until all flow fields achieve the required accuracy.

3.1. Integral Formulation of Governing Equations

All governing equations are expressed in integral form through the application of Green’s second identity. This involves the use of the fundamental solution $u^{\ast }$ of the diffusion operator for the two unknown field functions:

$$u^{\ast }={\displaystyle \frac{1}{4\pi \left|\stackrel{⃑}{\xi }-\stackrel{⃑}{r}\right|}},$$

(21)

here $\xi$ is a source or collocation point on the boundary $\mathsf{\Gamma }$ and $\stackrel{⃑}{r}$ integration point in the domain $\Omega$. The unknown boundary vorticity values are acquired by applying the single-domain BEM to the kinematics equation (Equation (6)), which needs to be expressed in its tangential form:

$$\begin{gathered} c\left(\stackrel{⃑}{\xi }\right)\stackrel{⃑}{n}\left(\stackrel{⃑}{\xi }\right)\times \stackrel{⃑}{v}\left(\stackrel{⃑}{\xi }\right)+\stackrel{⃑}{n}\left(\stackrel{⃑}{\xi }\right)\times {\int }_{\mathsf{\Gamma }}\stackrel{⃑}{v}\stackrel{⃑}{\nabla }{u}^{\ast }·\stackrel{⃑}{n}d\mathsf{\Gamma }= \\ \stackrel{⃑}{n}\left(\stackrel{⃑}{\xi }\right)\times {\int }_{\mathsf{\Gamma }}\stackrel{⃑}{v}\times (\stackrel{⃑}{n}\times \stackrel{⃑}{\nabla })u^{\ast }d\mathsf{\Gamma }+\stackrel{⃑}{n}\left(\stackrel{⃑}{\xi }\right)\times {\int }_{\mathsf{\Omega }}(\stackrel{⃑}{\omega }\times \stackrel{⃑}{\nabla }{u}^{\ast })d\mathsf{\Omega }, \end{gathered}$$

(22)

Here, $\mathsf{\Omega }$ represents the computational domain, and $\mathsf{\Gamma }=\partial \mathsf{\Omega }$ is the boundary of the domain. The geometric factor $c\left(\stackrel{⃑}{\xi }\right)$ is defined as $c\left(\stackrel{⃑}{\xi }\right)=\theta /4\pi$, where $\theta$ is the inner angle with the origin at $\stackrel{⃑}{\xi }$. If $\stackrel{⃑}{\xi }$ lies inside the domain, then $c\left(\stackrel{⃑}{\xi }\right)=1$; if $\stackrel{⃑}{\xi }$ lies on a smooth boundary, then $c\left(\stackrel{⃑}{\xi }\right)=1/2$. Additionally, $\stackrel{⃑}{n}$ is a vector normal to the boundary.

3.2. Kinematics Equation

The sub-domain BEM is applied to the kinematics equation (Equation (6)) to compute the velocity values within the domain. The integral equation is given by

$$c\left(\xi \right)\stackrel{⃑}{v}\left(\xi \right)+{\int }_{\mathsf{\Gamma }}\stackrel{⃑}{v}\left(\stackrel{⃑}{n}·\stackrel{⃑}{\nabla }\right)u^{\ast }d\mathsf{\Gamma }={\int }_{\mathsf{\Gamma }}\stackrel{⃑}{v}\times \left(\stackrel{⃑}{n}\times \stackrel{⃑}{\nabla }\right)u^{\ast }d\mathsf{\Gamma }+{\int }_{\mathsf{\Omega }}\left(\stackrel{⃑}{\omega }\times \stackrel{⃑}{\nabla }{u}^{\ast }\right)d\mathsf{\Omega }.$$

(23)

The main advantage of these formulations is that the resulting integral equation does not contain any derivatives of the velocity or vorticity fields. This property makes it possible to set the source point exclusively at the function nodes. The calculation of the values for the range velocity is based on the known boundary values of the velocity from the initial boundary conditions, where the values for the range and the boundary of the vorticity are known from the previous iteration.

3.3. Vorticity and Energy Equations

To obtain the integral representation of the equations for vorticity and energy, we use the same fundamental solution of Laplace’s equation as mentioned above. The final integral form of the vorticity transport equation is given by

$$\begin{gathered} c\left(\stackrel{⃑}{\xi }\right){\omega }_j\left(\stackrel{⃑}{\xi }\right)+{\int }_{\mathsf{\Gamma }}{\omega }_j\stackrel{⃑}{\nabla }{u}^{\ast }·\stackrel{⃑}{n} d\mathsf{\Gamma }={\int }_{\mathsf{\Gamma }}{u}^{\ast }{q}_j d\mathsf{\Gamma } \\ +{\displaystyle \frac{1}{Pr}}{\displaystyle \frac{1}{C_B}}{\displaystyle \frac{1}{\varphi }}{\int }_{\mathsf{\Gamma }}\stackrel{⃑}{n}·\left\{u^{\ast }\left(\stackrel{⃑}{v}{\omega }_j-\stackrel{⃑}{\omega }{v}_j\right)\right\} d\mathsf{\Gamma }-{\displaystyle \frac{1}{Pr}}{\displaystyle \frac{1}{C_B}}{\displaystyle \frac{1}{\varphi }}{\int }_{\mathsf{\Omega }}\left(\stackrel{⃑}{v}{\omega }_j-\stackrel{⃑}{\omega }{v}_j\right)·\stackrel{⃑}{\nabla }{u}^{\ast } d\mathsf{\Omega } \\ -{Ra}_T{\displaystyle \frac{C_A}{C_B}} \varphi {\int }_{\mathsf{\Gamma }}{\left(u^{\ast }T\stackrel{⃑}{g}\times \stackrel{⃑}{n}\right)}_j d\mathsf{\Gamma }-{Ra}_T{\displaystyle \frac{C_A}{C_B}} \varphi {\int }_{\mathsf{\Gamma }}{\left(T\stackrel{⃑}{\nabla }\times u^{\ast }\stackrel{⃑}{g}\right)}_j d\mathsf{\Omega } \\ +{\displaystyle \frac{1}{Da}}\varphi {\int }_{\mathsf{\Omega }}{\omega }_j{u}^{\ast }d\mathsf{\Omega }+{\displaystyle \frac{F}{\mathrm{Pr}\sqrt{Da}}}{\displaystyle \frac{1}{C_B}}\varphi \left|\stackrel{⃑}{v}\right|{\int }_{\mathsf{\Omega }}{\omega }_j{u}^{\ast }d\mathsf{\Omega }. \end{gathered}$$

(24)

Finally, the integral representation of the energy transport equation is expressed as

$$\begin{array}{c}c\left(\stackrel{⃑}{\xi }\right)T\left(\stackrel{⃑}{\xi }\right)+{\int }_{\mathsf{\Gamma }}T\stackrel{⃑}{\nabla }{u}^{\ast }·\stackrel{⃑}{n} d\mathsf{\Gamma }\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} ={\int }_{\mathsf{\Gamma }}{u}^{\ast }{q}_T d\mathsf{\Gamma }+{\displaystyle \frac{1}{C_C}}\left[{\int }_{\mathsf{\Gamma }}\stackrel{⃑}{n}·\left\{u^{\ast }\left(\stackrel{⃑}{v}T\right)\right\} d\mathsf{\Gamma }-{\int }_{\mathsf{\Omega }}\left(\stackrel{⃑}{v}T\right)·\stackrel{⃑}{\nabla }{u}^{\ast } d\mathsf{\Omega }\right].\end{array}$$

(25)

In the given equations, $q_j$ denotes a component of vorticity flux, and $q_T$ represents the heat flux.

Within the sub-domain BEM method, a mesh is created for the entire domain $\mathsf{\Omega }$, where each mesh element is referred to as a sub-domain. Equations are formulated for each of these sub-domains. Shape functions are employed to interpolate the field functions, as well as the flux across the boundary and within the domain. In this study, hexahedral sub-domains comprising 27 nodes were employed, allowing for continuous quadratic interpolation of field functions. The boundary of each hexahedron is composed of six boundary elements. Flux interpolation on each boundary element is carried out using a discontinuous linear interpolation scheme with four nodes. The use of discontinuous interpolation mitigates definition problems in corners and edges. The resultant discrete system of equations is over-determined and is addressed through a least squares solver [22].

4. Test Examples

The examined scenario involves examples of two- and three-dimensional cavities (Figure 1), filled with a porous medium and entirely saturated with nanofluid. The horizontal walls are treated as adiabatic, while the vertical walls are differentially heated with constant temperatures on opposite walls. Owing to the temperature contrast along the vertical walls, variations in fluid density occur, giving rise to buoyancy forces that induce convective motion. The fluid ascends along the hot wall, initiating heat transport toward the cold wall. The magnitude of heat flux is contingent upon the fluid type, the nature and quantity of the added nanoparticles, and the permeability of the porous medium. The assumptions made include the non-deformability, isotropy and homogeneity of the porous medium as well as the additional assumption that there is no heat transfer between the solid and liquid phases.

The wall heat flux is evaluated to quantify the overall heat transfer of nanofluids through porous media, which is characterized by the average Nusselt number. For nanofluids, the Nusselt number can be expressed as

$$Nu={\displaystyle \frac{k_{nf}}{k_f}}{\int }_{\mathsf{\Gamma }}\stackrel{⃑}{\nabla }T·\stackrel{⃑}{n} d\mathsf{\Gamma }$$

(26)

where $\mathsf{\Gamma }$ represents the surface through which the heat flux is computed, $\stackrel{⃑}{n}$ is the unit normal vector to this surface, $k_{nf}$ is the thermal conductivity of the nanofluid, and $k_f$ is the thermal conductivity of the base fluid.

The Cu nanoparticles were considered to be added to the water as a base fluid. The thermophysical properties of both are given in

Table 1. Thermophysical properties of water and Cu nanoparticles, adapted from [25].

	cp [J/kg K]	ρ [kg/m3]	K [W/m K]	β [×10−5 K−1]	α [×10−7 m2/s]
Water	4179	997.1	0.613	21	1.47
Cu	385	8933	400	1.67	1163

.

Validation Tests

To achieve a grid-independent solution, a grid sensitivity analysis was first carried out. Four different non-uniform grids (21 × 21, 41 × 41, 61 × 61 and 81 × 81) were tested for the Cu–water nanofluid and $Pr=6.2$, $Da={10}^{-2}, { 10}^{-4}$ and ${ 10}^{-6}$, respectively, ${Ra}_p=1000$, $\varphi =0.4$ and $\phi =0.05$, which are listed in

Table 2. Grid sensitivity analysis for 2D geometry; comparison of average Nusselt number for ${Ra}_p=1000$, $Pr=6.2$, $Da={10}^{-2}$, ${10}^{-4}$, and ${10}^{-6}$, $\varphi =0.4$, and $\phi =0.05$.

${\mathit{R}\mathit{a}}_{\mathit{p}}=1000$	${\mathit{R}\mathit{a}}_{\mathit{T}}={10}^5$	${\mathit{R}\mathit{a}}_{\mathit{T}}={10}^7$	${\mathit{R}\mathit{a}}_{\mathit{T}}={10}^9$
Nonuniform Gird	$\mathit{D}\mathit{a}={10}^{-2}$	$\mathit{D}\mathit{a}={10}^{-4}$	$\mathit{D}\mathit{a}={10}^{-6}$
21 × 21	3.416	9.076	13.185
41 × 41	3.400	9.132	12.911
61 × 61	3.401	9.132	12.992
81 × 81	3.400	9.131	12.991

together with the resulting average Nusselt numbers. Based on the grid independence test, the non-uniform 41 × 41 grid was selected for further analysis.

In the case of the 3D geometry example, four non-uniform grids (12 × 12 × 12, 20 × 8 × 20, 22 × 10 × 22, and 30 × 10 × 30) for the Cu–water nanofluid and $Pr=6.2$, $Da={10}^{-2}-{10}^{-5}$, ${Ra}_p=1000$, $\varphi =0.4$ and $\phi =0.0$ were tested. Based on the results presented in

Table 3. Grid sensitivity analysis for 3D geometry; comparison of average Nusselt number for Ra_p = 100, Pr = 6.2, ϕ = 0.4, φ = 0.0 and various Da.

${\mathit{R}\mathit{a}}_{\mathit{p}}=100$	${\mathit{R}\mathit{a}}_{\mathit{T}}={10}^3$	${\mathit{R}\mathit{a}}_{\mathit{T}}={10}^4$	${\mathit{R}\mathit{a}}_{\mathit{T}}={10}^5$	${\mathit{R}\mathit{a}}_{\mathit{T}}={10}^6$	${\mathit{R}\mathit{a}}_{\mathit{T}}={10}^7$
Grid	$\mathit{D}\mathit{a}={10}^{-1}$	$\mathit{D}\mathit{a}={10}^{-2}$	$\mathit{D}\mathit{a}={10}^{-3}$	$\mathit{D}\mathit{a}={10}^{-4}$	$\mathit{D}\mathit{a}={10}^{-5}$
12 × 12 × 12	1.0423	1.5428	2.3432	2.9784	3.3008
20 × 8 × 20	1.0394	1.5329	2.3313	2.9575	3.2950
22 × 10 × 22	1.0393	1.5327	2.3307	2.9552	3.2945
30 × 10 × 30	1.0393	1.5325	2.3303	2.9541	3.2934

, the 20 × 8 × 20 grid, with 28,577 nodes, demonstrated acceptable accuracy and was selected for further analysis.

The numerical code was validated using several test cases with different geometries and control parameters. A subset of the results is shown in

Table 4. Comparison of the average Nusselt number (Nu) for natural convection in porous media saturated with pure fluid ($Pr=1.0$) across various governing parameters with the data from [28,33,34].

$\mathit{D}\mathit{a}$	${\mathit{R}\mathit{a}}_{\mathit{p}}$	${\mathit{R}\mathit{a}}_{\mathit{T}}$	[28]	[33]	[34]	Present
			$\varphi =0.4$
	10	103	1.005	1.008	-	1.008
10−2	100	104	1.404	1.359	-	1.371
	1000	105	3.159	2.986	-	3.049
	10	105	1.064	1.062	-	1.067
10−4	100	106	2.580	2.702	-	2.671
	1000	107	7.677	8.903	-	8.377
	10	107	1.074	1.072	1.07	1.092
10−6	100	108	2.969	2.975	3.07	3.224
	1000	109	11.699	11.892	12.80	12.519
			$\varphi =0.6$
	10	103	1.012	1.010	-	1.012
10−2	100	104	1.489	1.533	-	1.503
	1000	105	3.430	3.602	-	3.499
	10	105	1.066	1.065	-	1.070
10−4	100	106	2.686	2.764	-	2.777
	1000	107	8.452	9.454	-	9.174
	10	107	1.074	1.072	-	1.093
10−6	100	108	2.982	2.980	-	3.241
	1000	109	12.098	11.924	-	12.895
			$\varphi =0.9$
	10	103	-	1.015	-	1.017
10−2	100	104	-	1.667	-	1.643
	1000	105	-	4.075	-	3.980
	10	105	-	1.066	-	1.073
10−4	100	106	-	2.817	-	2.867
	1000	107	-	9.947	-	9.917
	10	107	-	1.072	1.07	1.093
10−6	100	108	-	2.986	3.09	3.252
	1000	109	-	12.069	13.29	13.164

and

Table 5. Comparison of the average Nusselt number (Nu) for natural convection in porous media saturated with nanofluid ($Pr=6.2$) across various governing parameters with the data from [28].

${\mathit{R}\mathit{a}}_{\mathit{p}}$	$\mathit{D}\mathit{a}$	${\mathit{R}\mathit{a}}_{\mathit{T}}$	[28]	Present	[28]	Present	[28]	Present
$\phi =0.05$
			$\varphi =0.4$		$\varphi =0.6$		$\varphi =0.8$
1000	10−2	105	3.433	3.400	3.850	3.826	4.162	4.145
1000	10−4	107	9.117	9.132	9.590	9.743	9.901	10.154
1000	10−6	109	11.778	12.991	11.899	13.128	11.976	13.195
$\varphi =0.4$
			$\phi =0.0$		$\phi =0.025$		$\phi =0.05$
10	10−2	103	1.007	1.008	1.081	1.083	1.160	1.162
1000	10−2	105	3.302	3.282	3.370	3.345	3.433	3.400
1000	10−6	109	11.867	13.238	11.847	13.131	11.778	12.991

, in which the Nusselt number values (Nu) for different parameters are compared with values stated in the literature [28,33,34]. The high agreement between the calculated results and the published data confirms the accuracy of the developed numerical algorithm and it can be used for further calculations.

5. Results and Discussion

The heat transfer and flow properties were further investigated, focusing on the effects of the geometry, the inclination angle, the volume fraction of the nanoparticles and other governing parameters.

Table 6. Nusselt number values for the cases of 2D and 3D geometry, ${Ra}_p=1000$, $Pr=6.2$, $\varphi =0.4$, various $Da$ numbers, volume fraction of nanoparticles $\phi$ and inclination angle $\alpha$.

${\mathit{R}\mathit{a}}_{\mathit{p}}=1000$$, \mathit{P}\mathit{r}=6.2$$, \mathit{\varphi }=0.4$
$\alpha =0°$
${Ra}_T={10}^5$, $Da={10}^{-2}$		2D	3D
φ	0.0	3.282	3.072
	0.025	3.345	2.977
	0.05	3.400	2.884
${Ra}_T={10}^7$, $Da={10}^{-4}$
φ	0.0	9.072	8.528
	0.025	9.115	8.097
	0.05	9.132	7.727
${Ra}_T={10}^9$, $Da={10}^{-6}$
φ	0.0	13.238	14.270
	0.025	13.131	13.664
	0.05	12.991	13.078
α = 15°
${Ra}_T={10}^5$, $Da={10}^{-2}$		2D	3D
φ	0.0	2.847	2.653
	0.025	2.909	2.574
	0.05	2.966	2.497
${Ra}_T={10}^7$, $Da={10}^{-4}$
φ	0.0	7.466	6.935
	0.025	7.494	6.176
	0.05	7.503	6.313
${Ra}_T={10}^9$, $Da={10}^{-6}$
φ	0.0	10.726	11.684
	0.025	10.626	11.185
	0.05	10.496	10.686
α = 30°
${Ra}_T={10}^5$, $Da={10}^{-2}$		2D	3D
φ	0.0	2.266	2.123
	0.025	2.335	2.070
	0.05	2.401	2.017
${Ra}_T={10}^7$, $Da={10}^{-4}$
φ	0.0	5.137	4.749
	0.025	5.175	4.537
	0.05	5.206	4.337
${Ra}_T={10}^9$, $Da={10}^{-6}$
φ	0.0	7.143	7.637
	0.025	7.076	7.304
	0.05	6.998	6.980
α = 60°
${Ra}_T={10}^5$, $Da={10}^{-2}$		2D	3D
φ	0.0	1.276	1.249
	0.025	1.269	1.240
	0.05	1.261	1.231
${Ra}_T={10}^7$, $Da={10}^{-4}$
φ	0.0	1.522	1.474
	0.025	1.501	1.462
	0.05	1.487	1.451
${Ra}_T={10}^9$, $Da={10}^{-6}$
φ	0.0	1.568	1.563
	0.025	1.554	1.557
	0.05	1.513	1.548

shows the results of the average Nu number for ${Ra}_p=1000$, $Pr=6.2$, $\varphi =0.4$ and different values of Darcy number, nanoparticle volume fraction φ and inclination angle α for two- and three-dimensional geometry.

Furthermore, Figure 2 and Figure 3 present the three-dimensional temperature and velocity fields for various values of $Da$ and ${Ra}_p$ at $Pr=6.2$, $\phi =0.05$, and $\varphi =0.4$. Figure 4 presents isotherms and streamlines with maximum values of stream function $\left|{\psi }_{max}\right|$ for ${Ra}_p=1000$, $Pr=6.2$, $Da={10}^{-6}$, and $\phi =0.05$, at inclination angles of $0°$, $15°$, $30°$, and $60°$ at the horizontal midsection of the 3D cavity. Figure 5 and Figure 6 also show graphically the dependence of the angle of inclination and the volume fraction of the nanoparticles on the average heat transfer.

The results in Table 6 demonstrate that, within 2D geometry, the overall heat transfer increases with higher nanoparticle volume fractions at low inclination angles (α < 30°) and in a conduction-dominated regime ($Da\ge {10}^{-4}$). Conversely, in 3D geometry, changes in nanoparticle volume fraction do not significantly influence heat transfer under comparable conditions.

From the presented temperature and velocity fields it can be observed that at low Rayleigh numbers (${Ra}_p=10$, ${Ra}_T={10}^7$) heat transfer is primarily governed by conduction, resulting in a nearly linear temperature distribution across the enclosure or cavity. Temperature variations exhibit a steady gradient from the hot wall to the cold wall, with minimal deviation throughout the domain. In this regime, the fluid remains predominantly stationary, leading to a stable velocity field characterized by near-zero or very low velocities, primarily consisting of slight movements near the boundaries.

Conversely, at higher Rayleigh numbers (${Ra}_p>100, {Ra}_T={10}^8$), the convective motion becomes more pronounced, leading to a temperature field characterized by non-linear gradients, well-defined boundary layers, and curvature of isotherms that follow the flow patterns. This regime also exhibits potential temperature stratification, where warmer fluid rises and cooler fluid descends, resulting in active convective currents throughout the domain. The velocity field in this case reflects the dynamic movement of the fluid, illustrating the influence of buoyancy-driven convection on heat transfer processes within the enclosure.

Figure 3 illustrates the influence of the Darcy number on the temperature and velocity fields at a fixed porous Rayleigh number. In all three cases, the high Rayleigh number indicates that convection is the dominant heat transfer mechanism, with its effects becoming more pronounced at lower Darcy numbers. As the Darcy number decreases, the strength of the convective flow intensifies, leading to the distortion or bending of isotherms. Correspondingly, the velocity field is characterized by significant fluid movement, highlighting the enhanced convective activity in this regime.

Figure 4 clearly shows that increasing the angle of inclination of the cavity suppresses the convective movement. At an inclination angle of 60°, the isotherms approach linearity, indicating that conductive heat transfer predominates. Nevertheless, two opposing flows can be observed in the flow field, which illustrates the complex interplay between the angle of inclination and fluid movement in the cavity.

Figure 5 shows that increasing the Rayleigh number at low inclination angles improves the overall heat transfer. However, as the angle of inclination increases, the effectiveness of convective heat transfer is suppressed. In particular, the heat transfer is more pronounced in the case of the three-dimensional cavity than for low-dimensional configurations.

Figure 6 shows that in the case of 2D geometry at high Darcy numbers ($Da={10}^{-2}$), the total heat transfer increases with an increase in the volume fraction of the nanoparticles at inclination angles below 30°. In the 3D geometry, on the other hand, the increase in the nanoparticle volume fraction does not lead to a higher heat transfer rate for any of the Darcy numbers investigated.

5.1. Effect of the Darcy Number

Based on the average values of the Nusselt number ($Nu$) observed in both two- and three-dimensional geometries, a decrease in the Darcy number ($Da$) improves heat transfer within the cavity, which can be seen from the numerical results in Table 5 and from the temperature fields in Figure 2 and Figure 3. The Darcy number influences the value of the Darcy term in Equation (2). As the Darcy number increases, the flow regime gradually shifts towards a Darcy flow regime and thus more closely matches the properties described by Darcy’s law. At low Darcy numbers, the flow is primarily controlled by the resistance of the porous medium, resulting in a regime with minimal convective effects and significantly restricted fluid movement. This flow behavior is primarily determined by conduction rather than convection. As a result, the overall permeability of the medium is low, which promotes conductive heat transfer and suppresses advective heat transport. In this state, the system approximates the behavior of classical Darcy flow, where the effects of inertia and convection are negligible, and the flow is primarily determined by the pressure gradient and the permeability of the medium.

5.2. Effect of the Nanoparticle Volume Fraction

The solid volume fraction of nanoparticles ($\phi$) significantly influences the fluid velocity in nanofluid flows. With increasing $\phi$, the effective viscosity of the nanofluid usually increases due to the higher particle concentration, which leads to increased flow resistance. This increase in viscosity generally leads to a decrease in fluid velocity, especially in convection-dominated regions where the higher resistance dampens convective flows.

In conduction-dominated regimes, the influence of φ on velocity is less pronounced, as the flow is mainly driven by conduction rather than convection. However, in convection-dominated flows, higher φ values can lead to slower fluid movement, lower momentum transfer and weaker convection flows, which affect the overall heat transfer dynamics.

The results in Table 5 and Figure 6 show that in a 2D geometry, increasing the volume fraction of nanoparticles improves the overall heat transfer in the conduction-dominated regime ($Da>{10}^{-4}$) and at inclination angles $\alpha \le 30°.$ However, in the convection-dominated region ($Da \le {10}^{-4}$), the addition of nanoparticles reduces convective heat transfer, resulting in lower Nusselt numbers at higher nanoparticle volume fractions. In 3D geometry, increasing the nanoparticle volume fraction suppresses convective motion in both the conduction- and convection-dominated regimes.

5.3. Effect of the Inclination Angle

The inclination angle significantly influences fluid flow and heat transfer characteristics in both 2D and 3D geometries. As the inclination angle increases, the buoyancy-driven flow is altered, which can weaken or enhance convective currents depending on the specific conditions.

In the present case, the flow tends to be more robust at lower inclination angles (e.g., $\alpha \le 30°$), with strong convective motions contributing to higher heat transfer rates, especially in convection-dominated regimes. This is characterized by well-defined thermal boundary layers and increased Nusselt numbers (Figure 4). As the angle of inclination increases further (e.g., $\alpha >30°$), convective motion is increasingly suppressed, and heat transfer becomes more conduction dominated. The weakening of the convection currents at larger inclination angles leads to almost linear isotherms, as observed in Figure 4, and a lower overall heat transfer, which is reflected in lower Nusselt numbers. In extreme cases, e.g., near-vertical configurations, convection can be almost completely suppressed, resulting in a significant reduction in heat transfer efficiency.

5.4. Insights into 2D and 3D Heat Transfer Mechanisms

In the 2D configuration, an increase in the volume fraction of nanoparticles improves heat transfer in the areas dominated by thermal conduction, especially at inclination angles of α ≤ 30°. This enhancement indicates an improvement in thermal conductivity due to the incorporation of nanoparticles, which facilitates efficient heat transfer in these scenarios. Conversely, under convection-dominated conditions, the introduction of nanoparticles has a detrimental effect on convective heat transfer, resulting in a decrease in Nusselt number at higher volume fractions. This decrease can be attributed to the increased viscosity of the fluid, which hinders fluid movement and reduces the efficiency of convective heat transfer.

In contrast, the 3D geometry shows a consistent suppression of convective motion as the volume fraction of nanoparticles increases, affecting both conduction- and convection-dominated regimes. This suppression indicates a more pronounced effect of nanoparticle concentration in the 3D context, likely due to enhanced fluid interactions and increased drag caused by the presence of nanoparticles. As a result, the ability of the fluid to facilitate convective heat transfer is reduced, regardless of the predominant heat transfer mechanism.

Overall, 2D systems show improved heat transfer in conduction-dominated scenarios, with increased volume fractions of nanoparticles but a decrease in convective efficiency in convection-dominated conditions. In contrast, 3D geometries consistently suppress convective motion due to the influence of nanoparticles, highlighting the complex interplay between geometric configuration, nanoparticle volume fraction and heat transfer mechanisms, which should be further explored. This comparison emphasizes the need to consider geometric configurations when evaluating heat transfer performance in nanofluids.

6. Conclusions

A numerical analysis of natural convection in two- and three-dimensional cavities fully filled with nanofluid-saturated porous media was performed utilizing the boundary element method. A single-phase mathematical nanofluid model was used in the study, assuming a low nanoparticle concentration of up to 5% and that the nanoparticles behave analogously to water molecules. The conservation of momentum was described using the Brinkman–Forchheimer momentum equation, which takes inertial effects into account. The numerical approach integrated both the single-domain and sub-domain boundary element methods to solve the velocity–vorticity formulation of the governing equations.

The proposed numerical code was validated by comparing its results with those available in the literature over a broad spectrum of governing parameters. The study further investigated the influence of different types of nanoparticles on heat transfer enhancement, with a particular emphasis on nanoparticle volume fraction and various properties of the porous media. The results indicate that, while the addition of nanoparticles enhanced the thermal conductivity of the fluid, it generally suppressed natural convection phenomena within the porous media, particularly evident in the results obtained from the three-dimensional cavity analysis.

In the non-Darcy regime (high $Da$ values), the effect of the Brinkman viscous term in the momentum equation was significant, leading to enhanced overall heat transfer through the nanofluid-saturated porous cavity compared to the pure fluid. As the Da number decreased, the model approached the Darcy regime, and the addition of nanoparticles resulted in a reduction in overall heat transfer.

The three-dimensional analysis reveals complex interactions between the flow structure and heat transfer characteristics within the nanofluid-saturated porous cavity. The 3D temperature and velocity fields demonstrate how variations in parameters such as the Darcy number ($Da$), inclination angle (α) and nanoparticle volume fraction ($\phi$) affect convection patterns and overall thermal performance. At higher Da values in three-dimensional geometries, convective motion is more pronounced due to the increased permeability of the porous medium, which promotes stronger flow circulation and enhanced heat transfer. However, the addition of nanoparticles generally reduces the overall heat transfer in the three-dimensional cavity.

Inclination angles further modify the convection dynamics; lower angles sustain strong convective loops, whereas higher angles (e.g., α = 60°) suppress convection, leading to more stratified temperature distributions. This is characterized by nearly linear isotherms, indicating a predominance of conductive heat transfer.

The extended boundary element method presented in this study has proven to be an efficient alternative for solving complex non-linear diffusion-convection problems.

This study offers important insights into natural convection in nanofluid-saturated porous media, a field of significant industrial relevance for understanding heat transfer mechanisms. In geothermal reservoirs, nanofluids enhance thermal conductivity, improving heat recovery and efficiency. Additionally, natural convection in these materials is critical for developing thermal insulation and phase change materials for energy storage, thereby increasing energy efficiency in buildings and transportation. Furthermore, it aids in thermal management for groundwater remediation, where accurate temperature control is vital for process effectiveness.

Future challenges and directions in the field of natural convection in nanofluid-saturated porous media include the validation of theoretical and computational models through experimental studies and the evaluation of their effectiveness in practical applications. Research will focus on optimizing the properties of nanofluids, such as particle size and concentration, to improve thermal conductivity and stability. Advances in multiscale computational models will improve the prediction of heat transfer behavior. In addition, research into hybrid nanofluids and the development of innovative porous media structures with optimized pore sizes will further improve natural convection and heat transfer efficiency. These efforts aim to deepen the understanding and practical applications in this field.