International Journal of Thermal Sciences 47 (2008) 1020–1025

www.elsevier.com/locate/ijts

Influence of Darcy number on the onset of convection in a porous layer

with a uniform heat source
Ali Nouri-Borujerdi a , Amin R. Noghrehabadi a,b,c , D. Andrew S. Rees b,∗
a School of Mechanical Engineering, Sharif University of Technology, Tehran, Iran
b Department of Mechanical Engineering, University of Bath, Bath BA2 7AY, UK
c Department of Mechanical Engineering, Chamran University of Ahvaz, Ahvaz, Iran

Received 11 September 2006; received in revised form 21 June 2007; accepted 27 July 2007
Available online 31 August 2007

Abstract
This note considers the effect of different Darcy numbers on the onset of natural convection in a horizontal, fluid-saturated porous layer with
uniform internal heating. It is assumed that the two bounding surfaces are maintained at constant but equal temperatures and that the fluid and
porous matrix are in local thermal equilibrium. Linear stability theory is applied to the problem, and numerical solutions obtained using compact
fourth order finite differences are presented for all Darcy numbers between Da = 0 (Darcian porous medium) and Da → ∞ (the clear fluid limit).
The numerical work is supplemented by an asymptotic analysis for small values Da.
© 2007 Elsevier Masson SAS. All rights reserved.

1. Introduction Gasser and Kazimi [5] conducted a comprehensive study of

the onset of thermal convection in a horizontal porous layer us-
Natural convection has been a subject of intensive research ing a linear stability analysis of the basic nonlinear temperature
in porous media in view of its wide range of application in many distribution which is caused by both internal heat generation
engineering and technological areas. Applications include high and heating from below. Therefore two Rayleigh numbers ap-
performance insulation for buildings and cold storage, the insu- pear, one corresponding to internal heating, the other to the
lation of high temperature gas-cooled reactor vessels, the bury- external temperature gradient. They determined how the criti-
ing of drums containing heat-generating chemicals in the earth, cal internal Rayleigh number varies with the size of the external
regenerative heat exchangers containing porous materials and Rayleigh number and vice versa. When the external Rayleigh
exothermic chemical reactions in packed-bed reactors. number is zero, then the internal Rayleigh number is approxi-
Many authors have considered the conditions for instability mately 470. Rudraiah et al. [6] studied the same problem sub-
in a porous layer heated either from below or by means of in- sequently using trial functions to solve the linearised stability
ternal volumetric heat generation. Horton and Rogers [1] and equations. In a short work, Selimos and Poulikakos [7] ex-
Lapwood [2] were the first to establish analytically the critical tended the analyses of [5] and [6] by including a second diffus-
Rayleigh number for onset of convection in a fluid-saturated ing component and by adopting the Darcy–Brinkman momen-
porous layer heated from below without internal heat gener- tum equations. A comprehensive set of results are given in [7]
ation. Their analysis has since been extended substantially to for a small selection of values of the Darcy number. Vasseur
include other types of modeling of porous media, and to moder- and Robillard [8] considered convection in a layer with unequal
but constant heat fluxes imposed on the boundaries. On sub-
ately and strongly nonlinear situations. The reviews by Rees [3]
tracting out the mean temperature rise, they obtained a stability
and Tyvand [4] may be consulted for further details.
problem similar to that of Gasser and Kazimi [5] except that the
disturbance temperatures satisfy Neumann rather than Dirichlet
* Corresponding author. boundary conditions. Further work has appeared in the litera-
E-mail address: ensdasr@bath.ac.uk (D.A.S. Rees). ture detailing nonlinear effects, the effect of different boundary
1290-0729/$ – see front matter © 2007 Elsevier Masson SAS. All rights reserved.
doi:10.1016/j.ijthermalsci.2007.07.014
 A. Nouri-Borujerdi et al. / International Journal of Thermal Sciences 47 (2008) 1020–1025 1021

Nomenclature

C specific heat μ dynamic viscosity

