Open access peer-reviewed chapter

Simulation of Natural Convection in Porous Media by Boundary Element Method

Written By

Janja Kramer Stajnko, Renata Jecl and Jure Ravnik

Submitted: 23 August 2017 Reviewed: 25 September 2017 Published: 20 December 2017

DOI: 10.5772/intechopen.71230

From the Edited Volume

Numerical Simulations in Engineering and Science

Edited by Srinivas P. Rao

Chapter metrics overview

1,204 Chapter Downloads

View Full Metrics


In this chapter, the boundary element method (BEM) is introduced for solving problems of transport phenomena in porous media domains, which is an important topic in many engineering and scientific branches as well as in fields of practical interest. The main objective of the present work is to find a numerical solution of the governing set of equations written for fluid flow in porous media domains, representing conservation of mass, momentum, and energy. The momentum equation is based on the macroscopic Navier-Stokes equations and is coupled with the energy equation. In order to use BEM for the solution of the obtained set, the governing equations are transformed by the velocity-vorticity formulation, which separates the computational scheme into kinematic and kinetic computational parts. A combination of single- and sub-domain BEM is used to solve the obtained set of partial differential equations. Solution to a problem of natural convection in porous media saturated with pure fluid and nanofluid, respectively, for examples of 2D and 3D geometries, is shown. Results are compared to published work in order to estimate the accuracy of developed numerical algorithm. Based on the results, the applicability of the BEM for solving wide range of various problems is stated.


  • boundary element method
  • porous media
  • velocity-vorticity formulation
  • natural convection
  • nanofluids

1. Introduction

Problems of transport phenomena in porous media have been widely investigated over the last few decades, mainly because of several important applications, which could be found in industry and environment, e.g., building insulation systems, dispersion of contaminants through water saturated soil, protection of groundwater resources, combustion technology. Buoyancy driven flows in porous enclosures have been simulated using different mathematical models and numerical techniques. Most commonly used mathematical model of governing momentum equation is the Darcy’s law, which is valid for the laminar flow regime (Re < 10), where the velocities are low and the viscous forces are predominant over inertia forces [1]. Extensions of the governing momentum equations have been made by analogy with the Navier-Stokes equations with addition of Brinkman term in order to consider the viscous diffusion and Forchheimer term to study the inertia effects on the free convection [2].

Problems of natural convection in porous media were studied intensively in last few decades, mainly for the cases of two-dimensional geometries. Two types of geometries are commonly investigated: porous enclosures where temperature gradient is imposed horizontally [3, 4, 5, 6, 7] or vertically [8, 9, 10, 11]. Studies considering three-dimensional geometries are rare and are usually confined on using a simplified mathematical model, e.g., Darcy model or to conditions of heating from below [12, 13, 14, 15, 16, 17, 18]. Researches considering three-dimensional cavities with the condition of heating from the side were published in [19, 20, 21].

Recently, several researchers have been investigated buoyant flow in porous media domains saturated with nanofluids [22, 23, 24]. Nanoscale particles are often added to working fluids in order to enhance heat transfer or cooling processes. A comprehensive review of the studies considering convection heat transfer in porous media saturated with nanofluid was published in [25].

The solutions of the problems of transport phenomena in porous media have been obtained using different numerical methods, where the most commonly used methods are the finite element method (FEM), the finite difference method (FDM), and the finite volume method (FVM). As an alternative to others, in engineering practice widely used methods, the BEM was developed mainly because it was very efficient for solving potential problems of fluid mechanics (inviscid fluid flow, heat conduction, etc.), where the mathematical transformation of the governing set of partial differential equation results in boundary integral equations only. To rewrite the partial differential equation into an equivalent integral representation, the known fundamental solutions of the differential operator [26] and the Green’s theorem are used. The discretized system contains only a fully populated system of integrals over boundary elements, which represent the main advantage over the volume-based methods.

When dealing with nonhomogenous and nonlinear problems, e.g., diffusion-convection problems, the domain integrals occur in the integral representation as well, which demands the extension of the classical BEM in order to additionally deal the problem within the domain. The main issue in this case is the evaluation of the domain matrices, which are full and unsymmetrical and require a lot of storage space. Several techniques have been developed in order to eliminate the domain integrals or transform them into the boundary integrals. One of the possibilities is the dual reciprocity boundary element method (DRBEM), which transforms domain integrals into a finite series of boundary integrals [27, 28, 29]. The nonhomogenous term is expanded in terms of radial basis functions. Since the discretization of the domain is represented only by grid points and the discretization of the geometry and fields on the boundary is piecewise polygonal, the DRBEM is still more flexible and efficient against other numerical methods, e.g., FDM.

