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

$$
\varphi = \frac{V\_\*}{V\_\* + V\_f} \tag{1}
$$

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:

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 tem-

*ρnf* = *ρ*0(1 − *βnf* (*T* − *T*0)), (8)

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

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

⃑ *<sup>p</sup>* <sup>−</sup> *<sup>β</sup>nf*(*<sup>T</sup>* <sup>−</sup> *<sup>T</sup>*0) *<sup>g</sup>*⃑ <sup>+</sup> \_\_<sup>1</sup>

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

> 2 \_\_\_\_\_\_

*<sup>a</sup>* (<sup>1</sup> <sup>−</sup> *<sup>ϕ</sup>*)2, *<sup>F</sup>* <sup>=</sup> \_\_\_*<sup>b</sup>*

√ \_\_\_\_

*<sup>a</sup> <sup>ϕ</sup>*<sup>3</sup> , (11)

*φ μnf* \_\_\_ *ρnf* <sup>∇</sup><sup>2</sup> *<sup>v</sup>*⃑ <sup>−</sup> \_\_<sup>1</sup> *K μnf* \_\_\_ *ρnf <sup>v</sup>*⃑ <sup>−</sup> *<sup>F</sup> <sup>v</sup>*⃑ |*v* \_\_\_\_\_⃑ |

perature variations and can be given with the Boussinesq approximation as:

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

conservation of mass is given with the continuity equation:

⃑ ) *<sup>v</sup>*⃑ <sup>=</sup> \_\_\_<sup>1</sup> *ρnf* ∇

∇

where *v*⃑ is volume averaged velocity vector.

\_\_<sup>1</sup>

Ergun model as ([1]):

*<sup>K</sup>* <sup>=</sup> *<sup>ϕ</sup>*<sup>3</sup> *dp*

*ϕ* ∂*v*⃑ \_\_\_ <sup>∂</sup>*<sup>t</sup>* <sup>+</sup> \_\_<sup>1</sup> *ϕ*2 (*v*⃑ ∙ ∇ *ks* + 2 *kf* − 2*φ*(*kf* − *ks*)

\_\_\_\_\_\_\_\_\_\_\_\_\_ *ks* <sup>+</sup> <sup>2</sup> *kf* <sup>+</sup> *<sup>φ</sup>*(*kf* <sup>−</sup> *ks*) (7)

http://dx.doi.org/10.5772/intechopen.71230

83

Simulation of Natural Convection in Porous Media by Boundary Element Method

⃑ ∙ *v*⃑ = 0 (9)

*<sup>K</sup>*<sup>1</sup>/2 , (10)

*knf* = *kf*

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:

$$
\rho\_{\text{uf}} = \begin{pmatrix} 1 - \varphi \end{pmatrix} \rho\_f + \varphi \lrcorner \rho\_{\varphi} \tag{2}
$$

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:

$$
\mu\_{\eta^f} = \frac{\mu\_f}{(1-\phi)^{2\mathcal{N}}} \tag{3}
$$

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

$$\left(\rho c\_{p}\right)\_{n\!\!/} = \left(1-\phi\right)\left(\rho c\_{p}\right)\_{\!\!/} + \phi\left(\rho c\_{p}\right)\_{\!\!/} \tag{4}$$

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

$$(\rho \emptyset)\_{\eta\zeta} = \left(1 - \wp\right)(\rho \emptyset)\_{\zeta} + \wp \left(\rho \emptyset\right)\_{\zeta} \tag{5}$$

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

$$\beta\_{nf} = \beta\_f \left[ \frac{1}{1 + \frac{(1 - \rho)\rho\_f}{\rho\rho\_s} \frac{\beta\_s}{\beta\_f}} + \frac{1}{1 + \frac{\rho}{1 - \rho} \frac{\rho\_s}{\beta\_f}} \right]. \tag{6}$$

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 \frac{k\_\ast + 2\,k\_\circ - 2\,\rho(k\_\circ - k\_\circ)}{k\_\ast + 2\,k\_\circ + \rho(k\_\circ - k\_\circ)}\tag{7}$$

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:

$$\rho\_{\ast \circ'} = \rho\_o (1 - \beta\_{\ast \circ'} (T - T\_o))\_{\prime} \tag{8}$$

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:

$$
\dot{\nabla} \cdot \ddot{\boldsymbol{v}} = \mathbf{0} \tag{9}
$$

where *v*⃑ is volume averaged velocity vector.

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

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 viscos-

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:

*Vs* + *Vf*

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.

where index *f* stands for the fluid phase and index *s* for the solid phase. The effective dynamic

where the effective viscosity does not depend on the nanoparticle type. The heat capacitance

<sup>=</sup> (<sup>1</sup> <sup>−</sup> *<sup>φ</sup>*) (*ρcp*)*<sup>f</sup>*

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