Da Darcy number (Eq. (6)) μe effective viscosity
g gravity ν kinetic viscosity
k wavenumber of disturbance ε porosity
K permeability ψ streamfunction
L depth of the convection layer Ψ streamfunction disturbance
P pressure θ scaled temperature
q  rate of heat generation Θ temperature disturbance
Ra Darcy–Rayleigh number (Eq. (6)) λ amplification rate of disturbance
t time
Superscripts and subscripts
T temperature
u, v horizontal and vertical velocity  dimensional
x, y horizontal and vertical Cartesian coordinate basic basic state
PM porous media
Greek symbols f fluid
α diffusivity s solid
β coefficient of cubical expansion 0 wall temperature
ρ density . k-derivative
σ heat capacity ratio  y-derivative

conditions and flow within finite cavities; the reader is referred

to Nield and Bejan [9] for further discussion.
On the other hand, for layers filled with a clear fluid, Spar-
row et al. [10] performed an analytical study of the thermal
instability of an internally heated fluid layer both with and with-
out heating from below. They showed that with increasing heat
generation rate the fluid becomes more prone to instability, that
is, the critical Rayleigh number decreases. Takashima [11] ap-
plied linear stability theory to the problem of the stability of
natural convection that occurs in an inclined fluid layer with
uniformly distributed internal heat sources and with constant
and equal boundary temperatures. The marginal stability crite- Fig. 1. Definition sketch of the horizontal porous layer with the coordinate sys-
ria for different ranges of Prandtl number and angles of inclina- tem.
tion are reported.
The purpose of this short note is to determine how the critical model, and, subject to the Boussinesq approximation, the full
Rayleigh number and the corresponding wavenumber vary with two-dimensional governing equations take the form,
Darcy number for an internally heated porous layer. As such
it gives the full transition between convection when Darcy’s ∂ û ∂ v̂
law applies and when the full Navier–Stokes equations apply. + =0 (1a)
∂ x̂ ∂ ŷ
Comparison is then made between our results and those for the  2 
Darcy-flow and clear fluid limits. We provide a highly accurate μ ∂ P̂ ∂ û ∂ 2 û
û = − + μe + (1b)
set of numerical results and supplement this with an asymptotic K ∂ x̂ ∂ x̂ 2 ∂ ŷ 2
theory for small values of the Darcy number.  2 
μ ∂ P̂ ∂ v̂ ∂ 2 v̂
v̂ = − + μe + + ρgβ(T − T0 ) (1c)
K ∂ ŷ ∂ x̂ 2 ∂ ŷ 2
2. Governing equations and basic solution  
∂T ∂T ∂T
(ρC)PM + (ρC)f û + v̂
We consider an infinite porous layer confined between two ∂ tˆ ∂ x̂ ∂ ŷ
parallel rigid plates which are separated by a distance L, as  2 2 
∂ T ∂ T
depicted in Fig. 1. It is assumed that the fluid layer is heated = kPM 2
+ 2 + q  (1d)
∂ x̂ ∂ ŷ
internally by a uniform heat sources of strength q  and the two
bounding surfaces are each maintained at the constant temper- where x and y are the horizontal and vertical coordinates and u
ature T0 . The governing equations of motion of fluid in a ho- and v are the corresponding velocity components. All the other
mogeneous and isotropic porous medium follow the Brinkman terms have their usual meaning for porous medium convection,
1022 A. Nouri-Borujerdi et al. / International Journal of Thermal Sciences 47 (2008) 1020–1025