Another possible extension of BEM is the boundary domain integral method (BDIM), which enables solving of strong nonlinear problems, where the domain-based effects are dominant, typical for the examples of the diffusion-convection problems [30, 31, 32]. The numerical algorithm solves the velocity-vorticity formulation of the Navier-Stokes equations, which separate the numerical scheme into the kinematic and kinetic computational parts. Consequently, the pressure is removed from the field functions conservation equations, and the calculation of the boundary pressure values is eliminated. Further advantage of BDIM is efficient dealing with boundary conditions on the solid boundaries in case of solving the vorticity equation. The vorticity is calculated explicitly from the kinematic computational part without using any approximate formulae. Using the subdomain technique [33], the problem of fully populated system matrices and corresponding memory requirements can be importantly reduced. A very stable and accurate numerical description of coupled diffusion-convection problems follows the use of Green’s functions of the appropriate linear differential operators instead of upwinding schemes of different orders, as this is the case in other domain-type numerical techniques, which also eliminate the oscillations in the numerical solutions. A wavelet compression method for a single-domain BEM in 2D was introduced in [34].

In this chapter presented numerical algorithm is based on the combination of single- and sub-domain BEM. The main advantage of the used method is that it enables an accurate prediction of the vorticity fields, which are in general defined as a curl of the velocity field. The vorticity is generated on the walls of the domain and influences the development of the flow field and furthermore the heat transfer. The single-domain BEM is used to solve the kinematics equation. The method is based on the fast multipole algorithm (FMM), which was introduced by Greengard and Rokhlin [35] for particle simulations and was later used for a wide variety of problems, e.g., for acceleration of the boundary integral Laplace equation by [36] and for coupling with BEM for the boundary matrices by [37]. The sub-domain BEM is used to solve the equations of the diffusion-advection type. A mesh of entire domain is made, where the integral equations are written for each of the sub-domains separately ([38, 39, 40, 41]). Functions are discretized using the continuous quadratic boundary elements, whereas flux is discretized using discontinuous linear boundary elements, which enable to avoid flux definition problems in corners and edges. An over-determined system of linear equations is obtained, which is solved by a least squares manner.

A numerical approach based on the BEM has been used to solve a problem of buoyancy driven flows in porous media domain, saturated with pure fluid or nanofluid. The mathematical model is based on the Navier-Stokes equations, which are averaged over the representative elementary volume and rewritten into the velocity-vorticity formulation. The influence of several governing parameters, e.g., Rayleigh number, Darcy number, and volume fraction of nanoparticles, on the heat transfer and fluid flow characteristics is analyzed.


2. Mathematical model

2.1. Governing equations

The most general mathematical model for the transport phenomena in porous media is based on the volume-averaged Navier-Stokes equations, which are primarily written on the microscopic level for the problems of pure fluid flow. Defining sufficiently large representative elementary volume (REV), with the restriction that only one part of the volume is available for the fluid flow, one can write the macroscopic set of governing equations. The REV is sufficiently large in case when it contains both, solid and fluid phases, irrespective of its position in porous media.

The presented governing equations are written for the case, when porous media are saturated with the nanofluid. The formulation enables considering both types of fluid flow by choosing correct parameter values. The properties describing nanofluids are density ρnf, dynamic viscosity μnf , heat capacitance (cp)nf, thermal expansion factor βnf, and thermal conduction knf, where index nf stands for the nanofluid. Nanofluid properties are given in relation to pure fluid and pure solid properties linked with solid volume fraction of nanoparticles ϕ, which is given as:

φ = V s V s + V f E1

where Vs is volume of solid particles and Vf the volume of fluid. Relationships between nanofluid and pure fluid properties are described with models. A comprehensive review of different models can be found in [42]. In this chapter, macroscopic modeling of nanofluids is restricted to spherical nanoparticles, and it is suitable for small temperature gradients. Density of the nanofluid ρnf is given as:

ρ nf = 1 φ ρ f + φ ρ s , E2

where index f stands for the fluid phase and index s for the solid phase. The effective dynamic viscosity μnf can be given with according to [43] as:

μ nf = μ f 1 φ 2.5 , E3

where the effective viscosity does not depend on the nanoparticle type. The heat capacitance of nanofluid can be given as:

ρ c p nf = 1 φ ρ c p f + φ ρ c p s . E4

The nanofluid thermal expansion coefficient can be written in a similar way:

ρβ nf = 1 φ ρβ f + φ ρβ s , E5

taking into account the definition of ρnf, it follows:

β nf = β f 1 1 + 1 φ ρ f φ ρ s β s β f + 1 1 + φ 1 φ ρ s ρ f . E6

The effective thermal conductivity knf is given with the Wasp model [44], which is valid only for the spherical particles, since it does not take into account the shape of the particles:

k nf = k f k s + 2 k f 2 φ k f k s k s + 2 k f + φ k f k s E7

Further assumptions for the used model are: the nanoparticles are in thermal equilibrium with the base fluid and the nonslip boundary condition is considered. The fluid flow is assumed to be laminar, steady, Newtonian, and incompressible. The density depends only on the temperature variations and can be given with the Boussinesq approximation as:

ρ nf = ρ 0 1 β nf T T 0 , E8

where T is temperature and index 0 refers to a reference state.