[

1 +

\_\_\_\_\_\_\_ 1

(<sup>1</sup> <sup>−</sup> *<sup>φ</sup>*) *<sup>ρ</sup>* \_\_\_\_\_\_*<sup>f</sup> φρ<sup>s</sup>*

*β* \_\_*s βf*

()*nf* = (1 − *φ*) ()*<sup>f</sup>* + *φ* ()*<sup>s</sup>*

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

<sup>+</sup> *<sup>φ</sup>* (*ρcp*)*<sup>s</sup>*

<sup>+</sup> \_\_\_\_\_\_\_ <sup>1</sup> <sup>1</sup> <sup>+</sup> *<sup>φ</sup>* \_\_\_\_ 1 − *φ ρ*\_\_*s ρf* ]

*nf* , thermal expansion factor *βnf* , and thermal conduction *knf* , where

the volume of fluid. Relationships between

(<sup>1</sup> <sup>−</sup> *<sup>φ</sup>*)2.5, (3)

, (2)

. (4)

, (5)

. (6)

(1)

position in porous media.

82 Numerical Simulations in Engineering and Science

ity *μnf* , heat capacitance (*cp*

where *Vs*

)

*<sup>φ</sup>* <sup>=</sup> *<sup>V</sup>* \_\_\_\_\_*<sup>s</sup>*

Density of the nanofluid *ρnf* is given as:

of nanofluid can be given as:

(*ρcp*)*nf*

*βnf* = *β<sup>f</sup>*

is volume of solid particles and *Vf*

*ρnf* = (1 − *φ*) *ρ<sup>f</sup>* + *φ ρ<sup>s</sup>*

viscosity *μnf* can be given with according to [43] as:

*μnf* <sup>=</sup> *<sup>μ</sup>* \_\_\_\_\_\_ *<sup>f</sup>*

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

$$\frac{1}{\Phi} \frac{\partial \bar{\boldsymbol{\upsilon}}}{\partial t} + \frac{1}{\phi} \{\bar{\boldsymbol{\upsilon}} \cdot \bar{\boldsymbol{\mathfrak}}\} \boldsymbol{\bar{\upsilon}} = \frac{1}{\rho\_{\boldsymbol{\eta}\boldsymbol{\prime}}} \bar{\boldsymbol{\chi}} \boldsymbol{p} - \beta\_{\boldsymbol{\eta}\boldsymbol{\prime}} (\boldsymbol{T} - \boldsymbol{T}\_{\boldsymbol{\eta}}) \, \bar{\boldsymbol{\varrho}} + \frac{1}{\phi} \frac{\mu\_{\boldsymbol{\eta}\boldsymbol{\prime}}}{\rho\_{\boldsymbol{\eta}\boldsymbol{\prime}}} \nabla^{2} \boldsymbol{\bar{\upsilon}} - \frac{1}{K} \frac{\mu\_{\boldsymbol{\eta}\boldsymbol{\prime}}}{\rho\_{\boldsymbol{\eta}\boldsymbol{\prime}}} \boldsymbol{\bar{\boldsymbol{\upsilon}}} - \frac{F\_{\bar{\boldsymbol{\mathcal{U}}}} |\boldsymbol{\bar{\boldsymbol{\upsilon}}}|}{K^{32}}\tag{10}$$

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 = \frac{\phi^3 d\_p^2}{a \left(1 - \phi\right)^{2\nu}} F = \frac{b}{\sqrt{a \phi^{3\nu}}} \tag{11}$$

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:

$$
\sigma \frac{\partial T}{\partial t} + \left(\bar{\mathbf{q}} \cdot \bar{\mathbf{p}}\right) T = \frac{k\_\star}{\left(\mathcal{P}^\circ \mathbb{C}\_\eta\right)\_{\eta\sharp}} \nabla^2 T,\tag{12}
$$

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

Simulation of Natural Convection in Porous Media by Boundary Element Method

http://dx.doi.org/10.5772/intechopen.71230

85

) *f* /*μf kf* ;

The results are presented in terms of the porous Rayleigh number *RaP*, which links the ther-

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

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.

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

**Step 2**. Vorticity values on the boundary are calculated by a fast single-domain

The nondimensional parameters appearing in the momentum equation are:

are shown in the present chapter.

• Prandtl number Pr = *μf cp*

• Darcy number *Da* = *K*/*L*<sup>2</sup>

**2.2. Boundary element method**

• *RaP* = *RaT* ∙ *Da*.