and are given in the Nomenclature. The appropriate boundary which are to be solved subject to the boundary conditions,
conditions are,
1
L ψ = ψy = θ = 0 on y = ± (10)
û = v̂ = T = 0 on ŷ = ± (2) 2
2
The basic steady state now consists of no flow and the following
Here we have taken fixed temperature boundary and no-slip
parabolic temperature profile,
boundary conditions. Other choices of boundary condition have
been made in the published literature, including having one 1 1
θbasic = − y 2 + (11)
surface cooled with the other insulated and/or stress free con- 2 8
ditions. Our choice of thermal conditions means that the upper
half of the layer is unstably stratified, and we have an exam- 3. Linear stability theory
ple of penetrative convection, where disturbances may penetrate
into the lower, stably stratified region below. We may assess the stability characteristics of the evolving
Eqs. (1a–d) may be nondimensionalised using the following basic state using a straightforward perturbation theory. There-
substitutions, fore we set,
L2 σ αPM ψ = Ψ (y)eλt cos kx
tˆ = t, (x̂, ŷ) = L(x, y), (û, v̂) = (u, v) and
αPM L
θ = θbasic (y) + Θ(y)e sin kx
λt
(12)
αPM μ q  L2
P̂ = P, T = T0 + θ (3) where Ψ and Θ are both of a sufficiently small amplitude that
K kPM
where nonlinear terms may be neglected. The value, k, is the hori-
zontal wavenumber of the disturbances. We assume that the
(ρC)PM kPM
σ= , αPM = principle of exchange of stabilities applies (see, for example,
(ρC)f (ρC)f Drazin and Reid [12], and the discussion in Appendix A), in-
kPM = εkf + (1 − ε)ks , q  = εqf + (1 − ε)qs (4) dicating that the onset of convection is stationary, or, in other
These transformations yield the following system of equa- words, that neither travelling nor standing waves appear. There-
tions, fore the following system of linearised disturbance equations
are obtained,
∂u ∂v
+ =0 (5a)
∂x ∂y −Da(Ψ IV − 2k 2 Ψ  + k 4 Ψ ) + Ψ  − k 2 Ψ = Ra kΘ (13a)
 2 
∂P ∂ u ∂ 2u 
Θ −k Θ 2 
+ kΨ θbasic =0 (13b)
u=− + Da + (5b)
∂x ∂x 2 ∂y 2
 2  and this system is subject to the boundary conditions,
∂P ∂ v ∂ 2v
v=− + Da + + Ra θ (5c) 1
∂y ∂x 2 ∂y 2 Ψ = Ψ  = Θ = 0 on y = ± . (14)
2
∂θ ∂θ ∂θ ∂ 2θ ∂ 2θ
+u +v = 2 + 2 +1 (5d) In Eqs. (13a,b) primes denote differentiation with respect to y.
∂t ∂x ∂y ∂x ∂y
These equations cannot be solved analytically and therefore nu-
In above equations the nondimensional parameters, Da, and Ra merical methods must be employed.
are defined according to,
μe K gβKq  L3 4. Numerical simulations
Da = and Ra = (6)
μL2 ναPM kPM
The boundary conditions are now, Eqs. (13a,b) form an ordinary differential eigenvalue prob-
lem for Ra as a function of Da and the wavenumber, k. When
1
u = v = θ = 0 on y = ± (7) Da = 0 we recover the Darcy-flow case considered by Gasser
2 and Kazimi [5]. In this paper we solve the full system (13)
From the continuity equation, (5a), a streamfunction ψ may be using a direct method related closely to the one described
defined according to, in Rees [13]. Eqs. (13a,b) were reduced to a set of three
∂ψ ∂ψ second-order equations by introducing a vorticity-like variable,
u=− and v = (8)
∂y ∂x and then discretised using fourth order compact differences
After the elimination of the pressure P between Eqs. (5b) (Spotz [14]) on a uniform grid in the y-direction. The zero nor-
and (5c), Eqs. (5a–d) reduce to the system, mal flow, tangential flow and temperature conditions provide a
 4  sufficient number of boundary conditions for these equations.
∂ ψ ∂ 4ψ ∂ 4ψ ∂ 2ψ ∂ 2ψ ∂θ However, the eigenvalue, Ra, also needs to be found, and this
−Da 4
+ 2 2 2
+ 4
+ + = Ra
∂x ∂x ∂y ∂y ∂x 2 ∂y 2 ∂x requires one more condition; this is provided by the following
(9a) normalization condition,
∂θ ∂ψ ∂θ ∂ψ ∂θ 2
∂ θ 2
∂ θ
+ − = + +1 1
(9b) Θ = 1 on y = (15)
∂t ∂x ∂y ∂y ∂x ∂x 2 ∂y 2 2
 A. Nouri-Borujerdi et al. / International Journal of Thermal Sciences 47 (2008) 1020–1025 1023

The resulting discretised system is then solved using a standard