The problem of natural convection in saturated porous media can be described with the conservation equations for mass, momentum, and energy, written on the macroscopic level. The conservation of mass is given with the continuity equation:

v = 0 E9

where v is volume averaged velocity vector.

The momentum equation is also known as Brinkman-Forchheimer equation and reads as:

1 ϕ v t + 1 ϕ 2 v · v = 1 ρ nf p β nf T T 0 g + 1 φ µ nf ρ nf 2 v 1 K µ nf ρ nf v F v v K 1 / 2 , E10

where φ is porosity, t time, p pressure, T temperature, g gravitational acceleration, K permeability, and F Forchheimer coefficient. There are two viscous and two inertial terms in the momentum equation. The third term on the r. h. s. of the equation is the Brinkman viscous term, which is analogous to the Laplacian term in the classical Navier-Stokes equations for pure fluid flow and expresses the viscous resistance or viscous drag force exerted by the solid phase on the flowing fluid at their contact surfaces. With the Brinkman term, the nonslip boundary condition on a surface which bounds porous media is satisfied [1]. The fourth term on the r. h. s. of the momentum equation is the Darcy term, where K is the permeability, which in general depends on the geometry of the porous medium and is a second-order tensor. When assuming an isotropic porous media, the permeability is a scalar. The last term in the momentum equation is the Forchheimer inertia term, which describes the nonlinear influences at higher velocities. The Forchheimer coefficient F is a dimensionless form-drag constant and is varying with the nature of porous medium. It can be written according to Ergun model as ([1]):

K = ϕ 3 d p 2 a 1 ϕ 2 , F = b a ϕ 3 ,E11

where a and b are Ergun’s constants with values a = 150 and b = 1.75, while dp is the average particle size of the bed.

Finally, the energy equation can be written as:

σ T t + v · T = k e ρ c p nf 2 T , E12

where σ is the specific heat ratio, σ = ϕ + (1 − ϕ)(ρcp)p/(ρcp)nf , (ρcp)p and (ρcp)nf are heat capacitances of solid and fluid phase, respectively. Furthermore, ke is the effective conductivity of porous medium. It is assumed that the thermal properties of solid matrix and the nanofluid are identical [24, 45], resulting in σ = 1 and ke = knf .

Governing equations (9), (10), and (12) can be converted into a nondimensional form by introduction of the following dimensionless variables:

v v v 0 , r r L , t v 0 t L , g g g 0 , p p p 0 , T T T 0 T E13

The parameters in the above expressions are v0 characteristic velocity given as v0 = kf /(ρcp)f L, which is common choice for buoyant flow simulations, kf is the fluid thermal conductivity, (ρcp)f is the heat capacity for the fluid phase, and L is the characteristic length. Moreover, T0 is characteristic temperature T0 = (T2 − T1)/2 and ΔT is characteristic temperature difference ΔT = T2 − T1, p0 is the characteristic pressure p0 = 1bar, while gravitational acceleration is g0 = 9.81 m/s2.

In addition, the velocity-vorticity formulation of the governing equations is proposed by introduction of the vorticity vector, which is by the definition a curl of the velocity field ω = × v . The governing set of equations in nondimensional velocity-vorticity formulation can now be written in terms of kinematics equation, the vorticity transport equation, and the energy equation as:

2 v + × ω = 0 E14
v ω = ω v C A Pr Ra T ϕ 2 × T g + C B Pr ϕ 2 ω C B Pr Da ϕ 2 ω F Da ϕ 2 v ω , E15
v · T = C C 2 T E16

In the above equations, parameters CA, CB, and CC are presenting the nanofluid properties and are given with expressions:

C A = μ nf μ f ρ f ρ nf , C B = β nf β f , C C = α nf α f ,E17

where αnf is thermal diffusivity of nanofluid, αnf = knf /(ρcp)nf and αf thermal diffusivity of pure fluid αf = kf /(ρcp)f . The nanofluid properties are obtained using the expressions (2)–(8). For the simulation of the pure fluid flow, the parameters are CA = CB = CC = 1. The time derivatives in the vorticity and energy equations ω / t , ∂T/∂t are omitted, since only steady flow simulations are shown in the present chapter.

The nondimensional parameters appearing in the momentum equation are:

  • Fluid Rayleigh number RaT = g βT ΔT L3 ρf (ρcp)f /μf kf ;

  • Prandtl number Pr = μf cp/kf ;

  • Darcy number Da = K/L2;

The results are presented in terms of the porous Rayleigh number RaP, which links the thermal Rayleigh number and Darcy Number:

  • RaP = RaT ∙ Da.

2.2. Boundary element method

