double diffusive natural convection in a horizontal porous layer with the boundary domain integral method
RENATA JECL, JANJA KRAMER and LEOPOLD ŠKERGET
About the authors
Renata Jecl
University of Maribor,
Faculty of Civil Engineering
Smetanova ulica 17, 2000 Maribor, Slovenia
E-mail: renata.jecl@uni-mb.si
Janja Kramer
University of Maribor,
Faculty of Civil Engineering
Smetanova ulica 17, 2000 Maribor, Slovenia
E-mail: janja.kramer@uni-mb.si
Leopold Škerget
University of Maribor,
Faculty of Mechanical Engineering
Smetanova ulica 17, 2000 Maribor, Slovenia
E-mail: leo@uni-mb.si
Abstract
We present the boundary-domain integral method, one of the numerical methods for solving the transport phenomena in porous media. The results for the case of double diffusive natural convection in a porous horizontal layer, which is fully saturated with an incompressible fluid, are obtained. Modified Navier-Stokes equations were used to describe the fluid motion in porous media in the form of conservation laws for mass, momentum, energy and species. Several results for different cases of double diffusive natural convection in a porous horizontal layer are presented and compared with some published studies in which calculations with other numerical methods were performed.
Keywords
porous media, boundary domain integral method, double diffusive natural convection, Darcy-Brinkman equation
1 introduction
The main purpose of this article is to present the Boundary Domain Integral Method as a numerical method for solving problems encountered with the flow through porous media. Fluid dynamics in porous media is an important topic in many branches of engineering and science, as well as in many fields of practical interest. It plays a fundamental role in ground-water hydrology and soil mechanics, as well as petroleum, sanitary and chemical engineering. New computational methods and techniques have allowed us to model and simulate various transport phenomena in porous media, which means our understanding of them is improving constantly. The aim of this work is to obtain a numerical solution for the governing equations describing the flow of a viscous incompressible fluid in a porous medium, using an appropriate extension of the boundary-element method. The numerical scheme was tested on a natural convection problem within a porous square cavity as well as in the porous layer, where different temperatures and concentration values are applied on the horizontal walls. The results for the different governing parameters (the Rayleigh number, the Darcy number, the buoyancy ratio and the Lewis number) are presented, commented on and compared with previously published work.
1.1 transport phenomena in porous media
Fluid transport phenomena in porous media refer to those processes related to the transport of fluid momentum, mass and heat, through the given porous media. These processes, which are encountered in many different branches of science and technology, are commonly the subject of theoretical treatments that are based on methods traditionally developed in classical fluid dynamics. Natural convection is the most commonly studied transport phenomena in porous media and also a process to which we pay full attention. There have been several reported studies dealing with natural convection due to thermal buoyancy forces, mainly because of its
ACTA GEOTECHNICA SLOVENICA, 2009/1 51.
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
importance in industrial and technological applications, such as geothermal energy, fibrous insulation, etc. Less attention, however, has been dedicated to the so-called double diffusive problems, where density gradients occur due to the effects of combined temperature and compositional buoyancy. There are some important engineering applications for this phenomenon, such as the transport of moisture in fibrous insulations or grain-storage insulations and the dispersion of contaminants through water-saturated soil.
Three main configurations are studied when dealing with double diffusive natural convection in an enclosure filled with porous media [1], [2]:
-	temperature and species concentration or their gradi ents are imposed horizontally along the enclosure, either aiding or opposing each other,
-	temperature and species concentration or their gradients are imposed vertically, either aiding or opposing each other,
-	temperature or its gradient is imposed vertically and species concentration or its gradient is imposed horizontally or vice versa.
The present analysis focuses on the use of the Boundary Element Method for solving problems of fluid flows in porous media driven by coupled thermal and solutal buoyancy forces. The Darcy-Brinkman formulation is used for modeling fluid flow in porous media, thus enabling a satisfactory non-slip boundary condition on those surfaces that bound the porous media.
1.2 boundary element method for fluid dynamics in porous media
Fluid flows in porous media have been studied both experimentally and theoretically over recent decades. Different numerical methods have been introduced to obtain the solutions for some transport phenomena, e.g., the finite-difference method (FDM), the finite element method (FEM), the finite volume method (FVM), as well as the boundary element method (BEM). ). The main comparative advantage of the BEM, the application of which requires the given partial differential equation to be mathematically transformed into the equivalent integral equation representation, which is later to be discretized, over the discrete approximative methods, is demonstrated in cases where this procedure results in boundary integral equations only. This turns out to be possible only for potential problems, e.g., inviscid fluid flow, heat conduction, etc. In general, this procedure results in boundary-domain integral equations and,
therefore, several techniques have been developed to extend the classical BEM. The dual reciprocity boundary element method (DRBEM) represents one of the possibilities for transforming domain integrals into a finite series of boundary integrals, see for example [3] and [4]. The key point of the DRBEM is an approximation of the field in the domain by a set of global approximation functions and the subsequent representation of the domain integrals of these global functions by boundary integrals. The discretization of the domain is only represented by grid points (poles of global approximation functions) in contrast to FDM meshes. However, the discretization of the geometry and the fields on the boundary is piecewise polygonal, which gives this method greater flexibility compared to the FDM methods in coping with boundary quantities. In the DRBEM all the calculations reduce to the evaluation of the boundary integrals only. Another more recent extension of the BEM is the so-called boundary domain integral method (BDIM), see [5], [6], [7] and [8]. Here, the integral equations are given in terms of variables on the integration boundaries, as well as within the integration domain. This situation arises when we are dealing with strongly nonlinear problems that have prevailing domain-based effects, for example, diffusion-convection problems. Navier-Stokes equations are commonly used as a framework for the solution of such problems, since they provide a mathematical model of the physical conservation laws of mass, momentum and energy. The velocity-vorticity formulation of these equations results in a computational decoupling of the kinematics and kinetics of the fluid motion from the pressure computation, see [9]. Since the pressure does not appear explicitly in the field functions' conservation equations, the difficulty associated with the computation of the boundary pressure values is avoided. The main advantage of the BDIM, compared to the classical domain-type numerical techniques, is that it offers an effective way of dealing with boundary conditions on the solid walls when solving the vorticity equation. Namely, the boundary vorticity in the BDIM is computed directly from the kinematic part of the computation and not through the use of some approximate formulae. One of the few drawbacks of the BDIM is the considerable CPU time and memory requirements, but they can be drastically reduced by the use of a subdomain technique, see [10]. Convection-dominated fluid flows suffer from numerical instabilities. In domain-type numerical techniques upwinding schemes of different orders are used to suppress such instabilities, while in the BDIM the problem can be avoided by the use of Green's functions of the appropriate linear differential operators, which results in a very stable and accurate numerical description of the coupled diffusion-convective problems. There
18. ACTA GeOTeCHNICA SLOVENICA, 2009/1
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
are no oscillations in the numerical solutions, which would have to be eliminated by using some artificial techniques, e.g., upwinding, artificial viscosity, as is the case of other approximation methods.
2 governing equations
Due to the general complexity of the fluid transport process in a porous medium, our work is based on a simplified mathematical model in which it is assumed that:
-	the solid phase is homogeneous, non-deformable, and does not interact chemically with respect to the fluid,
-	the fluid is single phase and Newtonian; its density does not depend on pressure variation, but only on the variation of the temperature,
-	the two average temperatures, Ts for the solid phase and Tf for the fluid phase, are assumed to be identical - the porous media is in thermodynamic equilibrium, which means it is described by a single equation for the average temperature T = Ts = Tf ,
-	no heat sources or sinks exist in the fluid-saturated porous media; the thermal radiation and viscous dissipation are negligible,
-	the porosity and permeability are assumed to be constant throughout the whole cavity,
-	the density of the fluid depends on the temperature and concentration variations and can be described by p = p0 [1 -[3t(T-T0)-fiC(C-C0)] ,where the subscript 0 refers to a reference state, (3T is the volumetric thermal expansion coefficient ¡3T = —1/p[dp/df\ , and ¡3C is the volumetric concentration expansion coefficient 13C =-1/ p [dp/ 8C]t .
The transport phenomenon in the porous media is described using modified Navier-Stokes equations in the form of conservation laws for mass, momentum, energy and species. The equations are written at the macroscopic level, derived by averaging the microscopic equations for the pure fluid over the porous representative elementary volume and expressed by continuity, momentum, energy and species equations [11]:
dvi dx,
0
dC dv,C
-+ ——:
dt dx.
dx.
dC i-
dx.
D-
(4)
— = 0, (1)
1 dv, 1 dvv --L +
0 dt 0 dx.
1 dp	v d
=---- + Fg,--v, +-
p„ dx,	K dx.
v
2 t"
d r	i dvT d
-0 Cf + (1 - 4K ] T + c,= -
dT dx.
, (2)
(3)
The parameters used above are as follows: v, volume-averaged velocity, x, the i-th coordinate, 0 porosity, t time, p density, v kinematic viscosity, dp/dxi the pressure gradient, g, gravity, K permeability of the porous media and è,, the strain rate tensor è.. = l/2 (dv¡dx + dv /dxt ) andfinally F is the normalized density difference function and is given as F = (P - Po V Po =[ (T -T0) + I3C (C - C0 )] .Furthermore, cf = (pc)f and cs = (pc)s are the isobaric specific heat per unit volume for the fluid and solid phases, respectively, T is the temperature, Ae is the effective thermal conductivity of the porous media given as Ae = 0Af + (l — 0)As, where Af and As are the thermal conductivities for the fluid and solid phases, respectively. In the final equation C stands for the concentration, and D for the mass diffusivity.
Equation (2) is known as the Darcy-Brinkman equation and consists of two viscous terms: the Darcy viscous term (third on the r.h.s.) and the Brinkman viscous term, which is analogous to the Laplacian term in the Navier-Stokes equations for a pure fluid (fourth on the r.h.s.). A non-slip boundary condition on a surface that bounds the porous media is satisfied by using the additional Brinkman term. There are several situations where it is convenient to use the Brinkman equation, e.g., when one wishes to compare the flows in porous media with those in a pure fluid or to investigate the convective flows in the context of solidification, where permeability and porosity are variables in space and time. If the included parameter K (the permeability) tends towards infinity ( K ) the Brinkman equation transforms into the classical Navier-Stokes equation for a pure fluid. On the other hand, if the permeability tends to zero ( K ^ 0 ), the Brinkman term becomes negligible and the equation reduces to a classical Darcy equation [1].
A convective flow in a horizontal porous layer is possible above the critical Rayleigh number. In the case of double-diffusive convection, where the density differences are the result of combined temperature and concentration gradients, the critical Rayleigh number is a function of the Darcy number Da, the Lewis number Le and the buoyancy coefficient N [1].
In vertical cavities maintained to horizontal temperature and concentration gradients the flow in the cavity is always unicellular. In the case of a horizontal porous layer, where the temperature and concentration differences are imposed on the horizontal walls, the flow structure becomes multi-cellular and is also called a Rayleigh-Benard flow structure [2].
18. ACTA GeOTeCHNICA SLOVENICA, 2009/1
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
3 numerical method
The Boundary Domain Integral Method (BDIM), an extension of the classical Boundary Element Method (BEM), is used in order to solve the general set of equations. The discretization of the surface and the domain is required since the boundary and domain integrals are represented in the obtained set of integral equations. The above given equations (1), (2), (3), (4) should be modified in order to use the BDIM. Firstly, the kinematic viscosity in the momentum equation is partitioned into constant and variable parts, such as v = v + V, so the Brinkman term is also divided into two parts, and the equation becomes:
dv ' dv ' ,v ' 1 dp	vé
-L H--=---— + Fg¡--— v'¡-
dt dx	p0 dx¡	K
d2
dx ,dx.
-dx-, (5)
di¡ ' ~3t~

ew
?¡ikSk
dF dx,
dv ' ¡ vé ' dx, K
d I vdi'j + f öx, I dx, I dx,
, (10)
where w '¡ is the modified vorticity w \ = wjand fJ is any contribution arising as a result of non-linear material properties, given as fj = j/(V x è ).
The mathematical description of the transport phenomena in fluid flow is completed by providing suitable Dirichlet and Neumann boundary conditions [12], as well as some initial conditions for the energy transport and species transport equations (see Fig. 1):
-	dT _
T = T on the boundary r ; d— nj = q on the
boundary r2 ; T = T0 in the dornain Ü , (11)
where the term vj is now the modified velocity vi = vt . In the same way as the kinematic viscosity, the thermal diffusivity ap , which is defined as ap = \/cf , and the mass diffusivity D are divided into the constant and variable parts ap = ap + ap , D = D + D . By including the expression for the heat-capacity ratio a = $ + (1 — $)cs / cf in the energy equation, the formulations (3) and (4) can be rewritten as:
a dT_
é dt
dC_ dt ~
dvJT dx,
dv 'C
d T
dx,
dxjdxj
D d2c
é dxpx j
_d_ dx,
dx,
aP dT é dx,
D dC_
é dx,
(6)
(7)
-	dC _
C = C on the boundary r ; -n. = q on the
_	dx 1
boundary r2 ; C = C0 in the domain Ü , (12)
In the present analysis, a two-dimensional problem is considered, therefore all subsequent equations will be written for the case of planar geometry. The linear elliptical Laplace differential operator can be used for the velocity equation (9):
L\-\ =
d 2()
dxdx,
(13)
and the following relationship for the kinematics can be obtained:
In the next step the above-stated governing equations are transformed by the use of the velocity-vorticity formulation (VVF), and consequently the computational scheme is partitioned into its kinematic and kinetic parts [6]. The vorticity vector — is introduced, which
represents the curl of the velocity field: -=* % <8>
where ejjk is the unit permutation tensor. The kinematic part is represented by the ellipticalvelocity vector equation:
d2 v'
L iv'i+b=!;+b¡=0 -
(14)
where bi stands for the pseudo-body source term. An integral representation of the velocity vector can be formulated by using the Green theorems for scalar functions or the weighting residuals technique, rendering the following vector integral formulation:
c(&V.(Q + fv\q'dr= fdv-t-u'dT+ fbudQ , (15)
t	t	dn
d2v\ dxjdxj
yk
di\ dx,
= 0 ,
(9)
and the kinetics is governed by the vorticity, energy and species transport equation. The vorticity transport equation is obtained as a curl of the Brinkman momentum equation (5):
where £ is the source point, u is the elliptical Laplace fundamental solution and q is its normal derivative, e.g., q = du*/dn . The fundamental solution u* is given by the expression:
u = — ln 2^
r(£,s)
(16)
8. ACTA GeOïeCHNICA SLOVENICA, 2009/l
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
where r is the vector from the source point £ to the reference field point s. Equating the pseudo-body force as:
=tr - <17>
we obtain the following integral formulation:
«K (O+J v\qd r = J'dnru'd r+ev Sd7u'd ° m
Q 9Xj
The derivative of the vorticity in the pseudo-body source terms can be eliminated by using the Gauss divergence theorem. The integral representation of the kinematics is now [13]:
c(Z)«(!) X v) + n(i) xf v' qd r =	, (24)
1 = n(0 x f (q* x n) x v 'dr + n(Qx f u X q'dO r	o
in order to obtain an appropriate non-singular implicit system of equations for the unknown boundary vorticity or tangential velocity component values to the boundary
[13].
The formulations for the vorticity, temperature and concentration can generally be written as a non-homogeneous elliptical diffusion-convective equation [8]:
d2u dv'u
dxj dxj dxj At
-JL + bt = 0 , (25)
c(î)v\(î) + J v\qd r = J--u d r + et] J u' n]u*dY-ej J u' q*
r	r dn	r	a
,(19)
d a
By using these expressions for the vorticity definition, unit tangent, and normal vector, e.g.,
dv j/dn = dvj/OxjUj,
io' = e v dvi /dxi = dvy jdx — dvxjdy, n = (nx, ny) and t = (—ny, nx) for i, j = 1, 2 and applying the continuity equation (1), the following relationship can be derived:
+ ej^Vt" ' ^
The boundary integrals on the right-hand side of equation (19) can be rewritten as:
c(0v\(0 + f v'i qdT = —ev f^Vfu'dr — ei} f <j' qjda ,(21)
r	r dt	a
With further mathematical reformulations and by applying the Gauss theorem, the resulting integral representation can be written as follows:
c(Ov\(D + fv\q'dT = ejf v^q'dT — egf u'qjdQ ,(22)
r	r	a
where q* is the tangential derivative of the fundamental solution q* = du*/dt, or in the form of the elliptical integral vector formulation:
c(C)v '(£) + JVq'dr = J(q xn)xv'dr + JuXq*dQ ,(23)
The tangential form of equation (22) can be written as follows:
where u is taken as the vorticity u ', the temperature T and the concentration C, respectively, p is defined by considering the conservation laws and constitutive hypothesis always divided as p = p + p , and bt stands for the pseudo-body source term. Since the fundamental solution exists only for steady diffusion-convective Partial Differential Equations (PDEs) with constant coefficients, the velocity field is decomposed into an average constant vector vj ' and a variable vector Vi ', such that v\ = v V + V \. Thus, the following integral formulation can be obtained:
c(Ç)u(0 + pJ uq'd r = J(p q -uv'n)u*d r + J btud Q ,(26)
where u* is now the fundamental solution of the steady diffusion-convective PDE with a first-order reaction term [5], in the form of:
1
2nP
K o(pr )exp
2p
(27)
du *	du
Furthermore, q = — and q = — represents the
dn	dn
normal derivative of the fundamental solution and a suitable field function. The parameter ^ is defined as
V = i
— \2 v.
2 P
+ koL =
P
— \2 2P
+ 3 , where v'2 = v'.
/3 = 1/ p At, K0 is the modified Bessel function of the second kind of order 0 , and r is the magnitude of the vector from the source to the reference point, i.e., r = \xi(0 — xi. The following integral representations for the vorticity, temperature and concentration kinetics are obtained according to equation (26):
a
a
ACTA GeOTeCHNICA SLOVENICA, 2009/l 17.
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
) + f w 'Q'd r =1 f (q — w 'v'+ ejgjFnj + fjnj )U * d r +
r	v r
+1 f (w Vj — e jgjF— V qj — fj )Q*dn+1 f V^w U' dn + v V	v i K
— f wF ,u*dn
7At J F—1
vAt'
, (28)
c(0T(0 + f T Q'd r = fl a-^q — Tv'n
J	an J l m
r
I aaP
—=~ I \ ~irqj— Tv
a I m	■
aaP r
Qjd n
u dr—
p n
j ' aapAt n
f Tf—U'
, (29)
d n
)C(^) + f CQ*dr = D f \D q — cv
D-
—m flDqt—cv' Q*dn+^- fCF ,U' D J l m ^j j j D At J F—1
D-
u dr—
, (30)
D At-
d n
where Q* = dUYdn represents a normal and Q* = dU Ydej a space derivative of the modified elliptical diffusion-convective solution U* defined as U = vu in the momentum equation, U* = apjmu* in the energy equation and U* = D/mu* in the species equation. Furthermore, q = du/dn represents a normal and q. = du/dxj a space derivative of the field function (vorticity, temperature, concentration).
For an approximate numerical solution, the integral equations are written in a discretized manner, where the integrals over the boundary and the domain are approximated by the sum of the integrals over all the boundary elements and internal cells, respectively. The variation of the field functions within each boundary element or internal cell is approximated by the use of the appropriate interpolation polynomials [14]. After applying the discretized integral equations to the boundary and internal nodes, the implicit matrix systems can be obtained, firstly for the kinematic computational part as:
[H]{v',} = e, [H,][v
— e„ D». (31)
Furthermore, we obtain the implicit systems for the vorticity kinetics:
-1 [D,]{w V — ejgF — V qj— f} +1 [B]j V^w'} + ^ [B][w'r_,}
,(32)
then for the heat energy kinetics
[h]{t} = -a-[g]j^q — TVn [ —D l|a^q] — TV
1 Jl 1 aap j m	j aaP 1 j j[ m
m
,(33)
+-§- [b][T },—i
aapAt
and finally for the species kinetics:
[H ]{C} = t [G]j D q — C/}— t j D q; — CV' }+ Da [B][C}f—, .(34)
In the above equations the matrices [H], [Ht], [G], [Dj] and [5] are the influence matrices composed of those integrals that took over the individual boundary elements and internal cells. This system of discretized equations is solved by coupling these kinetic and kinematic equations and considering the corresponding boundary and initial conditions. Since the implicit set of equations is written simultaneously for all the boundary and internal nodes, this procedure results in a very large and fully-populated system matrix, influenced by diffusion and convection. The consequence of this approach is a very stable and accurate numerical scheme with substantial computer time and memory demands. The subdomain technique is used to improve the economics of the computation, where the entire computational domain is partitioned into subdomains, to which the same described numerical procedure can be applied. The
n
18. ACTA GeOTeCHNICA SLOVENICA, 2009/1
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
final system of equations for the entire domain is then obtained by adding the sets of equations for each subdomain and considering the compatibility and equilibrium conditions between their interfaces, resulting in a sparser matrix system, suitable for solving using iterative techniques. In the present case, each subdomain consists of four discontinuous 3-node quadratic boundary elements and one 9-node corner continuous quadratic internal cell [6].
3.1 solution procedure
The obtained set of discretized equations (31)-(34) is solved by coupling the kinetic and kinematic equations. In this way the solenoidality of the velocity field for an arbitrary vorticity distribution is assured. The following numerical algorithm has to be performed in order to obtain a final solution:
(1)	Start with some initial values for the vorticity distribution.
(2)	Kinematic computational part:
-	solve the implicit set for the boundary vorticity values (equation (31)),
-	determine the new domain velocity values,
(3)	Energy kinetic computation part:
-	solve the implicit set for the boundary and domain values (equation (33)),
(4)	Species kinetic computation part:
-	solve the implicit set for the boundary and domain values (equation (34),
(5)	Vorticity kinetic computational part:
-	solve the implicit set for the unknown boundary vorticity flux and internal domain vorticity values (equation (32)),
(6)	Relax new values and check the convergence. If the convergence criterion is satisfied, then stop; otherwise go to step 2.
4 test examples - results and discussion
The convective flow in a horizontal porous layer is possible above the critical Rayleigh number. In the case of double-diffusive convection, where the density differences are a result of combined temperature and concentration gradients, the critical Rayleigh number is a function of the Darcy number Da, the Lewis number Le and the buoyancy coefficient N [14]. In vertical cavities maintained at horizontal temperature and concentration gradients, the flow in the cavity is always
unicellular. In the case of a horizontal porous layer, where the temperature and concentration differences are imposed on the horizontal walls, the flow structure becomes multi-cellular and is also called a Rayleigh-Benard flow structure [15]. Most of the studies regarding double-diffusive convection or thermohaline convection (the case where the constituent is salt) in a horizontal porous layer are focused on the problem of convective instability. There are many studies dealing with the onset of convection on the basis of the linear stability theory [16], [17] or nonlinear perturbation theory [18]. In these studies, the critical Rayleigh numbers for the onset of convective flows are predicted. The theoretical and numerical study of heat and mass transfer affected by a high Rayleigh number Benard convection in a porous layer heated from below is obtained in [19]. The numerical results and a scale analysis of the flow in a porous medium are presented, where the buoyancy effect is due entirely to temperature gradients.
Some further numerical results for a double-diffusive convection in a horizontal porous layer with two opposing buoyancy sources can be found in [20]. The influence of the governing parameters (the Rayleigh number, the Lewis number, and the buoyancy ratio) on the overall heat and mass transfer is discussed for the case of a square cavity (aspect ratio equal to 1). Double-diffusive convection in a horizontal layer with some numerical results is also discussed in [21]. In this study the critical values of the Rayleigh numbers for the onset of convective motion are predicted on the basis of a nonlinear parallel flow approximation.
All the above-mentioned numerical results are obtained on the basis of the Darcy flow model, which is more convenient for porous media with a low permeability. The Brinkman extended Darcy model, which accounts for friction due to macroscopic shear is more appropriate when describing fluid flows in the porous matrix, when the inertia effects are negligible. It was used in [14] to investigate the onset and development of double-diffusive convection in a horizontal porous layer with uniform heat and mass fluxes specified at the horizontal boundaries. The obtained analytical solution is compared to some numerical results obtained for different values of Ra, N, Le and Da.
Only a few recently published studies deal with the cross relation between the temperature and the solutal gradient (the temperature or concentration are imposed horizontally and the concentration or temperature gradient imposed vertically) [22], [23], [24]. In those studies it is shown that the competition between the thermal and solutal force produces very complex flow patterns.
18. ACTA GeOTeCHNICA SLOVENICA, 2009/1
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
In our case the system under consideration is a horizontal layer of width L and height H, filled with homogenous, nondeformable porous media, which is fully saturated with a Newtonian fluid. In the first case, the horizontal walls are subjected to different temperature and concentration values (TB, CB at the bottom boundary and Tn, Cn at the upper boundary), while the vertical walls are adiabatic and impermeable (Figure 1). In addition, the second case - the example of cross gradients of the temperature and concentration - is presented in Figure 4, where the horizontal walls are maintained at different concentrations (CB at the bottom boundary and Cn at the upper boundary) and the vertical walls at different temperatures (TL at the left boundary and TR at the right boundary). The fluid, saturating the porous media, is modeled as a Boussinesq incompressible fluid, where the density depends only on the temperature and concentration variations P = Po(l- Pt (T - To) - Pc (C - Co)). Pt and ¡3C are the volumetric thermal and concentration expansion coefficients. The numerical results for the heat and mass transfer induced by the double-diffusive natural convection in a horizontal porous layer subjected to vertical gradients of temperature and concentration are presented.
The governing parameters of the problem are:
-	the porosity $ ,
-	the Darcy number Da = K / H2,
-	the aspect ratio A = L / H ,
-	the modified (porous) thermal Rayleigh number Ra = Kg PT AT H / a v , where PT is the thermal expansion coefficient, a the thermal diffusivity and v viscosity,
-	the Lewis number Le = a / D ,
-	the buoyancy ratio N = PC AC / PT AT ,
where K is the permeability, D and H are the width and height of the layer, respectively, g is the acceleration due to gravity, AT and AC are the temperature and concentration differences between the upper and lower boundaries, a is the thermal diffusivity and D is the mass diffusivity. For the aspect ratio A = 1, a non-uniform computational mesh 20x20 was used with a ratio between the longest and shortest elements of r = 6, and for A = 2 and A = 4 20x10 subdomains were used. The computational mesh was denser near the solid walls, where large gradients in the temperature, concentration and velocity profiles were expected. Time-steps ranging from At = 1016 s (steady state) to At = 10-4 s (time dependent) were employed, and the convergence criteria (representing the difference between two progressive values in each time step) were determined as e = 5 x 10-6 for all cases.
It should be noted that in the case of N = 0 the buoyancy effect is due entirely to temperature gradients. The mass transfer in this case is due to the temperature field and the concentration differences between the horizontal boundaries. In the case of positive values for the buoyancy ratio (N > 0), the thermal and solutal buoyancy forces aid each other (aiding convection) and for negative values of the buoyancy ratio (N < 0) the solutal and thermal effects have the opposite tendencies (opposing convection).
The obtained numerical model was tested for different values of the governing parameters. The results for total heat and mass transfer through the horizontal layer are given by the values of the Nusselt (Nu) and Sherwood numbers (Sh). Furthermore, some simulations of the temperature and concentration fields are presented.
Tu = o, a=o ,v, = o
H
Figure 1. Geometry of the first problem and the boundary conditions.
ACTA GeOTeCHNICA SLOVENICA, 2009/l 17.
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
The validation of the code was accomplished by a comparison with some published numerical experiments. However, it should be noted that there are not many published studies giving numerical results for the described problem.
Table 1 presents some results for A=1, Da = 10-5 and a different Rayleigh number, Lewis number and buoyancy ratio. The values of the overall heat and mass transfer are compared to the published results, where the numerical calculations based on the Darcy model are obtained [20]. The first result is for the case of Ra = 600, Le = 1 and N = 0 , which means that only the thermal buoyancy force is present. The overall heat and mass transfers, which are represented by Nu an Sh, are identical. The other two cases are for the Rayleigh number 100, the buoyancy ratio 0.2 and the Lewis numbers 10 and 30. In this case, both the thermal and solutal buoyancy forces are present and they aid each other. The values of the Sherwood numbers are now higher than those of the Nusselt numbers, which is a result of a higher Lewis number. The presented results are in agreement with the published ones.
Table 2 presents the influence of the Darcy number, where the governing parameters are as follows: an aspect ratio A = 1 , the Rayleigh number Ra = 100, the Lewis number Le = 10 and the buoyancy coefficient N = 0.2.
It is evident that with a decrease of the Darcy number the values of the Nusselt and Sherwood numbers increase. The Darcy number describes the influence of the additional Brinkman viscous term in the momentum equation. At higher values of the Darcy number, the effect of viscous forces is significant and generally slows down the convective motion. With a decrease of the Darcy number Da ^ 0 the viscous effect becomes smaller, so the convective motion becomes stronger, which can be observed from the presented results.
In Table 3 the results for the case of the opposing fluxes of heat and mass (N < 0) and the parameters A = 1 , Ra = 100, Le = 1, Le = 3 and N = -0.1 are presented. In the case where 10-3 < Da < 10-1 the values of the Nusselt and Sherwood numbers equal 1.0, which means the heat and solute transfer are governed by diffusion. The Rayleigh number, which depends on the Darcy number, is not high enough for the onset of convective motion, which begins in this from values of Da < 10 4. In the case of Le = 1 the values of the Nusselt and Sherwood numbers are identical. The Sherwood number increases with the Lewis number; however, on the other hand, the Nusselt number is almost independent of the Lewis number.
Table 1. Comparison of results with numerical experiments reported in the literature.
		Nu		Sh
	This study	Literature [20]	This study	Literature [20]
Ra=600, Le=1, N=0	7.01	6.6	7.01	-
Ra=100, Le=10, N=0.2	2.76	2.4	8.84	-
Ra=100, Le=30, N=0.2_250_25_14.80_15
Table 2. Nu and Sh numbers for different values of Da and A = 1 , Ra = 100, Le = 10, N = 0.2
Da	10-1	10-2	10-3	10-4	10-5	10-6
Nu	1.00	1.53	2.26	2.63	2.76	2.79
Sh	1.00	3.93	6.13	7.93	8.84	9.10
Table 3. Nu and Sh numbers for different values of Da and A = 1 , Ra = 100, N = -0.1
	Da	10-1	10-2	10-3	10-4	10-5	10-6
	Nu	1.00	1.00	1.00	2.30	2.40	2.42
	Sh	1.00	1.00	1.00	2.30	240	2.42
Le = 3
Nu
1.00
1.00
1.00
2.32
2.42
2.44
Sh
1.00
1.00
1.00
4.05
4.30
4.36
Table 4. Nu and Sh numbers for different values of Da and A = 4, Ra = 300, Le = 0.1, N = -2
Da	10-1	10-2	10-3	10-4	10-5
Nu	1.00	2.05	2.82	3.12	3.50
Sh	1.00	1.02	1.04	1.05	1.06
ACTA GeOTeCHNICA SLOVENICA, 2009/1 13 .
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
Table 4 presents the results for a shallow layer with the governing parameters as follows: the aspect ratio A = 4, the Rayleigh number Ra = 300, the Lewis number Le = 0.1 and the buoyancy ratio N = -2. In this case again the thermal and solutal buoyancy forces oppose each other. Because of the higher value of N, the combined buoyancy is dominated by the buoyancy due to concentration gradients. From the table it is evident that with any decrease in the Darcy number the value of the Nusselt number increases. In the cases where the Lewis number decreases Le ^ 0, the values of the Sherwood number tend to unity (Sh ^ 1), which implies that the mass transfer is dominated by diffusion. The same conclusions are also published in [14]. A direct comparison of the results is impossible, because of the different governing parameters (porosity, aspect ratio) in both studies.
The graphical simulations of the velocity, temperature and concentration fields are presented in Figure 2 for the square geometry (A = 1) with the parameters Ra = 100, Da = 10-4, Le = 10 and N = 0.2 in Figure 3 for the case where A = 2 for the same parameters.
The flow structure consists of a unicellular circulation filling up the entire layer. The temperature and concentration field present the classical stratified structures of the natural convective flows. The temperature and concentration boundaries are observed at the top and bottom walls.
The governing parameters for the second case are A = 2, Ra =100, Da = 10-3, Le = 10, and N = 1.
The flow in the horizontal layer where A = 2 becomes multi-cellular (with two cells) as it is obvious from Figure 3. The flow consists of rising hot fluid in the centre of the layer and colder fluid sinking along the vertical walls. In the centre of the domain, a higher solute concentration is found than along the adiabatic and impermeable side walls. Thin temperature and composition boundary layers are evident at the top and bottom walls.
Figure 3. Streamlines, isotherms and isoconcentrations in a horizontal layer for A = 2, Ra = 100, Da = 10-3, Le = 10, N=1.
The next example presents the convection under cross gradients of temperature and concentration. The geometry and boundary conditions are shown in Figure (4).
Figure 2. Streamlines, isotherms and isoconcentrations for A = 1, Ra = 100, Da = 10-4, Le = 10, N=0.2 .
18. ACTA GeOTeCHNICA SLOVENICA, 2009/1
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
TL = 1
— = o
dx u
v, =0
f)C
Cu=0 , f^=0 , v,=0
	porous media	TR = O
a y	g f	— = 0 dx u Vi = 0
X		
w	cB=i |Ç-O Vi = 0 ' dy ' L	
H
Figure 4. Geometry of the second problem and the boundary conditions.
Figure 5. Streamlines, isotherms and isoconcentrations for A = 2, Ra = 100, Da = 10-3, Le = 10 and N = 0 (left side) and N = -1.5 (right side).
18. ACTA GeOTeCHNICA SLOVENICA, 2009/1
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
In this case the temperature differences are imposed on vertical walls and the species concentration differences on horizontal walls.
Because a differentially heated cavity has a tendency to mix the fluid, while the solutal effect has a tendency to prevent such a mixing, it is expected that the competition between the solutal and thermal forces will produce complex flow patterns.
The graphical simulations for the streamlines, temperature and concentration fields are presented in Figures 5 and 6 for the parameters A = 2, Ra =100, Da = 10-3 and Le = 10 and the values of the buoyancy coefficient from N = 0 to N = -4 . In the case when N = 0 the flow remains unicellular, while with any decrease in the buoyancy coefficient the flow forms two slow cells. The temperature and concentration fields have typi-
cally stratified structures. With any decrease in N, the isotherms and concentration lines become more linear, which means the convective mechanism is suppressing. For the case when N = -4, heat and mass transfer take place mainly by conduction, which is evident from the temperature and concentration fields in Figure 6, where the isotherms and isoconcentration lines are almost parallel to the vertical and horizontal walls, respectively.
The values of the Nusselt and Sherwood numbers for the presented examples are given in Table 5. The results are compared to reference [24], where the same problem is solved for planar and spatial geometries with the use of the finite-difference method. Good agreement between the results can be obtained, which proves the accuracy of the mathematical model and also of the numerical scheme.
Figure 6. Streamlines, isotherms and isoconcentrations for A = 2, Ra = 100, Da = 10-3, Le = 10 and N = -2 (left side) and N = -4 (right side). l6. ACTA GeOTeCHNICA SLOVENICA, 2009/1
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
Table 5. Nu and Sh numbers for different values of N and A = 2, Ra = 100, Le = 10, Da = 10-3.
N		-4	-2	-1.5	-1	0
Nu	This study	0.49	0.50	0.50	0.52	0.81
	Literature [24]	-	0.50	0.50	0.55	-
Sh	This study	1.04	1.17	1.30	1.71	3.39
	Literature [24]	-	1.20	1.30	1.90	-
4 conclusion
A numerical approach, based on the boundary domain integral method (BDIM), which is an extension of the boundary element method (BEM), was applied for the solutions of the transport equations in porous media. The modified Navier-Stokes equations, Brinkman-extended Darcy formulation with the inertial term included, were employed to describe the fluid motion in porous media. The porous media is saturated with a viscous, incompressible fluid. A velocity-vorticity formulation of the governing equations is adopted, resulting in the computational decoupling of the kinematics and kinetics of the fluid motion from the pressure computation. Since the pressure does not appear explicitly in the field function conservation equations, the difficulty connected with the computation of boundary pressure values is avoided. The proposed numerical procedure is applied to the case of natural convection in a porous layer where the horizontal walls are subjected to different values of temperature and concentration and in a porous layer, where the horizontal walls are subjected to different values of concentration and the vertical walls to different values of temperature (convection under the cross gradients of temperature and concentration). The examples for different values of the modified Rayleigh number, the Darcy number, the Lewis number and buoyancy coefficient are presented and compared with the published studies. Very good agreement between the results obtained with alternative numerical methods (the finite difference method, the finite volume method, and the finite element method) can be observed. Since the final set of equations results in a very large and fully populated system matrix, the computer time and memory demands are large. The need to improve the economics of the computation, when using the boundary domain integral method, is still one of the challenges for researchers. One of the efficient mathematical tools, developed especially for saving computational time and computer storage with the boundary element method, is the wavelet transform [25].
However, it can be stated that the boundary domain integral method, extended in a way that also enables
an investigation of the fluid-transport phenomena in a porous media, appears to possess the potential to become a very powerful alternative to existing numerical methods, e.g., finite differences or finite elements, as a means to obtain solutions to the most complex systems of nonlinear partial differential equations, when attacking some unsolved problems in engineering practice.
references
[1]	Nield, D.A. and Bejan, A. (2006). Convection in Porous Media. 3rd.ed., Springer, Berlin.
[2]	Vafai, K. (2005). Handbook of Porous Media. Taylor&Francis, Boca Raton, London, New York, Singapore.
[3]	Pérez-Gavilán, J.J. and Aliabadi, M.H. (2000). A Galerkin Boundary Element Formulation with Dual Reciprocity for Elastodynamics. International Journal of Numerical Methods in Engineering, 48, 1331-1344.
[4]	Blobner, J., Hribersek, M. and Kuhn, G. (2000). Dual Reciprocity BEM-BDIM Technique for Conjugate Heat Transfer Computations. Computer Methods in Applied Mechanics and Engineering, 190, 1105-1116.
[5]	Skerget, L., Alujevic, A., Brebbia, C.A. and Kuhn, G. (1989). Natural and Forced Convection Simulation using the Velocity-vorticity Approach. Topics in Boundary Element Research, 5, 49-86.
[6]	Skerget, L., Hribersek, M. and Kuhn, G. (1999). Computational Fluid Dynamics by Boundary-Domain Integral Method. International Journal of Numerical Methods in Engineering, 46, 1291-1311.
[7]	Jecl, R., Skerget, L. and Petresin, E. (2001). Boundary Domain Integral Method for Transport Phenomena in Porous Media. International Journal of Numerical Methods in Engineering, 35, 39-54.
[8]	Jecl, R. and Skerget, L. (2003). Boundary element method for natural convection in non-Newtonian fluid saturated square porous cavity. Engineering Analysis in Boundary Elements, 23, 963-975.
ACTA GeOTeCHNICA SLOVENICA, 2009/l 17.
R. 1ECL ET AL.: DOUBLE DIFFUSE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD
[9]	Wu, J.C. (1982). Problem of General Viscous Flow. Chapter in Developments in BEM, Vol. 2., P.K. Banerjee, R.P. Shaw (eds.), Elsevier Applied Science Publication, London, UK.
[10]	Hribersek, M. and Skerget, L. (1996). Iterative Methods in Solving Navier-Stokes equations by the Boundary Element Method. International Journal of Numerical Methods in Engineering. 39, 115-139.
[11]	Bear, J. and Bachmat, Y. (1991). Introduction to Modeling of Transport Phenomena in Porous Media. Kluwer Academic Publishers, Dordrecht, Boston, London.
[12]	Skerget, L., and Jecl, R. (2002). Boundary element method for transport phenomena in porous medium. Chapter in Transport phenomena in porous media II, Derek B. Ingham, Ioan Pop (eds.), Elsevier Science.
[13]	Skerget, L. and Samec, N. (2005). BEM for the two-dimensional plane compressible fluid dynamics. Engineering Analysis in Boundary Elements, 29, 41-57.
[14]	Brebbia, C.A. and Dominquez, J. (1992). Boundary elements. An Introductory Course. McGraw-Hill, New York.
[15]	Amahmid, A., Hasnaoui, M., Mamou, M. and Vasseur, P. (1999). Double-diffusive parallel flow induced in a horizontal Brinkman porous layer subjected to constant heat and mass fluxes: analytical and numerical studies. Heat and Mass Transfer, 35, 409-421.
[16]	Kladias, N. and Prasad, V. (1989). Natural convection in horizontal porous layers: Effects of Prandtl and Darcy numbers. Journal of Heat Transfer, 111, 926-925.
[17]	Nield, D.A. (1968). Onset of thermohaline convection in a porous medium. Water Resources Research, 4, 553-560.
[18]	Nield, D.A., Manole, D.M. and Lage, J.L. (1993). Convetion induced by inclined thermal and solutal gradients in a shallow horizontal layer of a porous medium. Journal of Fluid mechanics, 257, 559-574.
[19]	Rudraih, N., Srimani, P.K. and Friedrich, R. (1982). Finite amplitude convection in a two component fluid saturated porous layer. International Journal of Heat and Mass Transfer, 25, 715-722.
[20]	Trevisan, O.V. and Bejan, A. (1987). Mass and heat transfer by high Rayleigh number convection in a porous medium heated from below. International Journal of Heat and Mass Transfer, 30, 2341-2356.
[21]	Rosenberg, N.D. and Spera, F. J. (1992). Thermohaline convection in a porous medium heated from below. International Journal of Heat and Mass Transfer, 35, 1261-1273.
[22]	Mamou, M., Vasseur, P., Bilgen, E. and Gobin, D. (1995). Double-diffusive convection in an inclined
slot filled with porous medium. European Journal of Mechanics B/Fluids, 14, 629-652.
[23]	Mohamad, A. A. and Bennacer, R. (2001). Natural convection in a confined saturated porous medium with horizontal temperature and vertical solutal gradients. International Journal of Thermal Sciences, 40, 2001.
[24]	Kalla, L., Vasseur, P., Beji, H. and Duval, R. (2001). Double diffusive convection within a horizontal porous layer salted from the bottom and heated horizontally. International Communications in Heat and Mass Transfer, 28, 1-10.
[25]	Mohamad, A. A. and Bennacer, R. (2002). Double diffusion natural convection in an enclosure filled with saturated porous medium subjected to cross gradients; stably stratified fluid. International Journal of Heat and Mass Transfer, 45, 3725-3740.
[26]	Ravnik, J., Skerget, L. and Hribersek M. (2004). The wavelet transform for BEM computational fluid dynamics. Engineering Analysis in Boundary Elements, 28, 1303-1314.
18. ACTA GeOTeCHNICA SLOVENICA, 2009/1