multidimensional Newton–Raphson iteration technique. The it-
eration matrix takes a block tridiagonal form where there is one
further column and row of nonzero blocks, and therefore the
block-Thomas algorithm was modified to account for these ex-
tra blocks; see Eq. (11) in Rees [13] where the full procedure is
described in more detail.
Numerical experiments indicate that the neutral stability
curves always have the same qualitative form, namely that
there is one minimum value of Ra at a critical value of k with
monotonic growth towards infinity as k → 0 and as k → ∞.
Therefore we concentrate solely on these critical values, since
such a minimum value of Ra signifies the point above which (a)
we may expect convection to take place in an infinite layer. This
minimization was achieved by insisting that ∂Ra/∂k = 0 and by
supplementing Eqs. (13a,b) with their derivatives with respect
to k. If we define the variables,
∂Ψ ∂Θ
Ψ̇ = and Θ̇ = (16)
∂k ∂k
then differentiation of Eqs. (13a,b) with respect to k yields the
following system,
−Da(Ψ̇ IV − 2k 2 Ψ̇  + k 4 Ψ̇ ) + Ψ̇  − k 2 Ψ̇ − Ra k Θ̇
= Da(−4kΨ  + 4k 3 Ψ ) + 2kΨ + Ra Θ (17a)
  
Θ̇ − k Θ̇
2
+ k Ψ̇ θbasic = 2kΘ − Ψ θbasic (17b)
subject to the boundary conditions, (b)
1 Fig. 2. Disturbance streamlines and isotherms corresponding to the Da = 0.1 in
Ψ̇ = Θ̇ = 0 on y = ± (18)
2 the (x, y)-plane: (a) streamlines, (b) isotherms.
As the wave number is now a second eigenvalue, we need to im-
pose a second normalization condition that Θ̇  ( 12 ) = 1, although
any other value of this derivative yields precisely the same val-
ues of Ra and k. This new extended system now consists of six
second order equations and two normalisation conditions. The
block matrices which appear in the Newton–Raphson iteration
matrix are now 6 × 6.
There is now only one parameter to vary, namely Da, and
solutions are presented for the range, 10−6  Da  102 . Uni-
form grids were used in the computation and 200 intervals in
the range −0.5  y  0.5 were used as a basic grid. Grid re-
finement was then used and the accuracy was improved further
using Richardson’s extrapolation technique to obtain over 10
significant figures of accuracy, even for values of the Darcy
number as small as 10−6 .
Figs. (2a) and (2b) show the respective streamlines and
isotherms of the disturbance shapes, Ψ cos kx and Θ sin kx, for
Da = 0.1. Each frame displays contours corresponding to 20
equally-spaced subintervals between their respective maxima
and minima. The streamlines and isotherms are displayed in
(x, y) space. Given that it is only the top half of the layer that
is unstably stratified, the disturbance temperature field is con-
centrated within this half with a pair of weak cells in the lower Fig. 3. Variation of a normalized horizontal velocity, −Ψ  with y for different
values of Da.
half. The streamlines also display a bias towards the upper half
of the channel.
Fig. 3 shows normalized profiles of −Ψ  , which is related which corresponds to the clear fluid limit, the velocity profile
to the horizontal fluid velocity, in order to see how these pro- varies relatively slowly across the cavity. At such large val-
files vary with the value of Da. When Da is relatively large, ues of Da, the maximum absolute velocity occurs at y = 0.35,
1024 A. Nouri-Borujerdi et al. / International Journal of Thermal Sciences 47 (2008) 1020–1025

Fig. 5. Variation of critical wavenumber, k, with Log10 Da for 10−6 