• Fluid Rayleigh number *RaT* = *g β<sup>T</sup>* Δ*T L*<sup>3</sup> *ρ<sup>f</sup>* (*ρc<sup>p</sup>*

mal Rayleigh number and Darcy Number:

/*kf* ;

;

derivatives to calculate boundary vorticity values.

The outline of the numerical algorithm is as follows:

BEM from the kinematics equation (14).

where *σ* is the specific heat ratio, *σ* = *ϕ* + (1 − *ϕ*)(*ρc<sup>p</sup>* )*p* /(*ρc<sup>p</sup>* )*nf* , (*ρc<sup>p</sup>* )*p* and (*ρc<sup>p</sup>* )*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:

 $\bar{v} \to \frac{\bar{v}}{\bar{\nu}\_{0}}$   $\bar{r} \to \frac{\bar{r}}{L}$   $\ t \to \frac{v\_{0}t}{L}$   $\bar{g} \to \frac{\bar{g}}{\mathcal{S}\_{0}}$   $\ p \to \frac{p}{\bar{\nu}\_{0}}$   $T \to \frac{(T - T\_{0})}{\Delta T}$ 

The parameters in the above expressions are *v*<sup>0</sup> characteristic velocity given as *v*<sup>0</sup>  = *kf* /(*ρc<sup>p</sup>* ) *<sup>f</sup>L*, which is common choice for buoyant flow simulations, *kf* is the fluid thermal conductivity, (*ρc<sup>p</sup>* ) *f* is the heat capacity for the fluid phase, and *L* is the characteristic length. Moreover, *T*<sup>0</sup> is characteristic temperature *T*<sup>0</sup>  = (*T*<sup>2</sup>  − *T*<sup>1</sup> )/2 and Δ*T* is characteristic temperature difference Δ*T* = *T*<sup>2</sup>  − *T*<sup>1</sup> , *p*0 is the characteristic pressure *p*<sup>0</sup>  = 1bar, while gravitational acceleration is *g*<sup>0</sup>  = 9.81 m/s<sup>2</sup> .

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:

$$
\nabla^2 \bar{\upsilon} + \dot{\nabla} \times \bar{\omega}\_{\dot{\omega}} = 0 \tag{14}
$$

$$\left(\bar{\boldsymbol{\omega}} \cdot \boldsymbol{\nabla}\right) \bar{\boldsymbol{\omega}} = \left(\bar{\boldsymbol{\omega}} \cdot \boldsymbol{\nabla}\right) \bar{\boldsymbol{\upsilon}} - \boldsymbol{\mathsf{C}}\_{\boldsymbol{A}} \operatorname{Pr} \mathrm{R} a\_{\Gamma} \boldsymbol{\phi}^{2} \,\nabla \times \boldsymbol{T} \,\bar{\boldsymbol{\varrho}} + \boldsymbol{\mathsf{C}}\_{\mathbf{z}} \operatorname{Pr} \boldsymbol{\phi} \,\nabla^{2} \bar{\boldsymbol{\omega}} - \boldsymbol{\mathsf{C}}\_{\mathbf{z}} \frac{\mathrm{Pr}}{\mathrm{D} a} \boldsymbol{\phi}^{2} \,\bar{\boldsymbol{\omega}} - \frac{\mathrm{F}}{\mathrm{D} a} \boldsymbol{\phi}^{2} \,\left| \bar{\boldsymbol{\upsilon}} \right| \,\bar{\boldsymbol{\omega}},\tag{15}$$

$$(\dot{\psi} \cdot \dot{\nabla})T = \mathbb{C}\_{\mathbb{C}} \,\nabla^2 T\tag{16}$$

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

$$\mathbf{C}\_{A} = \frac{\mu\_{\eta'}}{\mu\_{f}} \frac{\rho\_{f}}{\rho\_{\eta'}} , \mathbf{C}\_{\mathbf{a}} = \frac{\beta\_{\eta'}}{\beta\_{f}} , \mathbf{C}\_{\mathbf{c}} = \frac{a\_{\eta'}}{a\_{\eta'}} \tag{17}$$

where *αnf* is thermal diffusivity of nanofluid, *αnf* = *knf* /(*ρc<sup>p</sup>* ) *nf* and *α<sup>f</sup>* thermal diffusivity of pure fluid *α<sup>f</sup>*  = *kf* /(*ρc<sup>p</sup>* )*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:


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

• *RaP* = *RaT* ∙ *Da*.

where *a* and *b* are Ergun's constants with values *a* = 150 and *b* = 1.75, while *dp*

<sup>∂</sup>*<sup>t</sup>* <sup>+</sup> (*v*⃑ <sup>∙</sup>

*<sup>L</sup>*, *<sup>t</sup>* <sup>→</sup> *<sup>v</sup>*<sup>0</sup> *<sup>t</sup>* \_\_\_

is the heat capacity for the fluid phase, and *L* is the characteristic length. Moreover, *T*<sup>0</sup>

∇

⃑ )*<sup>T</sup>* <sup>=</sup> *<sup>k</sup>* \_\_\_\_\_ *<sup>e</sup>* (*<sup>ρ</sup> cp*)*nf* ∇2

porous medium. It is assumed that the thermal properties of solid matrix and the nanofluid

 = *knf* . Governing equations (9), (10), and (12) can be converted into a nondimensional form by intro-

*<sup>L</sup>* , *<sup>g</sup>*⃑ <sup>→</sup> *<sup>g</sup>*⃑ \_\_

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:

*g*0

 = 1bar, while gravitational acceleration is *g*<sup>0</sup>

⃑ × *T g*⃑ + *CB Prϕ* ∇<sup>2</sup> *ω*⃑ − *CB*

, *CC* <sup>=</sup> *<sup>α</sup>nf* \_\_\_ *αf*

> ) *nf* and *α<sup>f</sup>*

. The nanofluid properties are obtained using the expressions (2)–(8). For the

⃑ )*T* = *CC* ∇<sup>2</sup>

In the above equations, parameters *CA*, *CB*, and *CC* are presenting the nanofluid properties and

, *CB* <sup>=</sup> *<sup>β</sup>nf* \_\_\_ *βf*

*μf ρf* \_\_\_ *ρnf* , *<sup>p</sup>* <sup>→</sup> *<sup>p</sup>*\_\_ *p*0

)/2 and Δ*T* is characteristic temperature difference Δ*T* = *T*<sup>2</sup>

)*p* /(*ρc<sup>p</sup>*

)*nf* , (*ρc<sup>p</sup>* )*p*

particle size of the bed.

*<sup>v</sup>*⃑ <sup>→</sup> *<sup>v</sup>*⃑ \_\_

teristic temperature *T*<sup>0</sup>

is the characteristic pressure *p*<sup>0</sup>

⃑ ) *ω*⃑ = (*ω*⃑ ∙ ∇

are given with expressions:

*p*0

(*v*⃑ ∙ ∇

fluid *α<sup>f</sup>*

 = *kf* /(*ρc<sup>p</sup>* )*f*

Finally, the energy equation can be written as:

where *σ* is the specific heat ratio, *σ* = *ϕ* + (1 − *ϕ*)(*ρc<sup>p</sup>*

duction of the following dimensionless variables:

, *<sup>r</sup>*⃑ <sup>→</sup> *<sup>r</sup>*⃑ \_\_

which is common choice for buoyant flow simulations, *kf*

are identical [24, 45], resulting in *σ* = 1 and *ke*

*v*0

The parameters in the above expressions are *v*<sup>0</sup>

 = (*T*<sup>2</sup>  − *T*<sup>1</sup>

∇<sup>2</sup> *v*⃑ + ∇

(*v*⃑ ∙ ∇

*CA* <sup>=</sup> *μnf* \_\_\_

where *αnf* is thermal diffusivity of nanofluid, *αnf* = *knf* /(*ρc<sup>p</sup>*

⃑ ) *v*⃑ − *CA PrRaT ϕ*<sup>2</sup> ∇

tances of solid and fluid phase, respectively. Furthermore, *ke*

*σ* \_\_\_ <sup>∂</sup>*<sup>T</sup>*

84 Numerical Simulations in Engineering and Science

is the average

)*nf* are heat capaci-

<sup>∆</sup>*<sup>T</sup>* (13)

 = *kf* /(*ρc<sup>p</sup>* ) *f L*,

 = 9.81 m/s<sup>2</sup>

) *f*

is charac-

.


 − *T*<sup>1</sup> ,

*T*, (12)

is the effective conductivity of

and (*ρc<sup>p</sup>*

, *<sup>T</sup>* <sup>→</sup> (*<sup>T</sup>* <sup>−</sup> *<sup>T</sup>*0) \_\_\_\_\_

⃑ × *ω*⃑ = 0 (14)

*Da <sup>ϕ</sup>*<sup>2</sup> *<sup>ω</sup>*⃑ <sup>−</sup> \_\_\_*<sup>F</sup>*

*T* (16)

*Da ϕ*<sup>2</sup> |*v*⃑

, (17)

thermal diffusivity of pure

\_\_\_ *Pr*

is the fluid thermal conductivity, (*ρc<sup>p</sup>*

characteristic velocity given as *v*<sup>0</sup>

### **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:

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

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

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

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

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

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

$$\mathcal{L}\left(\bar{\xi}\right)\dot{\boldsymbol{n}}\{\bar{\xi}\} \times \dot{\boldsymbol{v}}\{\bar{\xi}\} + \boldsymbol{n}\{\bar{\xi}\} \times \mathbf{f}\_{\mathbb{L}} \,\bar{\boldsymbol{v}} \,\nabla \boldsymbol{u}^\* \cdot \boldsymbol{n} \,d\boldsymbol{\varGamma} = \boldsymbol{n}\{\bar{\xi}\} \times \mathbf{f}\_{\mathbb{L}} \,\boldsymbol{v} \times \left(\boldsymbol{n} \times \nabla\right) \boldsymbol{u}^\* \,d\boldsymbol{\varGamma} + \boldsymbol{n}\{\bar{\xi}\} \times \mathbf{f}\_{\Omega} \{\boldsymbol{\omega} \times \nabla \ u^\*\} d\Omega \quad \text{(18)}$$

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*<sup>∗</sup> is the fundamental solution of the Laplace equation given as:

$$
\mu^{\*} = \frac{1}{4\pi \left| \bar{\xi} - r \right|} \tag{19}
$$

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

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

*φ* |⇀*v* |∫ Ω *Pr* \_\_\_1 *CB* \_\_1 *φ* ∫ Ω(

*ω<sup>j</sup> u*\* *d*Ω ,

Simulation of Natural Convection in Porous Media by Boundary Element Method

*<sup>Γ</sup> n*⃑ ∙ {*u*∗(*v*⃑*T*)} *d* − ∫

is a component of vorticity flux, whereas *qT* is a heat flux. In the sub-

*C*\_\_\_*A CB φ* ∫ Γ(*T* ∇ ⇀

⇀*<sup>v</sup> <sup>ω</sup><sup>j</sup>* <sup>−</sup> ⇀*<sup>ω</sup> vj*) <sup>⋅</sup> <sup>∇</sup>

http://dx.doi.org/10.5772/intechopen.71230

<sup>×</sup> *<sup>u</sup>*\*⇀*<sup>g</sup>* )*<sup>j</sup> <sup>d</sup>*<sup>Ω</sup>

*<sup>Ω</sup>*(*v*⃑*T*) ∙ ∇

⇀ *u*\* *d*Ω

⃑ *u*<sup>∗</sup> *d*] (22)

(21)

87

⇀*<sup>v</sup> <sup>ω</sup><sup>j</sup>* <sup>−</sup> ⇀*<sup>ω</sup> vj*)} *<sup>d</sup>*<sup>Γ</sup> <sup>−</sup> \_\_<sup>1</sup>

*CC* [∫

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.

The physical models where the above developed numerical scheme was tested are a twodimensional 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

(*u*\* *<sup>T</sup>*⇀*<sup>g</sup>* <sup>×</sup> ⇀*n*)*<sup>j</sup> <sup>d</sup>*<sup>Γ</sup> <sup>−</sup> *<sup>R</sup> aT*

*<sup>Γ</sup> <sup>u</sup>*<sup>∗</sup> *qT <sup>d</sup>* <sup>+</sup> \_\_\_<sup>1</sup>

from the previous iteration.

transport equation is:

*c*( *ξ* ⇀)*ωj*( *<sup>ξ</sup>*

*c*(*ξ*⃑

**3. Test cases**

*2.2.2. Integral representation of the vorticity and energy equations*

*u*\* ⋅ ⇀*n d*Γ = ∫

⇀*n* ⋅ {*u*\* (

⃑ *u*<sup>∗</sup> ∙ *n*⃑ *d* = ∫

**Figure 1**. The boundary conditions for the current problem are:

*C*\_\_\_*A CB φ* ∫ Γ Γ *u*\* *qj d*Γ

*<sup>ω</sup><sup>j</sup> <sup>u</sup>*\* *<sup>d</sup>*<sup>Ω</sup> <sup>+</sup> \_\_\_\_\_ *<sup>F</sup>* Pr √ \_\_\_ *Da* \_\_\_1 *CB*

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

*Da φ* ∫ Ω

⇀) <sup>+</sup> <sup>∫</sup> Γ *ω<sup>j</sup>* ∇ ⇀

 + \_\_<sup>1</sup> *Pr* \_\_\_1 *CB* \_\_1 *φ* ∫ Γ

    − *R aT*

    + \_\_\_<sup>1</sup>

)*T*(*ξ*⃑ ) + ∫ *<sup>Γ</sup> T* ∇

In the above equations, *qj*

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:

$$\mathcal{L}(\boldsymbol{\xi}) \, \_{\bar{\mathcal{U}}} (\boldsymbol{\xi}) + \int\_{I} \bar{\boldsymbol{\varphi}} \{ \bar{\boldsymbol{u}} \cdot \bar{\boldsymbol{\nabla}} \} \, \boldsymbol{u}^{\*} \, d\boldsymbol{\Gamma} = \int\_{I} \bar{\boldsymbol{\mathcal{U}}} \boldsymbol{\gamma} \times \{ \bar{\boldsymbol{u}} \times \bar{\boldsymbol{\nabla}} \} \, \boldsymbol{u}^{\*} \, d\boldsymbol{\Gamma} + \int\_{\bar{\boldsymbol{\Omega}}} \{ \bar{\boldsymbol{\omega}} \times \bar{\boldsymbol{\nabla}} \} \, \boldsymbol{u}^{\*} \, d\boldsymbol{\Omega}.\tag{20}$$

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*( *<sup>ξ</sup>* ⇀) <sup>+</sup> <sup>∫</sup> Γ *ω<sup>j</sup>* ∇ ⇀ *u*\* ⋅ ⇀*n d*Γ = ∫ Γ *u*\* *qj d*Γ + \_\_<sup>1</sup> *Pr* \_\_\_1 *CB* \_\_1 *φ* ∫ Γ ⇀*n* ⋅ {*u*\* ( ⇀*<sup>v</sup> <sup>ω</sup><sup>j</sup>* <sup>−</sup> ⇀*<sup>ω</sup> vj*)} *<sup>d</sup>*<sup>Γ</sup> <sup>−</sup> \_\_<sup>1</sup> *Pr* \_\_\_1 *CB* \_\_1 *φ* ∫ Ω( ⇀*<sup>v</sup> <sup>ω</sup><sup>j</sup>* <sup>−</sup> ⇀*<sup>ω</sup> vj*) <sup>⋅</sup> <sup>∇</sup> ⇀ *u*\* *d*Ω     − *R aT C*\_\_\_*A CB φ* ∫ Γ(*u*\* *<sup>T</sup>*⇀*<sup>g</sup>* <sup>×</sup> ⇀*n*)*<sup>j</sup> <sup>d</sup>*<sup>Γ</sup> <sup>−</sup> *<sup>R</sup> aT C*\_\_\_*A CB φ* ∫ Γ(*T* ∇ ⇀ <sup>×</sup> *<sup>u</sup>*\*⇀*<sup>g</sup>* )*<sup>j</sup> <sup>d</sup>*<sup>Ω</sup>     + \_\_\_<sup>1</sup> *Da φ* ∫ Ω *<sup>ω</sup><sup>j</sup> <sup>u</sup>*\* *<sup>d</sup>*<sup>Ω</sup> <sup>+</sup> \_\_\_\_\_ *<sup>F</sup>* Pr √ \_\_\_ *Da* \_\_\_1 *CB φ* |⇀*v* |∫ Ω *ω<sup>j</sup> u*\* *d*Ω , (21)

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

$$\mathcal{L}(\xi)T(\xi) + \mathfrak{f}\_{\Gamma}T \, \bar{\nabla} \, u^\* \cdot \mathfrak{n} \, d\Gamma = \mathfrak{f}\_{\Gamma} \boldsymbol{u}^\* \, q\_{\Gamma} \, d\Gamma + \frac{1}{\mathfrak{C}\_{\emptyset}} [\![\![ \boldsymbol{f}\_{\Gamma} \boldsymbol{\hbar} \cdot \{\boldsymbol{u}^\* \boldsymbol{\prime} \boldsymbol{T} \} \boldsymbol{\hbar} \, \Gamma - \not{f}\_{\Omega} \langle \boldsymbol{\hbar} \cdot \boldsymbol{T} \rangle \cdot \bar{\nabla} \, \boldsymbol{u}^\* \, d\Omega] \tag{22}$$

In the above equations, *qj* is a component of vorticity flux, whereas *qT* is a heat flux. In the subdomain 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**

**Step 3**. Velocity values within the domain are calculated by a sub-domain BEM

**Step 4**. Temperature values are calculated by a subdomain BEM from the energy

**Step 5**. Vorticity values within the domain are calculated by a subdomain BEM

**Step 6**. Convergence check; all steps from 2 until 5 are repeated until all flow fields

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

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:

) × ∫

where Ω is the computational domain and Γ = ∂Ω is the boundary of the domain, *c*(*ξ*⃑

<sup>Γ</sup> *v*⃑ × (*n*⃑ × ∇⃑ ) *u*<sup>∗</sup> *d*Γ + *n*⃑(*ξ*⃑

) = *θ*/4*π*, *θ* is the inner angle with origin in *ξ*⃑

lies on a smooth boundary, then *c*(*ξ*⃑

4*π*|*ξ* ⃑ − *r*⃑

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

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:

*<sup>Γ</sup> v*⃑ × (*n*⃑ ×

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

∇

⃑ ) *u*<sup>∗</sup> *d* + ∫

*<sup>Ω</sup>*(*ω*⃑ × ∇

) × ∫ Ω

<sup>|</sup> (19)

is the fundamental solution of the Laplace equation

(*ω*⃑ × ∇⃑ *u*∗)*d*Ω (18)

. If *ξ*⃑

)=1/2. Furthermore, *n*⃑ is

⃑ *u*∗)*d*. (20)

) is

