© Acta hydrotechnica 23/38 (2005), Ljubljana
ISSN 1581-0267
1
UDK / UDC: 519.61/.64:532.5 Prejeto / Received: 14. 7. 2005
Izvirni znanstveni članek – Original scientific paper Sprejeto / Accepted: 9. 1. 2006
PRIMERJA V A MED FORCHEIMERJEVIM IN BRINKMANOVIM
MODELOM KONVEKTIVNEGA TOKA V POROZNI KOTANJI Z
ROBNOOBMO ČNO INTEGRALSKO METODO
COMPARISON BETWEEN THE FORCHEIMER AND THE BRINKMAN
MODEL FOR CONVECTIVE FLOW IN POROUS CA VITY WITH
BOUNDARY DOMAIN INTEGRAL METHOD
Renata JECL, Leopold ŠKERGET, Janja KRAMER
V prispevku je predstavljena uporaba Forcheimerjevega modela za analizo konvektivnega toka v
porozni kotanji s poudarkom na dolo čitvi vpliva dodatnega (Forcheimerjevega) vztrajnostnega
člena v gibalni ena čbi na hitrost in skupni prenos toplote. Znano je, da je vpliv tega člena pri
inženirskih problemih v zasi čeni porozni snovi minimalen, kar pa še ni bilo dokazano z uporabo
metode robnih elementov (MRE) ali s katero od njenih razširitev. V prispevku je zato podana
numeri čna rešitev problema naravne konvekcije v porozni kotanji z uporabo robnoobmo čne
integralske metode (ROIM), pri čemer so hitrosti in skupni prenos toplote izra čunani za razli čna
modificirana Rayleighova in Darcyjeva števila, rešitve pa primerjane z objavljenimi rezultati na
osnovi Brinkmanovega modela.
Klju čne besede: Forcheimerjev model, robnoobmo čna integralska metoda, porozna snov, naravna
konvekcija
The use of the Forcheimer model for analyzing the convective flow in a porous cavity is
represented, emphasizing the investigation of the effect of the additional (Forcheimer) inertia term
in momentum equation on the velocity and global heat transfer through the cavity. The effects that
this term has upon the engineering parameters of interest for the case of a fluid saturated media,
are known to be minimal, however, this has not been confirmed yet with the use of the boundary
element method (BEM) or any of its extensions. Therefore the numerical solution of the problem of
natural convection in porous cavity using the boundary domain integral method (BDIM) is
presented, where the velocities and the total heat transfer are calculated for different modified
Rayleigh and Darcy numbers and the solutions are compared with published results obtained with
the Brinkman model.
Key words: Forcheimer model, boundary domain integral method, porous medium, natural
convection
1. UVOD
Z vzgonom povzro čen konvektivni tok v
zasi čeni porozni snovi je pomemben zaradi
mnogih aplikacij v inženirski praksi. V
literaturi najdemo kar nekaj študij, ki
obravnavajo ta problem v razli čnih
geometrijah in z razli čnimi toplotnimi robnimi
pogoji. Rezultati za kotanjo, pri kateri so
vertikalne stene razli čnih temperatur, dobljeni
z uporabo Brinkman-Darcyjeve formulacije,
so podani drugje (med drugim tudi v Lauriat &
1. INTRODUCTION
Buoyancy induced convection in a fluid
saturated porous media is of considerable
interest, due to its many applications in
energy-related engineering problems. Studies
have been reported dealing with different
geometries and a variety of heating conditions.
For example, a vertical cavity in which a
horizontal temperature gradient is induced by
side walls maintained at different temperatures
has been analyzed by using the Brinkman-

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
2
Prasad, 1987; Vasseur et al., 1990; Jecl et al.,
2001; Lauriat & Prasad, 1989), za primer
naravne konvekcije v horizontalnem poroznem
sloju, gretem od spodaj, pa v (Lage, 1992).
Numeri čne metode, uporabljene za reševanje
osnovnih ena čb, so v ve čini primerov metoda
končnih razlik (MKR) in metoda kon čnih
elementov (MKE), medtem ko je v pri čujo čem
prispevku uporabljena robnoobmo čna
integralska metoda (ROIM). Gibalna ena čba
toka teko čine skozi porozno snov je ekvivalent
Navier-Stokesovi gibalni ena čbi za čisto
teko čino. Pri študiju prenosnih pojavov v
porozni snovi najdemo v literaturi kar nekaj
razli čnih oblik gibalne ena čbe (Nield & Bejan,
1999). Namen našega prispevka je ugotoviti
vpliv Forcheimerjevega člena v gibalni ena čbi
na hitrost in skupni prenos toplote z uporabo
robnoobmo čne integralske metode.
Teorija laminarnega toka skozi porozno
snov temelji na eksperimentih, ki jih je pri
študiju enodimenzijskega toka vode skozi
homogen in nedeformabilen porozen material
izvajal francoski inženir Henry Darcy leta
1856. Rezultat njegovih poizkusov je linearni
zakon prepustnosti (Darcyjev zakon), ki
podaja proporcionalno zvezo med hitrostjo
toka in tla čno razliko kot () P v ∇ − =
 µ K ,
kjer je µ dinami čna viskoznost teko čine,
koeficient K pa prepustnost. Prepustnost je
neodvisna od teko čine, njena vrednost je
odvisna le od geometrije porozne snovi, v
splošnem je tenzor drugega reda. V primeru
izotropne porozne snovi je prepustnost skalar
in Darcyjev zakon lahko zapišemo kot:
extended Darcy formulation, given elsewhere
(among others in Lauriat & Prasad, 1987;
Vasseur et al., 1990; Jecl et al., 2001; Lauriat
& Prasad, 1989) and for the natural convection
in the porous layer heated from below by
(Lage, 1992). The numerical methods used for
the solution of governing equations are the
finite difference method (FDM) and the finite
element method (FEM) while in the present
contribution the boundary domain integral
method (BDIM) is used. The momentum
equation for flow of the fluid through the
porous media is similar to the Navier-Stokes
momentum equation for pure fluid. When
studying transport phenomena in porous media
different models of momentum equations
could be found (Nield & Bejan, 1999). The
aim of this paper is to determine the influence
of the Forcheimer term in momentum equation
to velocity and heat transfer using the
boundary domain integral method.
The theory of laminar flow through the
porous media is based on an experiment, of
unidirectional flow through a homogenous,
fixed and nondeformable porous medium,
originally performed by the French engineer
Henry Darcy (1856). The result of his
experiment was the linear law of permeability
(Darcy’s law), where the relationship between
velocity and pressure difference is given as
( ) P v ∇ − =
 µ K , where µ is the fluid
dynamic viscosity and K is permeability of
the porous media. Permeability depends only
on geometry of porous media and is generally
a second order tensor. In case of an isotropic
porous media the permeability is a scalar and
the Darcy law can be written as:
v
K
P
 µ