The governing set of equations (14) and (15) in (16) is solved using an algorithm based on the combination of single-domain and sub-domain BEM, primarily developed for pure fluid flow simulations [39, 40] and later adopted for nanofluids [46] and porous media flow simulations [21]. The algorithm solves the velocity-vorticity formulation of Navier-Stokes equations. The sub-domain BEM solves the vorticity and energy transport equations. It is based on the domain decomposition, which results in sparse matrices and improves the efficiency of the solution to become comparable to FVM or FEM [47]. The kinematics equation for the calculation of the boundary vorticities is solved by a single-domain BEM. This results in full system of equations and limits the maximum grid size due to memory constraints. This drawback can be mitigated using the fast BEM, where sparse approximation of full matrices is used [37]. The main advantage of using the single-domain BEM for the boundary vorticity values is that the algorithm conserves mass in complex geometries, which is not the case when using velocity derivatives to calculate boundary vorticity values.

The numerical algorithm is devised as follows. At the beginning, the boundary conditions for the velocity and temperature are required and have to be given in terms of Dirichlet and Neumann type. In addition, the temperature and temperature flux on the solid walls and the no-slip boundary conditions are prescribed. The boundary conditions for the vorticity are unknown at the beginning and are calculated later as a part of numerical algorithm. The known boundary conditions are used to solve the kinematics equation (14) for the domain velocity values and energy equation (16) for the domain temperature values. The boundary vorticity values are first obtained using the single-domain BEM on the kinematics equation; moreover, the domain vorticity values are obtained out of vorticity transport equation (15) using a sub-domain BEM.

The outline of the numerical algorithm is as follows:

  1. Step 1. Fluid/nanofluid and porous media properties are determined.

  2. Step 2. Vorticity values on the boundary are calculated by a fast single-domain BEM from the kinematics equation (14).

  3. Step 3. Velocity values within the domain are calculated by a sub-domain BEM from the kinematics equation (14).

  4. Step 4. Temperature values are calculated by a subdomain BEM from the energy equation (16).

  5. Step 5. Vorticity values within the domain are calculated by a subdomain BEM from the vorticity equation (17).

  6. Step 6. Convergence check; all steps from 2 until 5 are repeated until all flow fields achieve the required accuracy.

In order to apply the proposed algorithm, governing equations have to be written in integral form. The integral representation is obtained using the Green’s second identity for the unknown field function and for the fundamental solution of the Laplace equation as proposed in [48].

2.2.1. Integral representation of the kinematics equation

For the unknown boundary vorticity values, the single-domain BEM is used on the kinematics equation (14). The integral representation of the kinematics equation in its tangential form is:

c ξ n ξ × v ξ + n ξ × Γ v u n = n ξ × Γ v × n × u d Γ + n ξ × Ω ω × u d Ω E18

where Ω is the computational domain and Γ = ∂Ω is the boundary of the domain, c ξ is geometric factor defined as c ξ = θ / 4 π , θ is the inner angle with origin in ξ . If ξ lies inside the domain, then c ξ =1, and if ξ lies on a smooth boundary, then c ξ =1/2. Furthermore, n is a vector normal to the boundary, and u is the fundamental solution of the Laplace equation given as:

u = 1 4 π ξ r E19

The discretized system of equations is written for unknown boundary vorticities, while domain vorticity and velocity values are taken from the previous nonlinear iteration. The source point is set into every boundary node of the whole computational domain, which follows in full system matrix, where number of rows and columns are equal to number of boundary nodes. The system is solved using a LU decomposition method. The storage requirements are reduced with a kernel expansion approximation technique [49].

In addition, the kinematics equation (14) is used again in order to calculate domain velocity values with the sub-domain BEM. The following form of the integral equation is used:

c ξ v ξ + Γ v n u = Γ v × n × u + Ω ω × u . E20

The obtained integral kinematics equation is without the derivatives of the velocity or vorticity fields, which enables that the source point is set to function the nodes only. The domain velocity values are calculated based on the known boundary values of the velocity from the initial boundary conditions, while the domain and boundary values of the vorticity are known from the previous iteration.

2.2.2. Integral representation of the vorticity and energy equations

In order to derive the integral form of the vorticity and energy equations, the same fundamental solution of the Laplace equation is used [26]. The final integral form of the vorticity transport equation is:

c ( ξ ) ω j ( ξ ) + Γ ω j u * n   d Γ = Γ u * q j   d Γ + 1 P r 1 C B 1 φ Γ n { u * ( v ω j ω v j ) }   d Γ 1 P r 1 C B 1 φ Ω ( v ω j ω v j ) u *   d Ω     R a T C A C B   φ Γ ( u * T g × n ) j   d Γ R a T C A C B   φ Γ ( T × u * g ) j   d Ω     + 1 D a φ Ω ω j u * d Ω + F Pr D a 1 C B φ   | v | Ω ω j u * d Ω , E21

and finally the integral form of the energy transport equation reads as:

c ξ T ξ + Γ T u n = Γ u q T + 1 C C Γ n u v T Ω v T u E22

In the above equations, qj is a component of vorticity flux, whereas qT is a heat flux. In the sub-domain BEM method, a mesh of the entire domain Ω is made, each mesh element is named a subdomain. All equations are written for each of the subdomains. The filed functions and flux across the boundary and within the domain are interpolated using shape functions. The hexahedral subdomains with 27 nodes are used, enabling continuous quadratic interpolation of field functions. The field functions on each element are interpolated using continuous quadratic interpolation, while fluxes are interpolated using the discontinuous linear interpolation. With discontinuous interpolation, the definition problems in corners and edges are avoided.


