Thermophysical properties of water and Cu nanoparticles.
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
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 . 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 .
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 .
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  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 , 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 .
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  for particle simulations and was later used for a wide variety of problems, e.g., for acceleration of the boundary integral Laplace equation by  and for coupling with BEM for the boundary matrices by . 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:
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 . 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:
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  as:
where the effective viscosity does not depend on the nanoparticle type. The heat capacitance of nanofluid can be given as:
The nanofluid thermal expansion coefficient can be written in a similar way:
taking into account the definition of ρnf, it follows:
The effective thermal conductivity knf is given with the Wasp model , 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 temperature variations and can be given with the Boussinesq approximation as:
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:
where is volume averaged velocity vector.
The momentum equation is also known as Brinkman-Forchheimer equation and reads as:
where φ is porosity, t time, p pressure, T temperature, 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 . 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 ():
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:
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:
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 . 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:
In the above equations, parameters CA, CB, and CC are presenting the nanofluid properties and are given with expressions:
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 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  and porous media flow simulations . 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 . 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 . 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 .
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:
where Ω is the computational domain and Γ = ∂Ω is the boundary of the domain, is geometric factor defined as , θ is the inner angle with origin in . If lies inside the domain, then =1, and if lies on a smooth boundary, then =1/2. Furthermore, is a vector normal to the boundary, and u∗ is the fundamental solution of the Laplace equation given as:
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 .
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:
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 . The final integral form of the vorticity transport equation is:
and finally the integral form of the energy transport equation reads as:
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:
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:
where Γ is the surface through which the heat flux is calculated and 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.
|cp [J/kg K]||ρ [kg/m3]||K [W/m K]||β [×10−5 K−1]||α [×10−7 m2/s]|
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|
|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|
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  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 . According to the comparable study, the Cu nanoparticles were added to the water as a base fluid.
|RaP = 1000, Pr = 0.71, φ = 0.8|
|ϕ = 0.05|
|φ = 0.4||φ = 0.6||φ = 0.9|
|φ = 0.4|
|ϕ = 0.0||ϕ = 0.025||ϕ = 0.05|
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.
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.