− = ∇ . (1)
Tako kot pri katerem koli problemu,
definiranem v toku teko čine, je tudi v porozni
snovi režim toka definiran z Reynoldsovim
številom Re, ki dolo ča razmerje med
vztrajnostnimi in viskoznimi silami. V
splošnem je obveljalo prepri čanje, da Darcyjev
zakon velja, dokler je 10 Re < – hitrosti so
majhne in govorimo o laminarnem režimu toka
skozi porozno snov, pri katerem prevladujejo
As in any fluid flow problem, the flow
regime in porous media is given with
Reynolds number Re, defined as the ratio
between inertia and viscous forces. It is
generally believed that the Darcy law is valid
until 10 Re < , velocities are small, the flow
regime through porous media is laminar and
viscous forces are predominant. When the
Reynolds number is greater, the inertia forces

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
3
viskozne sile. Ko pa je Reynoldsovo število
ve čje, so ve čje tudi vztrajnostne sile in
Darcyjevemu zakonu je treba dodati člene, s
katerimi zajamemo vpliv pove čanja hitrosti.
Najenostavnejša dopolnitev linearnega
zakona je Darcyjev zakon z upoštevanjem
vztrajnostnih u činkov po analogiji z Navier-
Stokesovo ena čbo:
are higher than the viscous ones and some new
terms have to be added to the Darcy equation
to take into account the effect of the velocity
deviations.
First, a very simple extension of the Darcy
law is the one with the addition of the terms
covering the influence of inertia effects
analogous to the Navier-Stokes equation:
() v
K
P v v
1
t
v 1
2
 µ
φ ∂
∂
φ
ρ − ∇ − =
⎥
⎦
⎤
⎢
⎣
⎡
∇ + , (2)
kjer so ρ gostota teko čine, φ poroznost in t
čas.
Naslednja velikokrat uporabljena gibalna
ena čba je Brinkmanova ena čba:
where ρ is the fluid density, φ is the porosity
and t time.
Next, often used momentum equation is the
Brinkman momentum equation:
() v v
K
P v v
1
t
v 1
2
e 2
 ∇ + − ∇ − =
⎥
⎦
⎤
⎢
⎣
⎡
∇ + µ
µ
φ ∂
∂
φ
ρ ,( 3 )
kjer 
e
µ predstavlja dinami čno viskoznost
porozne snovi. Z ena čbo je mogo če zadostiti
brezzdrsnemu robnemu pogoju, ki pravi, da je
hitrost na trdnem robu (na površini, ki omejuje
porozno snov) enaka ni č. Pri tem velja
omeniti, da se Brinkmanova ena čba prevede v
Navier-Stokesovo gibalno ena čbo za čisto
nestisljivo teko čino, ko prepustnost naraš ča
preko vseh meja ( ∞ → K ), ko pa se
prepustnost zmanjšuje in približuje ni č
( 0 K → ), dobimo Darcyjevo ena čbo.
V literaturi zasledimo tudi uporabo
Forcheimerjeve gibalne ena čbe:
where 
e
µ is the effective dynamic viscosity.
The equation enables us to satisfy the no-slip
boundary condition on impermeable surfaces
that bound the porous media (the velocity on
the boundaries is set to zero). It is important to
stress out that the Brinkman equation reduces
to the classical Navier Stokes equation for
pure uncompressible fluid when the
permeability tends to infinity ( ∞ → K ). When
the permeability approaches zero ( 0 K → ) we
recover the Darcy equation.
In some references the Forcheimer
momentum equation is also introduced as:
() v v
K
c
v
K
P v v
1
t
v 1
2
1
F
2
 ρ µ