3. Test cases

The physical models where the above developed numerical scheme was tested are a two-dimensional rectangular enclosure and a three-dimensional cubical enclosure filled with fully saturated porous medium. Porous medium is assumed to be isotropic, homogenous, and in thermal equilibrium with the fluid phase. The simulation of fluid flow and heat transfer through porous media domain of pure fluid and nanofluid, respectively, is presented. Two opposite vertical walls are subjected to a temperature differences, while the rest of the walls is adiabatic and impermeable. Geometry with corresponding boundary conditions is shown in Figure 1. The boundary conditions for the current problem are:

v x = v y = 0 , ω = 0 , T = T H at x = 0 , v x = v y = 0 , ω = 0 , T = T C at x = L , v x = v y = 0 , ω = 0 , T y = 0 at y = 0 and z = 0 , v x = v y = 0 , ω = 0 , T y = 0 at y = L and z = L E23

Figure 1.

Two-dimensional and three-dimensional enclosures with corresponding boundary conditions.

Due to applied temperature gradient, density differences are induced, which result in appearance of thermal buoyancy force producing a large vortex in the main part of the cavity.

The overall heat transfer through porous media is expected to depend on several fluid and porous media properties, such as porosity, permeability, thermal conductivity, heat capacitance, solid volume fraction of nanoparticles, and types of nanoparticles. In order to compare different conditions on the heat transfer characteristics, the wall heat flux is calculated, which is expressed in terms of the dimensionless Nusselt number as:

Nu = k nf k f Γ T n , E24

where Γ is the surface through which the heat flux is calculated and n is the unit normal to this surface. The definition is valid for nanofluids as well as for pure fluids, since there the ratio of thermal conductivities is knf /kf = 1.

In the present study, the Cu nanoparticles are added to the water as a base fluid. The thermophysical properties of Cu nanoparticles and water are given in Table 1 [50].

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

Table 1.

Thermophysical properties of water and Cu nanoparticles.

In order to obtain a grid-independent solution, at the beginning, a grid sensitivity analysis was performed. Two nonuniform meshes for 2D geometry and four nonuniform meshes for 3D geometry have been tested. The results are shown for the case, when porous media are saturated with pure fluid and parameters RaP = 100, Pr  = 0.71, φ = 0.8, σ = 1 and various values of Da. The results are presented in Table 2.

Mesh Number of nodes RaP = 100, Pr = 0.71, φ = 0.8
Da 10−1 10−2 10−3 10−4 10−5
2D 20 × 20 1681 1.0639 1.6329 2.3697 2.8756 3.1656
30 × 30 3721 1.0638 1.6331 2.3680 2.8537 3.1503
3D 12 × 12 × 12 15,625 1.0423 1.5428 2.3432 2.9784 3.3008
20 × 8 × 20 28,577 1.0394 1.5329 2.3313 2.9575 3.2950
22 × 10 × 22 42,525 1.0393 1.5327 2.3307 2.9552 3.2945
30 × 10 × 30 78,141 1.0393 1.5325 2.3303 2.9541 3.2934

Table 2.

Variations of Nusselt number with different grid sizes and various Darcy numbers.

When comparing 2D and 3D results, it can be observed that for the case of low values of Darcy number, 2D simulation underestimates the heat transfer up to 4.5%. Based on the presented results, the 20 × 20 mesh for 2D simulations and 20 × 8 × 20 for 3D simulations were chosen.

The validation of numerical code has been primarily performed for the pure fluid saturating the porous media. The results for 2D geometry, compared to [24, 51], are presented in Table 3. Results for 3D geometry are compared to [20] and are presented in Table 4. Furthermore, in Table 5, the results for natural convection in 2D enclosure for a nanofluid saturated porous media are shown and compared to [24]. According to the comparable study, the Cu nanoparticles were added to the water as a base fluid.

Da RaP [51] [24] Present
10−2 10 1.015 1.010 1.012
100 1.530 1.533 1.503
1000 3.555 3.602 3.499
10−4 10 1.071 1.065 1.070
100 2.725 2.764 2.777
1000 9.183 9.454 9.174
10−6 10 1.079 1.072 1.093
100 2.997 2.980 3.241
1000 11.790 11.924 12.895

Table 3.

Validation of the numerical code by a comparison of average Nu for natural convection in porous media saturated with a pure fluid for Pr = 1.0, φ=0.6 and different Da and RaP, for 2D geometry.

RaP = 1000, Pr = 0.71, φ = 0.8
Da 10−2 10−3 10−4 10−5 10−6
[20] 3.99 6.95 10.14 12.78 13.72
Present 3.770 6.922 10.558 13.242 14.568

Table 4.