lies inside

⃑ *u*<sup>∗</sup> ∙ *n*⃑ *d* = *n*⃑(*ξ*⃑

from the kinematics equation (14).

86 Numerical Simulations in Engineering and Science

from the vorticity equation (17).

achieve the required accuracy.

*2.2.1. Integral representation of the kinematics equation*

) × ∫ *Γ v*⃑∇

)=1, and if *ξ*⃑

*<sup>u</sup>*<sup>∗</sup> <sup>=</sup> \_\_\_\_\_\_\_ <sup>1</sup>

are reduced with a kernel expansion approximation technique [49].

⃑ ) *u*<sup>∗</sup> *d* = ∫

*<sup>Γ</sup> v*⃑(*n*⃑ ∙ ∇

geometric factor defined as *c*(*ξ*⃑

a vector normal to the boundary, and *u*<sup>∗</sup>

the domain, then *c*(*ξ*⃑

*c*(*ξ*) *v*⃑(*ξ*) + ∫

equation (16).

 *c*(*ξ* ⃑ ) *n*⃑(*ξ* ⃑ ) × *v*⃑(*ξ* ⃑ ) + *n*⃑(*ξ*⃑

given as:

The physical models where the above developed numerical scheme was tested are a twodimensional 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:

**Figure 1.** Two-dimensional and three-dimensional enclosures with corresponding boundary conditions.

$$\begin{aligned} \boldsymbol{v}\_{\times} &= \boldsymbol{v}\_{\times} = \boldsymbol{0}, \; \boldsymbol{\omega} = \boldsymbol{0}, \; T &= T\_{\rm H} \quad \text{at} \quad \mathbf{x} = \mathbf{0}, \\ \boldsymbol{v}\_{\times} &= \boldsymbol{v}\_{\times} = \boldsymbol{0}, \; \boldsymbol{\omega} = \mathbf{0}, \; T &= T\_{\rm c} \quad \text{at} \quad \mathbf{x} = L, \\ \boldsymbol{v}\_{\times} &= \boldsymbol{v}\_{\times} = \boldsymbol{0}, \; \boldsymbol{\omega} = \mathbf{0}, \; \frac{\partial T}{\partial \boldsymbol{y}} = \mathbf{0} \quad \text{at} \quad \mathbf{y} = \mathbf{0} \; and \; \boldsymbol{z} = \mathbf{0}, \\ \boldsymbol{v}\_{\times} &= \boldsymbol{v}\_{\times} = \boldsymbol{0}, \; \boldsymbol{\omega} = \mathbf{0}, \; \frac{\partial T}{\partial \boldsymbol{y}} = \mathbf{0} \quad \text{at} \quad \mathbf{y} = L \; \text{and} \; \boldsymbol{z} = L \end{aligned} \tag{23}$$

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 = \frac{k\_{\eta'}}{k\_{\uparrow}} \int\_{\Gamma} \bar{\nabla} \cdot \mathbf{T} \cdot \mathbf{n} \, d\Gamma\_{\prime} \tag{24}$$

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

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

**RaP = 100, Pr = 0.71, φ = 0.8**

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

**]** *K* **[W/m K]** *β* **[×10−5 K−1]** *α* **[×10−7 m2**

Simulation of Natural Convection in Porous Media by Boundary Element Method

**Da 10−1 10−2 10−3 10−4 10−5**

**/s]**

89

http://dx.doi.org/10.5772/intechopen.71230

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 nanopar-

From the comparison, it can be observed that the results agree well with the data from the

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

published studies, which confirm accuracy of the obtained numerical algorithm.

chosen.

*cp*

**Mesh Number of** 

 **[J/kg K]** *ρ* **[kg/m3**

**Table 1.** Thermophysical properties of water and Cu nanoparticles.

**nodes**

Water 4179 997.1 0.613 21 1.47 Cu 385 8933 400 1.67 1163

ticles were added to the water as a base fluid.

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

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


**Table 1.** Thermophysical properties of water and Cu nanoparticles.

*vx* = *vy* = 0, *ω* = 0, *T* = *TH at x* = 0,

**Figure 1.** Two-dimensional and three-dimensional enclosures with corresponding boundary conditions.

*vx* = *vy* = 0, *ω* = 0, *T* = *TC at x* = *L*, *<sup>v</sup> <sup>x</sup>* <sup>=</sup> *vy* <sup>=</sup> <sup>0</sup>, *<sup>ω</sup>* <sup>=</sup> <sup>0</sup>, \_\_\_ <sup>∂</sup>*<sup>T</sup>*

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

> *kf* ∫ *<sup>Γ</sup>* ∇

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

In the present study, the Cu nanoparticles are added to the water as a base fluid. The thermo-

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

*vx* <sup>=</sup> *vy* <sup>=</sup> <sup>0</sup>, *<sup>ω</sup>* <sup>=</sup> <sup>0</sup>, \_\_\_ <sup>∂</sup>*<sup>T</sup>*

is expressed in terms of the dimensionless Nusselt number as:

 = 1.

physical properties of Cu nanoparticles and water are given in **Table 1** [50].

*Nu* <sup>=</sup> *knf* \_\_\_

88 Numerical Simulations in Engineering and Science

of *Da*. The results are presented in **Table 2**.

thermal conductivities is *knf* /*kf*

<sup>∂</sup>*<sup>y</sup>* <sup>=</sup> <sup>0</sup> *at <sup>y</sup>* <sup>=</sup> <sup>0</sup> *and <sup>z</sup>* <sup>=</sup> <sup>0</sup>,

<sup>∂</sup>*<sup>y</sup>* <sup>=</sup> <sup>0</sup> *at <sup>y</sup>* <sup>=</sup> *<sup>L</sup> and <sup>z</sup>* <sup>=</sup> *<sup>L</sup>*

⃑ *T* ∙ *n*⃑ *d*, (24)

(23)


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

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


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

*Da*, the flow regime is transited into the Darcy flow regime, which can be described by the

**Figure 2.** Temperature fields for different values of *RaP* for Cu nanofluid and different solid volume fractions for Pr = 6.2,

Simulation of Natural Convection in Porous Media by Boundary Element Method

http://dx.doi.org/10.5772/intechopen.71230

91

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

**Figure 3.** Nusselt number values depending on porous Rayleigh number for Pr = 6.2, Da = 10−6 and different values of

Dracy's law.

solid volume fraction of nanoparticles.


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


**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).

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 Simulation of Natural Convection in Porous Media by Boundary Element Method http://dx.doi.org/10.5772/intechopen.71230 91

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

*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

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

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

**Table 5.** Nusselt number values for a natural convection in porous media saturated with nanofluid in 2D enclosure for

**Da RaP [24] Present [24] Present [24] Present**

**Da** *Ra<sup>P</sup>* **[51] [24] Present** 10−2 10 1.015 1.010 1.012

10−4 10 1.071 1.065 1.070

10−6 10 1.079 1.072 1.093

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

saturated with a pure fluid for Pr = 1.0, *φ*=0.6 and different Da and *RaP*, for 2D geometry.

100 1.530 1.533 1.503 1000 3.555 3.602 3.499

100 2.725 2.764 2.777 1000 9.183 9.454 9.174

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

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

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

various governing parameters (Pr = 6.2).

*RaP* **= 1000, Pr = 0.71,** *φ* **= 0.8**

90 Numerical Simulations in Engineering and Science

with pure fluid.

*ϕ* = 0.05 *φ* = 0.4 *φ* = 0.6 *φ* = 0.9

**Table 4.** Nusselt number values for the 3D natural convection in a cubic enclosure filled with porous media saturated

*φ* = 0.4 *ϕ* = 0.0 *ϕ* = 0.025 *ϕ* = 0.05

**Author details**

**References**

Janja Kramer Stajnko<sup>1</sup>

Maribor, Maribor, Slovenia

Springer; 2013

Transfer. 1986;**10**:557-570

and Mass Transfer. 2006;**49**:1430-1441

1972;**15**:1377-1383

\*, Renata Jecl1

Address all correspondence to: janja.kramer@um.si

and Jure Ravnik<sup>2</sup>

Simulation of Natural Convection in Porous Media by Boundary Element Method

http://dx.doi.org/10.5772/intechopen.71230

93

1 Faculty of Civil Engineering, Transportation Engineering and Architecture, University of

[1] Nield DA, Bejan A. Convection in Porous Media. 4th ed. United States of America:

[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] 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

[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] 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] Baytas AC, Pop I. Free convection in a square porous cavity using a thermal nonequilib-

[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

[8] Prasad V, Kulacki FA. Natural convection in horizontal porous layers with localized

[9] Kladias N, Prasad V. Natural convection in horizontal porous layers: Effects of Darcy

[10] Rosenberg ND, Spera FJ. Thermohaline convection in a porous medium heated from

[11] Storesletten L. Natural convection in a horizontal porous layer with anisotropic thermal

[12] Beck JL. Convection in a box of porous material saturated with fluid. Physics of Fluids.

below. International Journal of Heat and Mass Transfer. 1992;**35**:1261-1273

rium model. International Journal of Thermal Sciences. 2002;**41**:861-870

heating from below. Journal of Heat Transfer. 1987;**109**:795-798

and Prandtl numbers. Journal of Heat Transfer. 1989;**111**:926-935

diffusivity. Transport in Porous Media. 1993;**12**:19-29

2 Faculty of Mechanical Engineering, University of Maribor, Maribor, Slovenia

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

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