Open access peer-reviewed chapter

# Simulation of Natural Convection in Porous Media by Boundary Element Method

By Janja Kramer Stajnko, Renata Jecl and Jure Ravnik

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

DOI: 10.5772/intechopen.71230

## Abstract

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.

### Keywords

• 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:

φ=VsVs+VfE1

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=μf1φ2.5,E3

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

ρcpnf=1φρcpf+φρcps.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=βf11+1φρfφρsβsβf+11+φ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:

knf=kfks+2kf2φkfksks+2kf+φkfksE7

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=ρ01βnfTT0,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=0E9

where vis volume averaged velocity vector.

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

1ϕvt+1ϕ2v·v=1ρnfpβnfTT0g+1φµnfρnf2v1KµnfρnfvFvvK1/2,E10

where φ is porosity, t time, p pressure, T temperature, ggravitational 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=ϕ3dp2a1ϕ2,F=baϕ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:

σTt+v·T=keρcpnf2T,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:

vvv0,rrL,tv0tL,ggg0,ppp0,TTT0TE13

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:

2v+×ω=0E14
vω=ωvCAPrRaTϕ2×Tg+CBPrϕ2ωCBPrDaϕ2ωFDaϕ2vω,E15
v·T=CC2TE16

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

CA=μnfμfρfρnf,CB=βnfβf,CC=α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ξ×Γvun=nξ×Γv×n×udΓ+nξ×Ωω×udΩ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, nis a vector normal to the boundary, and u is the fundamental solution of the Laplace equation given as:

u=14πξrE19

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ξ+Γvnu=Γ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(ξ)+Γωju*n dΓ=Γu*qj dΓ+1Pr1CB1φΓn{u*(vωjωvj)} dΓ1Pr1CB1φΩ(vωjωvj)u* dΩ   RaTCACB φΓ(u*Tg×n)j dΓRaTCACB φΓ(T×u*g)j dΩ   +1DaφΩωju*dΩ+FPrDa1CBφ |v|Ωωju*dΩ,E21

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

cξTξ+ΓTun=ΓuqT+1CCΓnuvTΩvTuE22

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:

vx=vy=0,ω=0,T=THatx=0,vx=vy=0,ω=0,T=TCatx=L,vx=vy=0,ω=0,Ty=0aty=0andz=0,vx=vy=0,ω=0,Ty=0aty=Landz=LE23

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=knfkfΓTn,E24

where Γ is the surface through which the heat flux is calculated and nis 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]
Water4179997.10.613211.47
Cu38589334001.671163

### 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.

MeshNumber of nodesRaP = 100, Pr = 0.71, φ = 0.8
Da10−110−210−310−410−5
2D20 × 2016811.06391.63292.36972.87563.1656
30 × 3037211.06381.63312.36802.85373.1503
3D12 × 12 × 1215,6251.04231.54282.34322.97843.3008
20 × 8 × 2028,5771.03941.53292.33132.95753.2950
22 × 10 × 2242,5251.03931.53272.33072.95523.2945
30 × 10 × 3078,1411.03931.53252.33032.95413.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.

DaRaP[51][24]Present
10−2101.0151.0101.012
1001.5301.5331.503
10003.5553.6023.499
10−4101.0711.0651.070
1002.7252.7642.777
10009.1839.4549.174
10−6101.0791.0721.093
1002.9972.9803.241
100011.79011.92412.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
Da10−210−310−410−510−6
[20]3.996.9510.1412.7813.72
Present3.7706.92210.55813.24214.568

### Table 4.

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

DaRaP[24]Present[24]Present[24]Present
ϕ = 0.05
φ = 0.4φ = 0.6φ = 0.9
10−210003.4333.4003.8503.8264.1624.145
10−410009.1179.1329.5909.7439.90110.154
10−6100011.77812.99111.89913.12811.97613.195
φ = 0.4
ϕ = 0.0ϕ = 0.025ϕ = 0.05
10−2101.0071.0081.0811.0831.1601.162
10−210003.3023.2823.3703.3453.4333.400
10−6100011.86713.23811.84713.13111.77812.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.

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.

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.

## 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.

## More

© 2017 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## How to cite and reference

### Cite this chapter Copy to clipboard

Janja Kramer Stajnko, Renata Jecl and Jure Ravnik (December 20th 2017). Simulation of Natural Convection in Porous Media by Boundary Element Method, Numerical Simulations in Engineering and Science, Srinivas P. Rao, IntechOpen, DOI: 10.5772/intechopen.71230. Available from:

### Related Content

#### Numerical Simulations in Engineering and Science

Edited by Srinivasa Rao

Next chapter

#### Numerical Simulations of a High-Resolution RANS-FVDM Scheme for the Design of a Gas Turbine Centrifugal Compressor

By Chellapilla V K N S N Moorthy and Vadapalli Srinivas

#### Applications of Monte Carlo Method in Science and Engineering

Edited by Shaul Mordechai

First chapter

#### Monte Carlo Simulations in NDT

By Frank Sukowski and Norman Uhlmann

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

View all books