Da  102 . The dashed line corresponds to the small-Da solution given by
Fig. 4. Variation of critical Rayleigh number with Log10 Da for 10−6  Eq. (19b).
Da  102 . Values are presented in terms of the porous Rayleigh number, Ra,
and the clear fluid Rayleigh number Ra/Da. Dashed lines correspond to the
small-Da solution given by Eq. (19a). Rac = 471.384663(1 + 3Da1/2 + 70.6226324Da + · · ·) (19a)
kc = 4.67518897(1 + Da1/2 + · · ·) (19b)
which is well within the upper half of the layer. However, as
where our computed data is correct to more than 10 significant
Da decreases towards zero, the position at which the maximum
figures. We note that the value, Ra ∼ 469, which may be de-
velocity occurs rises towards the upper surface. A very distinct
rived from the data presented in Gasser and Kazimi [5], is close
boundary layer is formed at this upper surface within which
to the leading term in (19a), but it was obtained using a small
the velocity changes rapidly to zero; a similar though weaker
number of terms in a Galerkin expansion. We see from Fig. 5
boundary layer is formed at the lower surface. It may be shown
that the above expression for kc is accurate only for Da < 10−4 .
using a straightforward order-of-magnitude argument that the
However, Fig. 4 shows that Eq. (19a) is very accurate indeed for
boundary layers are of thickness O(Da−1/2 ). A detailed asymp-
the whole range of values of Da, even though it is the result of
totic analysis of this phenomenon as it applies to Darcy–Bénard
a small-Da analysis. In particular, when Da is large, (19a) gives
convection was given by Rees [13], and we shall present the re-
Rac /Da ∼ 33290.43, which is very close indeed to the results of
sults of a similar analysis below.
both Sparrow et al. [4] and Takashima [7], Rac /Da = 37325.17,
The respective variations in the critical values of Ra and k an error of only 12%. Indeed, if we were to make an ad hoc
with Da are shown in Figs. 4 and 5. Fig. 4 shows the crit- modification to (19a) by taking the large-Da value of Ra into
ical Rayleigh number in two forms, Ra, the porous medium account:
Rayleigh number, and Ra/Da, the clear fluid Rayleigh num-
ber. The transition between the porous medium limit and the Rac = 471.384663(1 + 3Da1/2 + 79.1820Da) (20)
clear fluid limit is smooth with Ra increasing monotonically then this formula is in error by less than 1% over the whole
and Ra/Da decreasing monotonically. With regard to the for- range of Da values, with the largest error occurring at Da ∼ 0.1.
mer, it is to be expected that Ra should rise because viscous Eq. (20) may therefore be treated as a good correlation.
effects, as mediated by the Brinkman terms, increase in severity
as Da increases. On the other hand, the variation of the criti- 5. Conclusion
cal wavenumber is not monotonic, a property it shares with the
Darcy–Bénard problem (see [13]), but there is little overall vari- In this short paper we have determined how the presence
ation. of the Brinkman terms affects the onset criterion for the sta-
The approach to the Darcy limit is seen clearly in these fig- bility of natural convection in a horizontal porous layer with
ures, where the dashed curves represent a small-Da asymptotic uniform heat generation and the standard no-slip boundary con-
analysis. This analysis follows precisely the one described in ditions. A smooth monotonic variation in the critical Rayleigh
detail in Rees [13] for the Darcy–Bénard problem, except that is found, with the porous and clear fluid limits being reproduced
numerical solutions were required to solve the various ordinary very accurately. The variation in the critical wavenumber, k,
differential equations which arise. We find that is not large, but it is not monotonic. In terms of the critical
 A. Nouri-Borujerdi et al. / International Journal of Thermal Sciences 47 (2008) 1020–1025 1025