Nusselt number values for the 3D natural convection in a cubic enclosure filled with porous media saturated with pure fluid.

Da RaP [24] Present [24] Present [24] Present
ϕ = 0.05
φ = 0.4 φ = 0.6 φ = 0.9
10−2 1000 3.433 3.400 3.850 3.826 4.162 4.145
10−4 1000 9.117 9.132 9.590 9.743 9.901 10.154
10−6 1000 11.778 12.991 11.899 13.128 11.976 13.195
φ = 0.4
ϕ = 0.0 ϕ = 0.025 ϕ = 0.05
10−2 10 1.007 1.008 1.081 1.083 1.160 1.162
10−2 1000 3.302 3.282 3.370 3.345 3.433 3.400
10−6 1000 11.867 13.238 11.847 13.131 11.778 12.991

Table 5.

Nusselt number values for a natural convection in porous media saturated with nanofluid in 2D enclosure for various governing parameters (Pr = 6.2).

From the comparison, it can be observed that the results agree well with the data from the published studies, which confirm accuracy of the obtained numerical algorithm.

The isotherms for Cu-water nanofluid under different values of porous Rayleigh number and Darcy number at porosity φ = 0.4 and different values of solid volume fraction are shown in Figure 2. Solid lines correspond to ϕ = 0.0, dotted lines to ϕ = 0.025, and dashed lines to ϕ = 0.05. Heat transfer in porous medium is mostly affected by a Rayleigh and Darcy numbers. At RaP = 10, the heat transfer in horizontal direction is weak, and the main heat transfer mechanism in this case is conduction. Increase of RaP results in stronger convective motion, which is clearly evident from the temperature field; when RaP = 1000, thin boundary layers are created near the hot and cold walls, and the isotherms in the core region become almost horizontal and parallel to adiabatic and impermeable walls. According to the temperature fields, decrease of Da enhances the heat transfer through cavity. The Da number influences the magnitude of the Darcy term in the vorticity equation (10). With increase of Da, the flow regime is transited into the Darcy flow regime, which can be described by the Dracy’s law.

Figure 2.

Temperature fields for different values of RaP for Cu nanofluid and different solid volume fractions for Pr = 6.2, Da = 10−6; solid lines are for ϕ = 0.0, dashed lines for ϕ = 0.025 and dotted lines for ϕ = 0.05.

According to the results, the addition of nanoparticles into the water causes the attenuation of the convective motion. However, the overall heat transfer is enhanced with increase of solid volume fraction of nanoparticles in cases of conduction-dominated regimes (low values of Rayleigh numbers RaP < 100 and high values of Darcy numbers Da > 10−4). On the other hand in convection dominated regimes (RaP > 100, Da < 10−4), the addition of nanoparticles diminishes the convection, which results in lower values of Nusselt numbers. Figure 3 shows the dependence of Nu on porous Rayleigh number and solid volume fraction of nanoparticles. It can be observed that for any values of RaP with increase of ϕ, the heat transfer increases.

Figure 3.

Nusselt number values depending on porous Rayleigh number for Pr = 6.2, Da = 10−6 and different values of solid volume fraction of nanoparticles.

When observing the flow structure in 3D enclosure, it is obvious that the flow field is not far from being 2D, which is a consequence of the applied temperature difference between the opposite walls, which causes large two-dimensional vortex in the y plane. In order to observe the 3D nature of the phenomena, the iso-surfaces of the absolute value of y velocity component are shown in Figure 4. The extent of movement perpendicular to the plane of the main vortex is small, but it becomes more apparent in case of higher values of RaP and lower values of Da.

Figure 4.

Iso-surfaces for RaP = 500, Da = 10−3 and absolute value of velocity component |vy| = 3 (left) and RaP = 1000, Da = 10−3, |vy| = 7 (right). Contours of temperature are displayed on the iso-surfaces (−0.5 < T < 0.5). In addition, the velocity vectors on the plane y = 0.5 are displayed.


4. Conclusion

The numerical method based on the BEM is presented for solving the coupled set of partial differential equations, which describe the fluid flow and heat transfer in porous medium domain. The mathematical model is based on the Navier-Stokes equations, which are averaged over the representative elementary volume. The proposed numerical algorithm solves the velocity-vorticity formulation of the governing equations. The numerical scheme is split into the single-domain BEM, which solves the kinematic equation for unknown boundary vorticity values and sub-domain BEM for domain velocity, vorticity, and temperature values.

The numerical algorithm is tested on an example of natural convection phenomena in porous media domain for a 2D as well as 3D geometry. Porous media are fully saturated with pure fluid or water-based nanofluid with addition of Cu nanoparticles. Obtained numerical results were validated with available benchmark solutions.

The natural convection phenomena strongly depend on the parameters, e.g., Rayleigh number and Darcy number. Addition of nanoparticles to a base fluid enhances the heat transfer through porous media, when conduction is the dominant heat transfer mechanism. On the other hand, in convection dominated regime, the addition of nanoparticles reduces the magnitude of convective motion.