φ ∂
∂
φ
ρ − − ∇ − =
⎥
⎦
⎤
⎢
⎣
⎡
∇ + ,( 4 )
kjer je 
F
c koeficient odvisen od geometrije
por in prepustnosti in je definiran v
nadaljevanju v ena čbi (8).
where 
F
c is the coefficient, which varies with
the nature of porous media and permeability
and is defined below in equation (8).
2. OSNOVNE ENA ČBE 2. GOVERNING EQUATIONS
Makroskopske ena čbe ohranitve mase,
gibalne koli čine in energije v porozni snovi
dobimo s postopkom povpre čenja po
reprezentativnem elementarnem volumnu
The general set of macroscopic equations
for conservation of mass, momentum and
energy in porous media are obtained by
averaging the Navier-Stokes equation for pure

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
4
(Bear & Bachmat, 1991) iz Navier-Stokesovih
ena čb za čisto nestisljivo teko čino ob
upoštevanju dejstva, da je za tok teko čine na
voljo samo del obravnavanega volumna
(izražen s poroznostjo). Te ena čbe v literaturi
najdemo pod imenom Brinkmanov model, ki
ga sestavljajo kontinuitetna, Brinkmanova
gibalna in energijska ena čba (Lauriat &
Prasad, 1987). V našem primeru smo gibalno
ena čbo dopolnili še s Forcheimerjevim
členom, dobljeni sistem pa je v literaturi znan
kot Forcheimerjev model (Lauriat & Prasad,
1989) v obliki:
uncompressible fluid over a representative
elementary volume, (Bear & Bachmat, 1991),
considering that only a part of the volume
(expressed with porosity) is available for the
flow of the fluid. These equations are known
as the Brinkman model and consist of
continuity equation, Brinkman momentum
equation and energy equation, (Lauriat &
Prasad, 1987). In our case the Forcheimer term
is added to the momentum equation, and the
whole set, in the literature known as the
Forcheimer model (Lauriat & Prasad, 1987),
can be written as:
0
x
v
i
i
=
∂
∂
,( 5 )
i
F
j j
i
2
e i
f
i
i j
i j
2
i
v v
K
c
x x
v
v
K
g F
x
P 1
x
v v
1
t
v 1
−
∂ ∂
∂
+ − +
∂
∂
− =
∂
∂
+
∂
∂
ν
ν
ρ φ φ
,( 6 )
()() () [] ()
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
=
∂
∂
+ − +
∂
∂
j
e
j j
j
f s s f
x
T
x x
T v
c T c 1 c
t
λ ρ ρ φ ρ φ ,( 7 )
kjer je zadnji člen na desni strani ena čbe (6)
dodaten Forcheimerjev vztrajnostni člen zaradi
nelinearnih u činkov, ki so posledica pove čanja
hitrosti. Ostale oznake so: 
i
v i-ta komponenta
filtracijske hitrosti, v absolutna vrednost
hitrosti (
2
z
2
y
2
x
v v v v v + + = = ), 
i
x i-ta
koordinata, φ poroznost, 
f
ν kinemati čna
viskoznost teko čine, 
e
ν kinemati čna
viskoznost porozne snovi, K prepustnost
porozne snovi, 
i
x P ∂ ∂ tla čni gradient, ρ
gostota teko čine, 
i
g gravitacijski pospešek.
Povezava med gostoto in temperaturo je
podana s funkcijo F kot
() ()
0 T 0 0
T T F − − = − = β ρ ρ ρ , kjer je 
0
ρ
referen čna gostota teko čine pri temperaturi 
0
T
in 
T
β koeficient temperaturnega volumskega
raztezka. Nadalje so 
s
ρ in ρ gostote trdne in
teko če faze porozne snovi, 
s
c in 
f
c specifi čne
toplote pri konstantnem tlaku za trdno in
teko čo fazo, T je temperatura, 
e
λ pa
predstavlja dejansko toplotno prevodnost
zasi čene porozne snovi. Ta se obi čajno dolo či
z uporabo klasi čnega ’’mešalnega’’ pravila
where the last term on the r.h.s. of equation (6)
is the additional Forcheimer inertia term,
which includes non-linear effects due to the
increase of velocity. Other parameters are 
i
v
volume-averaged velocity, v absolute value of
velocity (
2
z
2
y
2
x
v v v v v + + = = ), 
i
x the i -th
coordinate, φ is porosity, 
f
ν the fluid
kinematic viscosity, 
e
ν the effective kinematic
viscosity, K permeability of porous media,
i
x P ∂ ∂ pressure gradient in the flow
direction, ρ the fluid density, 
i
g gravity. The
normalized density-temperature variation
function F is ( ) ( )
0 T 0 0
T T F − − = − = β ρ ρ ρ ,
with 
0
ρ denoting the reference fluid mass
density at temperature 
0
T and 
T
β being the
thermal volume expansion coefficient of the
fluid. Furthermore 
s
ρ and ρ are the solid and
fluid densities respectively, 
s
c and 
f
c the
solid and the fluid specific heats at constant
pressure, T stands for temperature, and 
e
λ
represents the effective thermal conductivity
of the saturated porous media, which is given
as ( )
s f e
1 λ φ λ φ λ − + = , where 
f
λ and 
s
λ are

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
5
()
s f e
1 λ φ λ φ λ − + = , kjer 
f
λ in 
s
λ
predstavljata toplotne prevodnosti teko čine in
trdne snovi. Forcheimerjev člen – zadnji člen
na desni strani ena čbe (6) povezuje kvadrat
hitrosti s koeficientom 
F
c , odvisnim od
geometrije por s prepustnostjo K. Oba
koeficienta lahko izra čunamo po Ergunu,
podano v (Nield and Bejan, 1999), v
odvisnosti od poroznosti z naslednjimi
relacijami:
thermal conductivities for fluid and solid
matter, respectively. Forcheimer term – the
last on the r.h.s. of equation (6) connects the
square of velocity with coefficient 
F
c ,
dimensionless form-drag constant varying with
nature of porous media, and permeability K .
Using the Ergun model, as described in (Nield
and Bejan, 1999), it is possible to calculate
both coefficients as:
()
2 1
3
F 2
3 2
150
75 . 1
c ;
) 1 ( 150
d
K
φ
φ
φ
=
−
= ,( 8 )
kjer je d srednja vrednost premera delca.
Gibalna ena čba (6), imenovana tudi
Brinkman-Forcheimerjeva ena čba, ima v
našem primeru dva viskozna in dva
vztrajnostna člena. Prvi viskozni je obi čajen
Darcyjev člen (tretji na desni), drugi pa je
Brinkmanov člen ( četrti na desni), ki je
analogen Laplaceovemu členu v Navier-
Stokesovi ena čbi za čisto teko čino. V tem
členu (Brinkmanov člen ali Brinkmanova
razširitev) nastopa dejanska viskoznost 
e
ν , ki
je odvisna od geometrije porozne snovi in ima
praviloma druga čno vrednost od viskoznosti
teko čine 
f
ν . Prvi vztrajnostni člen je
konvektivni vztrajnostni člen, podan kot
produkt hitrosti in divergence hitrosti (drugi
člen na levi), drugi pa je Forcheimerjev
vztrajnostni člen, zapisan kot produkt hitrosti
in njene absolutne vrednosti (zadnji na desni).
where d is the mean particle diameter.
The momentum equation (6), commonly
known as the Brinkman-Forcheimer equation,
consists of two viscous and two inertia terms.
The first viscous term is the usual Darcy term
(third on the r.h.s.), and the second is
analogous to the Laplacian term that appears
in the Navier-Stokes equations for pure fluid
(fourth on the r.h.s.). In the Laplace term
(Brinkman term or Brinkman extension) the
effective viscosity 
e
ν is used which normally
has a different value than the fluid viscosity
f
ν . The first inertia term is the convective
inertia term represented by the velocity times
its divergent (second on the l.h.s.) and the
second inertia term is the Forcheimer inertia
term (drag) represented by the velocity times
its absolute value (last term on the r.h.s.).
3. ROBNOOBMO ČNA
INTEGRALSKA METODA
3. BOUNDARY DOMAIN
INTEGRAL METHOD
Numeri čna metoda, izbrana za reševanje
problema, je robnoobmo čna integralska
metoda (ROIM), katere osnova je klasi čna
metoda robnih elementov (MRE, Škerget et
al., 1999). Če viskoznost razdelimo na
konstantni in spremenljivi del, tako da velja
f f f
ν ν ν
~
+ = , se Brinkmanov člen v gibalni
ena čbi razdeli na dva dela in ena čba (6) se
zapiše kot:
The numerical method chosen for this
investigation is the Boundary Domain Integral
Method (BDIM) based on the classical
Boundary Element Method (BEM, Škerget et
al., 1999). If the viscosity is partitioned into
constant and variable parts so that
f f f
ν ν ν
~
+ = , the Brinkman extension in
momentum equation is divided into two parts
and the equation (6) becomes:

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
6
()
i
F
ij f
j
j j
i
2
f i
f
i
i j
i j
2
i
v v
K
c
~
2
x
x x
v
v
K
g F
x
P 1
x
v v
1
t
v 1
−
∂
∂
+
∂ ∂
∂
+ − +
∂
∂
− =
∂
∂
+
∂
∂
ε ν Λ
ν Λ
ν
ρ φ φ
 ,( 9 )
kjer je 
ij
ε deformacijski tenzor,
()
i j j i ij
x v x v 2 1 ∂ ∂ + ∂ ∂ = ε . Parameter Λ
predstavlja razmerje med dejansko
viskoznostjo porozne snovi 
e
ν in viskoznostjo
teko čine 
f
ν , definirano kot 
f e
ν ν Λ = .
Nadalje se toplotna difuzivnost porozne snovi
p
a , definirana kot 
f e p
c a ρ λ = , podobno kot
kinemati čna viskoznost razdeli na konstantni
in spremenljivi del 
p p p
a
~
a a + = . Če razmerje
toplotnih kapacitet zapišemo z izrazom
() () ()
f s s f
c c 1 c ρ ρ φ ρ φ σ − + = , se ena čba
ohranitve energije (7) preoblikuje v naslednjo
obliko:
where 
ij
ε represent the strain rate tensor,
( )
i j j i ij
x v x v 2 1 ∂ ∂ + ∂ ∂ = ε . Parameter Λ
represents the ratio between the effective
viscosity of the porous media 
e
ν and the fluid
viscosity 
f
ν defined as 
f e
ν ν Λ = .
Furthermore, the thermal diffusivity of the
porous media 
p
a , defined as 
f e p
c a ρ λ = , is
similarly as the kinematic viscosity,
partitioned into constant and variable parts
p p p
a
~
a a + = . Introducing the heat capacity
ratio with ( )() ()
f s s f
c c 1 c ρ ρ φ ρ φ σ − + = ,
the heat energy equation (7) can be rewritten
in the following form:
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
+
∂ ∂
∂
=
∂
∂
+
∂
∂
j
p
j j j
2
p
j
j
x
T
a
~
x x x
T
a
x
T v
t
T
σ . (10)
Vodilne ohranitvene ena čbe (5), (9) in (10)
z uporabo hitrostno vrtin čne formulacije
zapišemo v obliki prenosnih ena čb za
kinematiko, kinetiko vrtin čnosti in
temperaturno kinetiko (Jecl et al., 2001). Z
vektorjem vtrin čnosti 
i
ω , ki predstavlja rotor
hitrostnega polja
The governing equations (5), (9) and (10)
are in the next step transformed with the use of
the velocity-vorticity variables formulation
into transport equations for kinematics and
kinetics (Jecl et al., 2001). With the vorticity
vector 
i
ω , representing the curl of the velocity
field
0
x
,
x
v
e
j
j
j
k
ijk i
=
∂
∂
∂
∂
=
ω
ω (11)
in je po definiciji solenoidni vektor, se
ra čunska shema razdeli na kinemati čni in
kineti čni del, tako da se kontinuitetna ena čba
in ena čba ohranitve gibalne koli čine ter
energije nadomestita z ena čbama kinematike
in kinetike.
Kinematika je podana v obliki vektorske
elipti čne Poissonove ena čbe za hitrostni vektor
which is a solenoidal vector function by its
definition the computational scheme is
partitioned into its kinematic and kinetic parts
so that the continuity and momentum
equations are replaced by the equations of
kinematics and kinetics.
The kinematics is given by vector elliptic
Poisson equation for the velocity  vector
0
x
e
x x
v
j
k
ijk
j j
i
2
=
∂
∂
+
∂ ∂
∂ ω
, (12)

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
7
ki izraža kompatibilitetni in omejitveni pogoj
med vektorskim poljem hitrosti in poljem
vektorja vrtin čnosti. Z namenom pospešitve
konvergence in pove čanja stabilnosti vezane
hitrostno-vrtin čne iterativne sheme, ena čbo
(12) zapišemo v nepravi nestacionarni obliki,
tako da dodamo umetni akumulacijski člen in
jo zapišemo v naslednji paraboli čni difuzijski
obliki:
representing the compatibility and restriction
conditions between velocity and vorticity field
functions. To accelerate the convergence and
the stability of the coupled velocity-vorticity
iterative scheme, the false transient approach
is used by adding the artificial accumulation
term and hence the equation (12) can be
written as a parabolic diffusion one:
0
x
e
t
v 1
x x
v
j
k
ijk
i
j j
i
2
=
∂
∂
+
∂
∂
−
∂ ∂
∂ ω
α
, (13)
kjer je α relaksacijski parameter. O čitno je,
da je vodilna hitrostna ena čba (12) natan čno
izpolnjena v stacionarnem stanju ( ∞ → t ), ko
umetni časovni odvod hitrosti odpade.
Vrtin čno kinetiko podamo z nelinearno
paraboli čno difuzivno-konvektivno prenosno
ena čbo vrtin čnosti, ki jo izpeljemo z
delovanjem operatorja rotor na gibalno ena čbo
(6), pri čemer vpeljemo še novo spremenljivko
v
τ , ki predstavlja modificirani vrtin čni časovni
korak, definiran kot φ τ t
v
= :
where α is a relaxation parameter. It is
obvious that the governing velocity equation
(12) is exactly satisfied at the steady state
( ∞ → t ), when the false time derivative
vanishes.
The vorticity transport is governed by the
nonlinear parabolic diffusion-convection
equation obtained as a curl of the momentum
equation (6), where the new variable 
v
τ , the so
called modified vorticity time step, is
introduced as φ τ t
v
= :
j
k ijk
F 2
i
F 2
j
ij 2
j
i
f
j
2
i
f 2
j
k ijk
2
j j
i
2
f
2
j
i j
j
i
j
v
x
v
v e
K
c
v
K
c
x
f
x
~
x K x
F
g e
x x x
v
x
v
∂
∂
− −
∂
∂
+
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
+ −
∂
∂
+
∂ ∂
∂
+
∂
∂
=
∂
∂
+
∂
∂
φ ω φ Λ φ
ω
ν Λ φ ω
ν
φ φ
ω
ν Λ φ
ω
ω
τ
ω
, (14)
V ena čbi (14) člen 
ij
f predstavlja prispevek
zaradi u činkov spremenljivih lastnosti snovi.
Forcheimerjev člen je razdeljen na dva dela.
Prvi del (sedmi člen desne strani ena čbe) bo
skupaj z Darcyjevim členom pridružen
sistemski matriki, drugi del (zadnji člen desne
strani ena čbe), ki vklju čuje prvi odvod hitrosti,
pa bo skupaj z vzgonskim členom in členom
nelinearnih lastnosti snovi modeliran kot del
izvornega člena. Zaradi uporabe hitrostno-
vrtin čne formulacije pri transformaciji
energijske ena čbe vpeljemo še modificiran
temperaturni časovni korak 
T
τ , σ τ t
T
= v
energijsko ena čbo (10) in jo zapišemo v
naslednji obliki:
In equation (14), the term 
ij
f represents a
contribution arising on account of non-linear
material properties. The Forcheimer term is
split into two parts. The first part (seventh
term on the r.h.s) will be combined with the
Darcy term and numerically treated in the
system matrix. The second part (last term on
the r.h.s), involving first-order derivatives of
the velocity modulus as well as the buoyancy
term and non-linear material properties term,
will be modelled as a source term. For the
same reason as with vorticity kinetics, we
introduce the new modified temperature time
step σ τ t
T
= into the energy equation which
permits us to rewrite equation (10) in the form:

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
8
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
∂
∂
+
∂ ∂
∂
=
∂
∂
+
∂
∂
j
p
j j j
2
p
j
j
T
x
T
a
~
x x x
T
a
x
T
v
T
τ
. (15)
Ena čbe (13), (14) in (15) predstavljajo
nelinearni sistem ena čb, ki jih z uporabo
metode utežnih ostankov (v kombinaciji z
ustreznimi Greenovimi osnovnimi rešitvami)
preoblikujemo v ustrezen sistem integralskih
ena čb. Vsaka od teh ena čb se lahko zapiše v
obliki naslednje splošne diferencialne
ohranitvene ena čbe:
The equations (13), (14) and (15) represent
the non-linear set of equations to which the
weighted residuals technique (in combination
with the suitable Green function) has to be
applied resulting in a set of integral equations.
Each of these equations can be written as the
following general differential conservation
equation:
[ ] 0 b u = + ℵ , (16)
kjer je linearni diferencialni operator [ ] ⋅ ℵ
paraboli čni difuzijski diferencialni operator,
zapisan v obliki:
where the linear differential operator [ ] ⋅ ℵ
will be the parabolic diffusion differential
operator of the form:
[]
t
u
x x
u
a u
j j
2
∂
∂
−
∂ ∂
∂
= ℵ , (17)
u je poljubna funkcija polja (hitrost,
vrtin čnost ali temperatura), člen b pa
predstavlja nehomogeni del nelinearnih
vplivov oziroma volumskih izvorov.
Spremenljivka a ustreza α za kinematiko,
f
2
ν Λ φ za vrtin čno in 
p
a za temperaturno
kinetiko. Ker so numeri čni rezultati v
pri čujo čem prispevku omejeni na
dvodimenzionalni problem, bodo vse ena čbe v
nadaljevanju zapisane samo za primer
ravninske geometrije. Integralska predstavitev
kinematike se preoblikuje in zapiše v naslednji
obliki:
u is an arbitrary field function (velocity,
vorticity or temperature), while the
nonhomogeneous term b is generally used for
the nonlinear transport effects or pseudo body
forces. The term a equals α for kinematics,
f
2
ν Λ φ and 
p
a for vorticity and temperature
kinetics, respectively. As the computational
results of the present work are limited to the
two-dimensional case, all the subsequent
integral equations will consequently be written
for the case of planar geometry only. The
integral representation of the kinematic
equation is given with the following equation:
()( )
Ω Ω ω α
Γ ω α α Γ α ξ ξ
Ω Ω
Γ Γ
d u v d t d
x
u
e
d t d u n e
n
v
d t d
n
u
v t , v c
*
1 F 1 F , i
t
t j
*
ij
t
t
*
j ij
i
t
t
*
i F i
F
1 F
F
1 F
F
1 F
∫ ∫∫
∫∫ ∫∫
− −
+
∂
∂
−
⎟
⎠
⎞
⎜
⎝
⎛
+
∂
∂
=
∂
∂
+
−
− −
. (18)

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
9
Parameter ) ( c ξ predstavlja geometrijski
koeficient, povezan z osnovno rešitvijo, in je
odvisen od lege izvorne to čke ξ , Γ pa je rob
obmo čja Ω . Ena čba (18) opisuje kinematiko
toka teko čine skozi porozno snov v integralski
obliki.
Transportna ena čba vrtin čnosti v integralski
obliki je podana z naslednjim izrazom:
Parameter ) ( c ξ denotes the fundamental
solution related coefficient depending on the
position of the source point ξ , Γ is the
boundary of domain Ω. Equation (18)
expresses the kinematic of fluid flow in the
porous medium in the integral form.
The vorticity transport equation is given
with the following integral representation:
∫ ∫∫
∫∫
∫∫ ∫∫
− −
+ ⎟
⎠
⎞
⎜
⎝
⎛
+ −
∂
∂
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
+ − −
∂
∂
− +
⎟
⎟
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎜
⎜
⎝
⎛
−
+ + −
∂
∂
=
∂
∂
+
−
−
− −
Ω Ω
τ
τ
Ω
τ
τ
Γ
τ
τ Γ
τ
τ
Ω ω
τ ∆
Ω τ ω φ
ν
φ
Ω τ φ Λ φ φ
ω
ν Λ φ ω
Γ τ
φ Λ φ
φ ω
ω
ν Λ φ
Γ τ ω ν Λ φ ξ ω ξ
d u
1
d d u v
K
c
K
d d
x
u
v v e
K
c
f F g e
x
~
v
d d u
n v v e
K
c
n f
n F g e v
n
d d
n
u
) ( ) ( c
*
1 F 1 F
v
v
* F 2 f 2
v
j
*
j ij
F 2
j
2
j ij
2
j
f
2
j
v
*
j j ij
F 2
j j
2
j j ij
2
n f
2
v
*
f
2
F
1 F
F
1 F
F
1 F
F
1 F
. (19)
Ob uporabi opisanega postopka
preoblikovanja parcialne diferencialne ena čbe
ohranitve energije (15) izpeljemo naslednjo
integralsko obliko:
Applying a similar procedure to the heat
transport equation (15) finally the next integral
statement is derived:
∫ ∫∫
∫∫ ∫∫
− −
+
∂
∂
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
∂
∂
− +
⎟
⎠
⎞
⎜
⎝
⎛
−
∂
∂
=
∂
∂
+
−
− −
Ω Ω
τ
τ
Γ
τ
τ Γ
τ
τ
Ω
τ ∆
Ω τ
Γ τ Γ τ ξ ξ
d u T
1
d d
x
u
x
T
a
~
T v
~
d d u v T
n
T
a d d
n
u
T a ) ( T ) ( c
*
1 F 1 F
T
T
j
*
j
p j
T
*
n p T
*
p
F
1 F
F
1 F
F
1 F
. (20)
V vseh integralskih ena čbah smo upoštevali
konstantni potek vseh funkcij polja () T , , v
i
ω
znotraj posameznega časovnega intervala
1 F F
t
−
− = τ τ ∆ , kar omogo ča, da lahko
časovne integrale rešimo analiti čno, pri čemer
velja:
In all integral equations the constant
variation of field functions () T , , v
i
ω is
assumed within the individual time increment
1 F F
t
−
− = τ τ ∆ , therefore the time integrals
may be evaluated analytically, for example:
∫
−
=
F
1 F
dt u a U
* *
τ
τ
. (21)
Prav tako smo pri preoblikovanju vseh
diferencialnih ena čb v integralske uporabili
paraboli čno difuzijsko osnovno rešitev 
*
u v
obliki:
Also in the transformation of all differential
equations into integral ones, the parabolic
diffusion fundamental solution 
*
u is used as
given by:

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
10
at 4
r
*
2
e
t a 4
1
u
π
= , (22)
kjer je pravilna uporaba spremenljivke a
opisana zgoraj in je r dolžina vektorja od
izvorne do referen čne to čke.
Za numeri čno aproksimativno rešitev
razli čnih obravnavanih funkcij polja, npr.
hitrosti, vrtin čnosti, temperature, moramo
pripadajo če robnoobmo čne integralske ena čbe
zapisati v diskretni obliki, tako da robne in
obmo čne integrale aproksimiramo z vsoto
integralov po posameznih robnih elementih in
notranjih celicah. Rezultirajo ča matri čna
oblika ena čb za kinematiko, vrtin čno kinetiko
in temperaturno kinetiko se zapiše v naslednji
obliki:
where the proper use of a is as expressed
above and r is the magnitude of the vector
from the source point to the reference point.
For the numerical approximate solution of
the field functions, namely the velocity,
vorticity, and temperature, the integral
equations are further written in a discretized
manner in which the integrals over the
boundary and domain are approximated by a
sum of the integrals over all boundary
elements and over all internal cells. The
resulting matrix forms of the equations for
kinematics, vorticity kinetics and heat energy
kinetics are represented here as:
[] {} [ ] { } [ ] { }
[] {}[] {}[] {} ω
ω
x x y
y y x
v v
v v
D H H
D H H
t
t
− =
− =
, (23)
[] {} [] []
[] {} [] {}
[]
⎭
⎬
⎫
⎩
⎨
⎧
−
∂
∂
− − − +
+
⎥
⎦
⎤
⎢
⎣
⎡
+ −
⎭
⎬
⎫
⎩
⎨
⎧
+ + + −
⎭
⎬
⎫
⎩
⎨
⎧
=
−
v v e
K
c
x
~
f F g e v
1
1
v
K
c
K
1
v v e
K
c
n f n F g n e v
1
n
j ij
F 2
j
f
2
j
2
j ij
2
j
f
2
1 F
f
2
F 2 f
2
f
2
n ij
F 2
j j
2
j j i ij
2
n
f
2
f
f
φ
ω
ν Λ φ Λ φ φ ω
ν Λ φ
ω β
ν Λ φ
ω φ
ν φ
ν Λ φ
φ Λ φ φ ω
ν Λ φ ∂
ω ∂
ν
ν
ω
j
D
B B
G G H
, (24)
[] {} [] [] [] {} [] [] {}
[] [] {}
1 F
p j p
p
j
p
n
p p
p
T
a
1
x
T
a
a
~
T v
a
1
T v
a
1
n
T
a
a
T
−
+
⎭
⎬
⎫
⎩
⎨
⎧
− + −
⎭
⎬
⎫
⎩
⎨
⎧
=
B D
D G G H
j
j
β
∂
∂
∂
∂
. (25)
V zgornjih ena čbah pomenijo [][] [ ] D , G , H
in [] B matrike, sestavljene iz integralov, ki
predstavljajo integracijo po robnih elementih
in notranjih celicah, zapisanih za vsa robna in
obmo čna vozliš ča. Vidimo, da je
Forcheimerjev člen v ena čbi (24) razdeljen na
tri dele, pri čemer sta prispevka, pomnožena z
matrikama [] G in [] B , pridružena sistemski
matriki, zadnji del člena z matriko []
j
D pa je
In the equations given above the matrices
[ ] [ ] [ ] D , G , H and [ ] B are the influence
matrices and they are composed of integrals
taken over the individual boundary elements
and over the internal cells. The Forcheimer
term in vorticity kinetic matrix equation,
equation (24) is split into three parts. The first
two contributions multiplied with matrices [ ] G
and [ ] B are numerically treated in the system

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
11
obravnavan kot del izvornega člena oziroma
člena, s katerim zajamemo vpliv nelinearnosti.
Sistem diskretiziranih ena čb je rešen kot vezan
sistem kineti čnih in kinemati čnih ena čb,
upoštevaje nelinearne osnovne konstitutivne
zveze in ustrezne robne in za četne pogoje. Ker
se implicitni sistem ena čb zapiše isto časno za
vsa robna in obmo čna vozliš ča, dobimo po
opisanem postopku zelo veliko, polno
sistemsko matriko, ki vsebuje vplive difuzije
in konvekcije. Rezultat je stabilna in natan čna
numeri čna shema, ki pa za reševanje potrebuje
veliko ra čunalniškega spomina in časa. Za
izboljšanje oziroma optimizacijo izra čuna se
uporablja tehnika podobmo čij, pri čemer se
celotno obmo čje rešitve razdeli na
podobmo čja, pri katerih uporabimo enak
numeri čni postopek, kot je opisan zgoraj.
Kon čni sistem ena čb za celotno obmo čje
dobimo z združevanjem sistemov ena čb za
vsako posamezno podobmo čje, upoštevaje
ustrezne pogoje enakosti in kontinuitetne
pogoje na vmesnih robovih. Rezultirajo ča
sistemska matrika je veliko bolj prazna, tako
da je dobljeni matri čni sistem primernejši za
reševanje z iterativnimi tehnikami. V našem
primeru je vsako podobmo čje sestavljeno iz
štirih nezveznih 3-to čkovnih kvadratnih robnih
elementov in ene zvezne 9-to čkovne kvadratne
notranje oziroma obmo čne celice (Jecl &
Škerget, 2003).
matrix and the last part, involving first-order
derivatives of the velocity modulus multiplied
with [ ]
j
D , is modelled as a source term. The
system of discretized equations is solved by
coupling kinetic and kinematic equations,
accounting for the non-linear constitutive
hypothesis and considering the corresponding
boundary and initial conditions. Since the
implicit set of equations is written
simultaneously for all boundary and internal
nodes, this procedure results in a very large
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. To improve the
economics of the computation, the subdomain
technique is used, where the entire solution
domain is partitioned into subdomains to
which the same described numerical procedure
can be applied. The final system of equations
for the entire domain is then obtained by
adding the sets of equations for each
subdomain considering the compatibility and
equilibrium conditions between their
interfaces, resulting in a much sparse system
matrix suitable to solve with iterative
techniques. In our case each subdomain
consists of four discontinuous 3-node
quadratic boundary elements, and 9-node
corner continuous quadratic internal cell, (Jecl
& Škerget, 2003).
4. TESTNI PRIMER
U činkovitost metode ter pravilnost
delovanja predlagane numeri čne sheme,
temelje če na dodanem Forcheimerjevem
členu, smo testirali na primeru naravne
konvekcije v porozni kotanji, greti od strani.
Kotanja je zapolnjena s homogeno,
nedeformabilno porozno snovjo, popolnoma
zasi čeno z Newtonsko teko čino. Leva stena
kotanje je greta, desna hlajena, horizontalni pa
toplotno izolirani. Za porozno snov
predpostavimo, da ima konstantno poroznost
in prepustnost ter da je hidrodinami čno in
toplotno izotropna.
4. TEST CASE
To test the validity and accuracy of the
proposed numerical scheme, which includes
the additional Forcheimer term, the problem of
natural convection in a porous cavity heated
from the side was investigated. The cavity is
filled with the homogeneous, nondeformable
porous media saturated with Newtonian fluid.
Thr left vertical wall of the cavity is
isothermally heated and the right is
isothermally cooled while the horizontal walls
are adiabatic. The porosity and permeability of
the porous media are assumed to be uniform
throughout the system.

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
12
Teko čina, s katero je porozna snov
zasi čena, je v toplotnem ravnovesju s trdno
fazo, kar pomeni, da lahko toplotno obnašanje
porozne snovi opišemo z eno samo ena čbo za
povpre čno temperaturo. Geometrija ter
hidrodinami čni in temperaturni robni pogoji
obravnavanega problema so podani na sliki 1.
Izra čune smo opravili za dva primera, in sicer
za pravokotno kotanjo z razmerjem stranic
( D H A = ) 1 A = in za vertikalno kotanjo pri
5 A = .
Furthermore, the solid and the fluid phase
are assumed to be in the local thermal
equilibrium, therefore the thermal behaviour of
the porous media can be described by a single
equation for the average temperature. Detailed
presentations of the geometry and boundary
conditions are given in Figure 1. The
calculations are performed for a square cavity
with aspect ratio ( D H A = ) 1 A = and for the
tall cavity with aspect ratio 5 A = .
Slika 1. Geometrija ter robni pogoji.
Figure 1. Geometry and boundary conditions.
Ostali parametri problema so še Darcyjevo
število ) D K ( Da
2
Λ = , teko činsko
Prandtlovo število 
f f f
a Pr ν = , modificirano
Rayleighovo število λ Da Ra Ra
*
= , kjer je
Ra teko činsko Rayleighovo število,
f f
3
T
a T D g Ra ν ∆ β = , λ razmerje toplotnih
prevodnosti, 
e f
λ λ λ = in σ razmerje
toplotnih kapacitet. Zaradi poenostavitve smo
upoštevali 1 Pr
f
= , 1 = λ in 1 = σ . Za
pravokotno kotanjo je poroznost φ enaka 0,5,
Forcheimerjev koeficient 
F
c , dolo čen po
Ergunovem modelu, je 0,404. Uporabili smo
dve neekvidistantni ra čunski mreži 20 x 20 m
in 40 x 40 m podobmo čij, pri čemer je bilo
razmerje med najdaljšo in najkrajšo stranico
6 r = za prvo in 12 r = za drugo mrežo.
Redkejša mreža daje namre č primerljive
Other relevant governing parameters for the
present problem are Darcy number
) D K ( Da
2
Λ = , fluid Prandtl number
f f f
a Pr ν = modified Rayleigh number
λ Da Ra Ra
*
= , where Ra represents fluid
Rayleigh number, 
f f
3
T
a T D g Ra ν ∆ β = , λ is
the ratio of heat conductivity, 
e f
λ λ λ = and
σ is the heat capacity ratio. For
simplification: 1 Pr
f
= , 1 = λ and 1 = σ . In
square cavity porosity φ equals 0.5 and
coefficient 
F
c , is 0.404 as calculated by using
the Ergun model. Two non-uniform
computational meshes of 20 x 20 m and 40 x
40 m subdomains were used with a ratio 6 r =
for the first and 12 r = for the second mesh
with r indicating the ratio between the longest
and the shortest element. This is due to the fact
0
y
T
=
∂
∂
,
y
0
y
T
=
∂
∂
,
D
x
H g
  porozna snov/
porous medium
1 T
h
= 0 T
c
=
0 v
i
= 0 v
i
=
0 v
i
=
0 v
i
=

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
13
rezultate za ve čja Darcyjeva števila,
3
10 Da
−
≥ , ko pa se Da zmanjšuje,
potrebujemo gostejšo mrežo zaradi velikih
hitrostnih gradientov v bližini vertikalnih sten.
Časovni koraki so se zmanjševali z
zmanjševanjem Darcyjevega števila od
s 10 t
16
= ∆ do s 10 t
5 −
= ∆ , konvergen čni
kriterij pa je bil 
6
10 5
−
× = ε . Za visoko
kotanjo smo upoštevali poroznost 8 0, = φ in
koeficient 2 0 c
F
, = . Neekvidistantna ra čunska
mreža je bila v tem primeru sestavljena iz
40 30 × podobmo čij z razmerji 30 r
x
= in
20 r
y
= .
that the first mesh provides results in good
agreement with the published values for
3
10 Da
−
≥ and the second mesh is required for
smaller values of Da because of the steep
velocity gradients near the vertical walls. Time
steps ranging from s 10 t
16
= ∆ to s 10 t
5 −
= ∆
have been employed and the convergence
criterion is determined to be 
6
10 5
−
× = ε . For
tall cavity the porosity is 8 . 0 = φ and
coefficient 2 . 0 =
F
c . A nonuniform
computational mesh of 40 30 × subdomains
was used in this case with 30 r
x
= and
20 r
y
= .
Preglednica 1. Skupni prenos toplote skozi kotanjo, Nu.
Table 1. Total heat transfer across the cavity, Nu .
razmerje stranic Da
10
-1
10
-2
10
-3
10
-4
10
-5
Aspect ratio
*
Ra
mreža/mesh
sedanji rezultati
20 x 20
1.084 1.647 2.322 2.762 2.588
Current results
40 x 40
1.085 1.655 2.330 2.783 2.902
Lauriat & Prasad (1987)
41 x 41
/ 1.70 2.41 2.84 3.02
Jecl & Škerget (2003)
20 x 20
1.086 1.685 2.402 2.823 2.67
100
Jecl & Škerget (2003)
40 x 40
1.088 1.695 2.414 2.847
2.995
sedanji rezultati
Current results 40 x 40
1.660 3.021 4.844 6.633 7.812
500
Jecl & Škerget (2003)
40 x 40
1.681 3.145 5.235 7.185 8.428
sedanji rezultati
Current results
40 x 40
2.07 3.79 6.57 9.31 10.61
A = 1
1000
Jecl & Škerget (2003)
40 x 40
2.08 4.03 6.98 9.89 11.45
sedanji rezultati
Current results
30 x 40
5.299 7.033 8.834 9.609 /
A = 5 100
Jecl & Škerget (2003)
30 x 40
5.312 7.254 9.134 9.953 /
Za prikaz vpliva dodatnega
Forcheimerjevega člena na prenos toplote
skozi kotanjo smo Nusseltova števila,
()
∫ =
∂ ∂ − =
1
0
0 x
dy x T Nu , ki predstavljajo skupni
prenos toplote skozi kotanjo in so izra čunana z
uporabo Brinkmanovega modela (Lauriat &
To demonstrate the effect of the additional
Forcheimer term on the heat transfer across the
cavity, the overall Nusselt number
representing the total heat transfer across the
cavity, ()
∫ =
∂ ∂ − =
1
0
0 x
dy x T Nu , calculated by
using the Brinkman model (Lauriat & Prasad,

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
14
Prasad, 1987; Jecl & Škerget, 2003), primerjali
z Nusseltovimi števili, izra čunanimi na podlagi
Forcheimerjevega modela. Primerjava je
podana v preglednici 1. Iz nje je razvidno, da
je razlika v skupnem prenosu toplote skozi
kotanjo, izra čunana po obeh modelih,
minimalna in nikoli ne preseže 7 %. Ve čje
razlike se pri čakovano pojavijo s
pove čevanjem modificiranega Rayleighovega
števila in z zmanjševanjem Darcyjevega
števila. Ugotovimo lahko tudi, da
Forcheimerjev model daje nižje vrednosti
skupnega prenosa toplote skozi kotanjo v
primerjavi z Brinkmanovim modelom.
Rezultati se ujemajo z ugotovitvami, podanimi
v Lauriat & Prasad (1989) za vertikalno
porozno kotanjo kot tudi z zaklju čki v Lage
(1992) za horizontalni porozni sloj.
Spremembe Nusseltovega števila za oba
modela glede na razli čna Darcyjeva in
modificirana Rayleighova števila za primer
kvadratne kotanje so grafi čno podane na sliki
2.
1987; Jecl & Škerget, 2003), is compared with
the Nusselt number gotten by using the above
described numerical procedure based on the
Forcheimer model. The comparison is given in
Table 1. It shows that the difference in the heat
transfer rate as calculated with both models is
minimal and never exceeds 7%. Higher
distinction is, as expected, obtained with an
increase in the modified Rayleigh numbers and
decrease in the Darcy number. Table 1 further
shows that the Forcheimer model predicts
lower total heat transfer across the cavity than
the Brinkman model. These results are in
agreement with the conclusions reported in
Lauriat & Prasad (1989) for the vertical porous
cavity and also with the conclusions reported
in Lage (1992) for the horizontal porous layer.
The variation of the overall Nusselt number
for both models considering different Darcy
and modified Rayleigh numbers for the case of
square cavity are represented in Figure 2.
0
5
10
15
0.00001 0.0001 0.001 0.01 0.1
Darcyjevo število - Da
Nusseltovo število - Nu
Ra*=100, Brinkman Ra*=100, Forcheimer Ra*=500, Brinkman
Ra*=500, Forcheimer Ra*=1000, Brinkman Ra*=1000, Forcheime
Slika 2. Sprememba Nusseltovega števila pri razli čnih 
*
Ra in Da za oba modela.
Figure 2: Variation of overall Nusselt number with different 
*
Ra and Da for both models.

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
15
01 02 03 0
0 0.1 0.2 0.3 0.4 0.5
x
v-y
Brinkman, Da=1E-2 Forcheimer, Da=1E-2
Brinkman, Da=1E-4 Forcheimer, Da=1E-4
Slika 3. Vertikalni hitrostni profili skozi vodoravno središ čnico
za polovico kvadratne kotanje pri 100 Ra
*
= za oba modela.
Figure 3. Vertical velocity profiles at the horizontal midplane for 100 Ra
*
= for both models.
Iz slike 2 lahko povzamemo, da skupni
prenos toplote skozi kotanjo vedno naraš ča z
modificiranim Rayleighovim številom in pada
glede na Darcyjevo število ter da vklju čitev
Forcheimerjevega člena v gibalno ena čbo ne
spremeni asimptoti čnega obnašanja krivulje
odvisnosti med Nusseltovim in Darcyjevim
številom.
Vpliv Forcheimerjevega vztrajnostnega
člena je prikazan še na sliki 3, na kateri so
podani vertikalni hitrostni profili skozi
vodoravno središ čnico leve polovice kotanje,
izra čunani po obeh modelih za dve razli čni
Darcyjevi števili 
2
10 Da
−
= in 
4
10 Da
−
= za
primer kvadratne kotanje, za 1 A = in
100 Ra
*
= . Pri majhnem Darcyjevem številu
ima porazdelitev hitrosti strme gradiente v
bližini vertikalne stene, z naraš čanjem
Darcyjevega števila se hitrost zmanjšuje,
mesto nastopa maksimalne hitrosti se pomakne
proti središ ču kotanje, pri čemer pa lahko
vidimo, da je razlika v hitrostih za oba modela
From Figure 2 we can observe that the total
heat transfer across the cavity always increases
with the modified Rayleigh number, but the
effect of the Darcy number is just the reverse.
The inclusion of the additional Forcheimer
term in the equation of motion does not change
the asymptotic behaviour of the Nusselt
number versus the Darcy number curve.
The effect of the Forcheimer additional
term is illustrated also in Figure 3, where the
vertical velocity profiles at the horizontal
midplane for the left half of the cavity are
presented for both Brinkman and Forcheimer
model for two different Darcy numbers,
2
10 Da
−
= and 
4
10 Da
−
= in the case of tall
cavity 5 A = and for 100 Ra
*
= . For a small
Darcy number the vertical velocity distribution
shows its largest gradient near the vertical
walls. But when Da is increased the
maximum velocity reduces and the velocity
peaks move away from the walls. It can be
seen that the difference between both models

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
16
minimalna.
Vsi izra čuni so bili opravljeni na delovni
postaji Hewlett - Packard B1000 s procesorjem
PARISC 300 MHz  in 1 GB pomnilnika.
Uporabljen je bil operacijski sistem HP-UX
10.20. Ra čunski časi za predstavljene primere
se gibljejo od nekaj sekund za kvadratno
kotanjo pri redkejši mreži 20 x 20 m in
relativno velikem Darcyjevem številu
1
10
−
= Da do treh dni in pol (81 ur) za
najzahtevnejši primer visoke kotanje pri mreži
30 x 40 m za 
4
10
−
= Da . Ti časi veljajo za
Brinkmanov model, medtem ko so časi ob
dodanem Forcheimerjevem členu daljši, in
sicer v povpre čju za 35 %.
is minimal.
The calculations were obtained by using the
Hewlett - Packard workstation B1000 with
PARISC 300 MHz processor and 1 GB of
memory. The operating system was HP-UX
10.20. The computational time for the
presented test cases ranged from several
seconds for square cavity with 20 x 20 m
subdomains and relatively high Darcy number
1
10
−
= Da to three and a half days (81 hours)
for the most demanding case of tall cavity with
the mesh consisting of 30 x 40 m subdomains
and for 
4
10
−
= Da . Those values are valid for
the Brinkman model while in the case of the
additional Forcheimer term included, the
computational times are greater by about 35
%.
5. ZAKLJU ČKI 5. CONCLUSIONS
Prikazana je uporaba robnoobmo čne
integralske metode (ROIM) za primer naravne
konvekcije v porozni kotanji, greti od strani.
Za prikaz vpliva Forcheimerjevega
vztrajnostnega člena je uporabljena s
Forcheimerjevim členom razširjena Darcy-
Brinkmanova ena čba, rezultati pa so
primerjani z vrednostmi, dobljenimi na osnovi
Brinkmanovega modela.
Vklju čitev Forcheimerjevega člena
povzro ča minimalno zmanjšanje prenosa
toplote skozi kotanjo in tudi minimalno
zmanjšanje hitrosti, ne spremeni pa se
asimptoti čna oblika krivulje odvisnosti med
Nusseltovim in Darcyjevim številom.
Numeri čna shema je zaradi dodanega
nelinearnega člena še bolj kompleksna, kar
vodi do povpre čno 35 % pove čanega
ra čunskega časa, potrebnega za konvergenco.
Na podlagi teh rezultatov je mogo če
potrditi, da v obravnavanem primeru
Forcheimerjev člen ne vpliva bistveno na
skupni prenos toplote, kot je bilo to že
objavljeno v literaturi, kjer so bile za izra čun
uporabljene druge numeri čne metode.
The problem of natural convection in
porous cavity heated from the side is
investigated utilizing the Boundary Domain
Integral Method (BDIM). The Brinkman-
Forcheimer extended Darcy momentum
equation is used to examine the influence of
the additional Forcheimer inertia term and the
results are compared with reported values
obtained based on the Brinkman model.
The inclusion of the Forcheimer term in the
momentum equation leads to a minimal
reduction of the heat transfer rate and velocity
but does not change the asymptotic behaviour
of the Nusselt number versus the Darcy
number curve. The numerical code becomes
more complex and also the computation time
required to achieve convergence is increased
by about 35 %.
Therefore it is possible to conclude that in
the range covered by the present investigation
the Forcheimer term is not really relevant to
the calculation of the global heat transfer
parameter, as reported in the literature in
which the calculations were performed with
other numerical methods.

Jecl, R., Škerget, L., Kramer, J.: Primerjava med Forcheimerjevim in Brinkmanovim modelom konvektivnega toka v
porozni kotanji z robnoobmo čno integralsko metodo – Comparison between the Forcheimer and the Brinkman model
for convective flow in a porous cavity with Boundary-Domain Integral Method
© Acta hydrotechnica 23/38 (2005), 1–17, Ljubljana
17
VIRI – REFERENCES
Bear, J., Bachmat, Y. (1991). Introduction to Modelling of Transport Phenomena in Porous Media,
Kluwer Academic Publishers, Dordrecht.
Jecl, R., Škerget, L., Petrešin, E. (2001). Boundary domain integral method for transport
phenomena in porous media, Int. Journal for Numerical Methods in Fluids 35, 39–54.
Jecl, R., Škerget, L. (2003). Boundary element method for natural convection in non-Newtonian
fluid saturated square porous cavity, Eng. Anal. with Boundary Elements 27, 963–975.
Lage, J. L. (1992). Comparison between the Forcheimer and the convective inertia terms for Benard
convection within a fluid saturated porous medium, HTD 193, Fundamentals of Heat Transfer in
Porous Media, ASME, 49–55.
Lauriat, G., Prasad, V. (1987). Natural convection in a vertical porous cavity: a numerical study for
Brinkman-extended Darcy formulation, Journal of Heat Transfer 109, 688–696.
Lauriat, G., Prasad, V. (1989). Non-Darcian effects on natural convection in a vertical porous
enclosure, International Journal of Heat and Mass Transfer 32, 2135–2148.
Nield, D. A., Bejan, A. (1999). Convection in porous media (Second ed.), Springer-Verlag, New
York.
Škerget, L., Hriberšek, M., Kuhn, G. (1999). Computational fluid dynamics by boundary domain
integral method, Int. Journal for Numerical Methods in Engineering 46, 1291–1311.
Vasseur, P., Wang, C.H., Sen, M. (1990). Natural convection in an inclined rectangular porous slot:
Brinkman-extended Darcy model, Journal of Heat Transfer 112, 507–511.
Naslovi avtorjev – Authors’ Addresses
izr. prof. dr. Renata Jecl
Univerza v Mariboru – University of Maribor
Fakulteta za gradbeništvo – Faculty of Civil Engineering
Smetanova 17, SI-2000 Maribor, Slovenia
E-mail: renata.jecl@uni-mb.si
prof. dr. Leopold Škerget
Univerza v Mariboru – University of Maribor
Fakulteta za strojništvo – Faculty of Mechanical Engineering
Smetanova 17, SI-2000 Maribor, Slovenia
E-mail: leo@uni-mb.si
Janja Kramer
Univerza v Mariboru – University of Maribor
Fakulteta za gradbeništvo – Faculty of Civil Engineering
Smetanova 17, SI-2000 Maribor, Slovenia
E-mail: janja.kramer@uni-mb.si