porous Rayleigh number, we may neglect the Brinkman terms 1/2

  
and safely use Darcy’s law when Da < 10−3.5 . On the other = −λ Ψ Ψ¯  + k 2 Ψ Ψ̄ dy (A.3)
hand, Da > 1 reproduces the clear fluid limit with a high de- −1/2
gree of accuracy. However, the formula given in (20) may be
All the integrals in Eq. (A.3) are strictly real, and therefore λ
used over the whole range of Da with less than 1% error.
must take real values. Therefore critical values of Ra corre-
spond only to zero values of λ, and the Principle of Exchange
Acknowledgements
of Stabilities applies to the Darcy case.
It is possible to extend this analysis easily to the more gen-
The second author (ARN) wishes to express his cordial
eral Darcy–Brinkman case given by Eqs. (13), but only for
thanks to the British Council for financial support for his visit
stress-free boundary conditions (also see the discussion of Her-
to the University of Bath, and to the University of Bath for their
ron [15] and cited references). Therefore we shall assume that
hospitality. The authors would like to thank the anonymous ref-
the Principle of Exchange of Stabilities is valid.
erees for their constructive comments.
References
Appendix A
[1] C.W. Horton, F.T. Rogers, Convective currents in a porous medium,
J. Appl. Phys. 16 (1945) 367–370.
If we were to replace ψ and θ in Eqs. (13) by Ψ eλt and
[2] E.R. Lapwood, Convection of a fluid in a porous medium, Proc. Camb.
Θeλt , where λ is the exponential growth rate of disturbances, Phil. Soc. 44 (1948) 508–521.
then Eqs. (13a) and (13b) become, [3] D.A.S. Rees, The stability of Darcy–Bénard convection, in: K. Vafai (Ed.),
Handbook of Porous Media, Marcel Dekker, 2000, pp. 521–558.
Ψ  − k 2 Ψ = Ra kΘ (A.1a) [4] P.A. Tyvand, The onset of Rayleigh Bénard convection in porous bodies,
 in: D.B. Ingham, I. Pop (Eds.), Transport Phenomena in Porous Media II,
Θ − k Θ = kyΨ + λΘ
2
(A.1b)
Pergamon, 2002, pp. 82–112.
When Da = 0 the appropriate boundary conditions are that [5] R.D. Gasser, M.S. Kazimi, Onset of convection in a porous medium with
internal heat generation, J. Heat Transfer (1976) 49–54.
1 [6] N. Rudraiah, B. Veerappa, S. Balachandra Rao, Convection in a fluid sat-
Ψ =Θ =0 on y = ± (A.1c) urated porous layer with non-uniform temperature gradient, Int. J. Heat
2
Mass Transfer 25 (1982) 1147–1156.
We note that these boundary conditions, when combined with [7] B. Selimos, D. Poulikakos, On double diffusion in a Brinkman heat gen-
Eq. (A.1a), yield the fact that Ψ  = 0 on the boundaries. Our erating porous layer, Int. Comm. Heat Mass Transfer 12 (1985) 149–158.
intention here is to show that λ takes only real values, so that [8] P. Vasseur, L. Robillard, The Brinkman model in a porous layer: effects of
the Principle of Exchange of Stabilities applies. non uniform thermal gradient, Int. J. Heat Mass Transfer 36 (1993) 4199–
4206.
We may eliminate Θ from Eqs. (A.1) to give, [9] D.A. Nield, A. Bejan, Convection in Porous Media, Springer-Verlag, New
Ψ  − 2k 2 Ψ  + (k 4 − Ra k 2 )Ψ = λ(Ψ  − k 2 Ψ ) (A.2a) York, 2006.
[10] E.M. Sparrow, R.J. Goldstein, V.K. Jonsson, Thermal instability in a hor-
which is subject to the boundary conditions, izontal fluid layer: effect of boundary conditions and non linear tempera-
ture profile, J. Fluid Mech. 18 (1964) 513–528.
1
Ψ = Ψ  = 0 on y = ±
[11] M. Takashima, The stability of natural convection in an inclined fluid layer
(A.2b) with internal heat generation, J. Phys. Soc. Japan 58 (1989) 4431–4440.
2
[12] P.G. Drazin, W.H. Reid, Hydrodynamic Stability, second ed., Cambridge
On taking Ψ to be complex in general, we may multiply University Press, Cambridge, UK, 2004.
Eq. (A.2a) by Ψ̄ and perform a sufficient number of integra- [13] D.A.S. Rees, The onset of Darcy–Brinkman convection in a porous layer:
tions by parts to obtain the following, an asymptotic analysis, Int. J. Heat Mass Transfer 45 (2002) 2213–2220.
[14] W.F. Spotz, High-order compact finite difference schemes for computa-
1/2 tional mechanics, PhD thesis, University of Texas at Austin, Austin, TX,
   
Ψ Ψ̄ + 2k 2 Ψ  Ψ¯  + (k 4 − Ra k 2 )Ψ Ψ̄ dy 1995.
[15] I.H. Herron, Onset of convection in a porous medium with internal heat
−1/2 source and variable gravity, Int. J. Eng. Sci. 39 (2001) 201–208.