The good agreement of the results with the published ones confirms the efficiency of the BEM-based methods as a powerful alternative to the existing numerical methods.


  1. 1. Nield DA, Bejan A. Convection in Porous Media. 4th ed. United States of America: Springer; 2013
  2. 2. Lauriat G, Prasad V. Natural convection in a vertical porous cavity: A numerical study for Brinkman-extended Darcy formulation. Journal of Heat Transfer. 1987;109:688-696
  3. 3. Beckermann C, Viskanta R, Ramadhyani S. A numerical study of non-Darcian natural convection in a vertical enclosure filled with a porous medium. Numerical Heat Transfer. 1986;10:557-570
  4. 4. Lauriat G, Prasad V. Non-Darcian effects on natural convection in a vertical porous enclosure. International Journal of Heat and Mass Transfer. 1989;32:2135-2148
  5. 5. Jecl R, Škerget L, Petrešin E. Boundary domain integral method for transport phenomena in porous media. International Journal for Numerical Methods in Fluids. 2001;35:39-54
  6. 6. Baytas AC, Pop I. Free convection in a square porous cavity using a thermal nonequilibrium model. International Journal of Thermal Sciences. 2002;41:861-870
  7. 7. Basak T, Roy S, Paul T, Pop I. Natural convection in a square cavity filled with a porous medium: Effects of various thermal boundary conditions. International Journal of Heat and Mass Transfer. 2006;49:1430-1441
  8. 8. Prasad V, Kulacki FA. Natural convection in horizontal porous layers with localized heating from below. Journal of Heat Transfer. 1987;109:795-798
  9. 9. Kladias N, Prasad V. Natural convection in horizontal porous layers: Effects of Darcy and Prandtl numbers. Journal of Heat Transfer. 1989;111:926-935
  10. 10. Rosenberg ND, Spera FJ. Thermohaline convection in a porous medium heated from below. International Journal of Heat and Mass Transfer. 1992;35:1261-1273
  11. 11. Storesletten L. Natural convection in a horizontal porous layer with anisotropic thermal diffusivity. Transport in Porous Media. 1993;12:19-29
  12. 12. Beck JL. Convection in a box of porous material saturated with fluid. Physics of Fluids. 1972;15:1377-1383
  13. 13. Holst PH, Aziz K. Transient three-dimensional natural convection in confined porous media. International Journal of Heat and Mass Transfer. 1972;15:73-90
  14. 14. Zebib a, Kassoy DR. Onset of natural convection in a box of water-saturated porous media with large temperature variation. The Physics of Fluids. 1977;20:4-9
  15. 15. Zebib A, Kassoy DR. Three-dimensional natural convection motion in a confined porous medium. Physics of Fluids. 1978;21:1-3
  16. 16. Horne RN. Three-dimensional natural convection in a confined porous medium heated from below. Journal of Fluid Dynamics. 1979;92:25-38
  17. 17. Neto HL, Quaresma JNN. Natural convection in three-dimensional porous cavities: Integral transform method. International Journal of Heat and Mass Transfer. 2002;45:3013-3032
  18. 18. Sezai I. Flow patterns in a fluid-saturated porous cube heated from below. Journal of Fluid Mechanics. 2005;523:393-410
  19. 19. Dawood AS, Burns PJ. Steady three-dimensiional convective heat transfer in a porous box via multigrid. Numerical Heat Transfer, Part A. 1992;22:167-198
  20. 20. Sharma RV, Sharma RP. Non-Darcy effects on three-dimensional natural convection in a porous box. Annals of the Assembly for International Heat Transfer Conference. 2006
  21. 21. Kramer J, Ravnik J, Jecl R, Škerget L. Simulation of 3D flow in porous media by boundary element method. Engineering Analysis with Boundary Elements. 2011;35:1256-1264
  22. 22. Shermet MA, Pop I. Conjugate natural convection in a sqare porous cavity filled by a nanofluid using Buongiorno's mathematical model. International Journal of Heat and Mass Transfer. 2014;79:137-145
  23. 23. Grosan T, Revnic C, Pop I, Ingham DB. Free convection heat transfer in a square cavity filled with a porous medium saturated by a nanofluid. International Journal of Heat and Mass Transfer. 2015;87:36-41
  24. 24. Nguyen MT, Aly AM, Lee SW. Natural convection in a non-Darcy porous cavity filled with Cu-water nanofluid using the characteristic-based split procedure in finite-element method. Numerical Heat Transfer, Part A. 2015;67:224-247
  25. 25. Mahdi RA, Mohammed HA, Munisamy KM, Saeid NH. Review of convection heat transfer and fluid flow in porous media with nanofluid. Renewable and Sustainable Energy Reviews. 2015;41:715-734
  26. 26. Wrobel LC. The Boundary Element Method. Chichester, England, New York: John Wiley & Sons Ltd; 2002
  27. 27. Partridge P, Brebbia C, Wrobel L. The Dual Reciprocity Boundary Element Method. Southampton: Computational Mechanics Publications; 1992
  28. 28. Perez-Gavilan JJ, Aliabadi MH. A Galerkin boundary element formulation with dual reciprocity for elastodynamics. International Journal of Numerical Methods in Engineering. 2000;48:1331-1344
  29. 29. Blobner J, Hriberšek M, Kuhn G. Dual reciprocity BEM-BDIM technique for conjugate heat transfer computations. Computer Methods in Applied Mechanics and Engineering. 2000;190:1105-1116
  30. 30. Škerget L, Hriberšek M, Kuhn G. Computational fluid dynamics by boundary domain integral method. International Journal for Numerical Methods in Engineering. 1999;46:1291-1311
  31. 31. Jecl R, Škerget L, Petrešin E. Boundary domain integral method for transport phenomena in porous media. International Journal for Numerical Methods in Engineering. 2001;35:39-54
  32. 32. Kramer J, Jecl R, Škerget L. Boundary domain integral method for study of double diffusive natural convection in porous media. Engineering Analysis with Boundary Elements. 2007;31:897-905
  33. 33. Hriberšek M, Škerget L. Iterative methods in solving Navier-stokes equations by the boundary element method. International Journal for Numerical Methods in Engineering. 1996;39:115-139
  34. 34. Ravnik J, Škerget L, Hriberšek M. The wavelet transform for BEM computational fluid dynamics. Engineering Analysis with Boundary Elements. 2004;28:1303-1314
  35. 35. Greengard L, Rokhlin V. A fast algorithm for particle simulations. Journal of Computational Physics. 1997;135:280-292
  36. 36. Popov V, Power H, Walker SP. Numerical comparison between two possible multipole alternatives for the BEM solution of 3D elasticity problems based upon Taylor series expansions. Engineering Analysis with Boundary Elements. 2003;27:521-531
  37. 37. Ravnik J, Škerget L, Žunič Z. Comparison between wavelet and fast multipole data sparse approximations for Poisson and kinematic boundary-domain integral equations. Computer Methods in Applied Mechnaics and Engineering. 2009;198:1473-1485
  38. 38. Ramšak M, Škerget L. 3D multidomain BEM for solving the Laplace equation. Engineering Analysis with Boundary Elelments. 2007;31:528-538
  39. 39. Ravnik J, Škerget L, Žunič Z. Velocity-vorticity formulation for 3D natural convection in an inclined enclosure by BEM. International Journal of Heat and Mass Transfer. 2008;51:4517-4527
  40. 40. Ravnik J, Žunič Z. Combined single domain and sudomain BEM for 3D laminar viscous flow. Engineering Analysis with Boundary Elements. 2009;33:420-424
  41. 41. Popov V, Power H, Škerget L, editors. Domain Decomposition Techniques for Boundary Elements: Applications to Fluid Flow. Southampton, Billerica: WIT Press; 2007
  42. 42. Haddad Z, Oztop HF, Abu-Nada E, Mataoui A. A review on natural convective heat transfer of nanofluids. Renewable and Sustainable Energy Reviews. 2012;16:5363-5378
  43. 43. Brinkman HC. The viscosity of concentrated suspensions and solutions. The Journal of Chemical Physics. 1952;20:571-581
  44. 44. Wasp FJ. Solid-liquid slurry pipeline transportation. Clausthal Trans Tech Publications; 1977
  45. 45. Bourantas GC, Skouras ED, Loukopoulos VC, Burganos VN. Heat transfer and natural convection of nanofluids in porous media. European Journal of Mechanics B/Fluids. 2014;43:45-56
  46. 46. Ravnik J, Hriberšek M. Analysis of three-dimensional natural convection of nanofluids by BEM. Engineering Analysis with Boundary Elements. 2010;34:1018-1030
  47. 47. Ravnik J, Hriberšek M. 2D velocity vorticity based LES for the solution of natural convection in a differentially heated enclosure by wavelet transform based BEM and FEM. Engineering Analysis with Boundary Elements. 2006;30:671-686
  48. 48. Škerget L, Žunič Z. Natural convection flows in complex cavities by BEM. International Journal of Numerical Methods for Heat & Fluid Flow. 2003;13:720-735
  49. 49. Ravnik J, Žunič Z. Fast single domain-subdomain BEM algorithm for 3D incompressible fluid flow and heat transfer. International Journal for Numerical Methods in Engineering. 2009;77:1627-1645
  50. 50. Oztop HF, Abu-Nada E. Natural convection of water based nanofluids in an inclined enclosure with heat source. International Journal of Heat and Fluid Flow. 2008;29:1326-1336
  51. 51. Nithiarasu P, Seetharamu KN, Sundararajan T. Natural convective heat transfer in a fluid saturated variable porosity medium. International Journal of Heat and Mass Transfer. 1997;40:3955-3967

Written By

Janja Kramer Stajnko, Renata Jecl and Jure Ravnik

Submitted: 23 August 2017 Reviewed: 25 September 2017 Published: 20 December 2017