Transport in Porous Media 17: 271-283, 1994. 9 1994 KluwerAcademic Publishers. Printed in the Netherlands.
271
Convection in Anisotropic Inclined Porous Layers MARK TREW 1 and R O B E R T M c K I B B I N 2 1Department of Engineering Science, School of Engineering, University of Auckland, Private Bag, Auckland, New Zealand 2Department of Mathematics, Massey University, Private Bag, Palmerston North, New Zealand (Received: 13 November 1993; in final form: 15 November 1994) Abstract. Two-dimensional convective flow in a system of fluid-saturated porous layers is considered for small fluid velocities. The layers are inclined at an angle to the horizontal and there is a temperature gradient across them. A mathematical model of the system is constructed and the resulting series solution is numerically summed. The unique features of this model include its ability to represent non-symmetric multi-layered systems and materials whose permeabilities are anisotropic. The model is applied to various physical configurations. System aspect ratios, system symmetry and permeablity anisotropies are found to influence the effect that layers of different materials have on the convective flow. Key words: Porous media, convection, layering, mathematical model.
Nomenclature
b g H
11/L, aspect ratio of the system gravitational body force per unit mass [ms -2] thickness of the sloping system [m]
K
(I~
L N
length of the sloping system [m] number of layers in the system pressure [kg m- t s-a]
P Ra T
To
IVT I U v x
Ifz0) 'permeabilitytens~
gI(zlaH 2
IV Tdl sin 7, Rayleigh number for the problem
temperature [K] a reference temperature [K] magnitude of the diffusive temperature gradient [Km- 1] non-dimensional component of mean specific discharge in the X direction (u, w), mean specific discharges [ms -1] (x, z), position vector [m]
272 X Z oz
7
r 0 # P p0 k9
MARK TREW AND ROBERT McKIBBIN
non-dimensionalised z non-dimensionalised z volumetric thermal expansion coefficient [K -1] angle between the isotherms and the horizontal permeability anisotropy ratio angle between the sloping system and the horizontal thermal diffusivity of the fluid phase [m2s - 1] dynamic fluid viscosity [kg m-Xs -l] kinematic viscosity [mZs- 1] density of the fluid [kg m -3] density at the reference temperature To [kg m -3 ] dimensional streamfunction [mas -1] non-dimensional streamfunction
Subscripts: z, z denote the z and z directions, respectively i denotes the ith layer of an N layer system 1. Introduction
Thermal convection in layered porous material is of interest in many engineering fields. These include chemical processes, the passage of air though insulating layers in house walls, oil recovery techniques, oil basin analysis and geothermal engineering. In sedimentary basins, reactions between groundwater and solids results in a cycle of dissolution and precipitation of various minerals. When coupled with fluid motion, this diagenetic process creates porosity changes in certain regions of the system. To gain an understanding of the occurrence and location of these porosity changes, it is necessary to first construct a model of the groundwater flow and then to analyse the solute flow in the system. Porosity patterns found from core samples may be linked to particular system geometries if the fluid flow leading to such patterns is known. This is of importance in oil basin analysis. Density gradients and gravitational body forces create buoyancy forces within a fluid and lead to convection. The density gradients are caused by temperature gradients (in the first instance, before solute transport) and geometric factors. Owing to fluid expansion, density normally decreases with increasing temperature (an exception to this is H20 at 0 to 4~ If fluid density variations in a saturated porous material depends only on temperature, then the occurrence of convective flow is determined by the value of the Rayleigh number, a dimensionless parameter based on system properties as well as the temperature gradient. The flow is polycellular with the number of cells formed
CONVECTION IN ANISOTROPIC INCLINED POROUS LAYERS
.... ~i i ~::~i~:~i~!~!~i ! i i ::i ::~ ~.
~
-
~
273 rN
Fig. 1. An inclined stratification plane and isotherms and the resulting system of sloping porous layers that is modelled. being dependent on the most unstable mode for the system (see for example McKibbin & O'Sullivan, 1980). This convection is frequently referred to as Rayleigh convection. However, in an inclined system density isolines are non-horizontal and this inherent instability leads to convective flow regardless of the temperature gradient. This problem was investigated experimentally by Bories & Combarnous (1973). The flow is two-dimensional and unicellular, in the z direction (see Figure 1), and remains unicellular provided the Rayleigh number does not exceed a value of approximately 47r2. Non-horizontal density isolines occur in geological systems where the stratification plane is inclined to the horizontal (see Figure 1). A non-horizontal stratification plane is a common natural occurrence, due to the folding of the earth's crust, so convection resulting from inclined density isolines is generally found more frequently than Rayleigh convection in geological applications. This paper outlines the construction of a model whose origins are based on work conducted by Ludvigsen et al. (1992). The model of low velocity convective flow presented here is generalised to non-symmetric multilayered porous systems and allows for variation of layer permeability in each direction. The geometry on which the model is based is shown in Figure 1.
2. Constructing the Mathematical Model Steady convection in a porous system is represented by equations of conservation of mass, conservation of momentum (Darcy's law), conservatiion of energy and a state equation describing the variation of density with temperature. Using the Boussinesq approximation, the equations are: V.v = 0
(conservation of mass)
(2.1)
v --- ~ ( - V p + pg) # v . V T = ~X72T
(conservation of momentum)
(2.2)
(conservation of energy)
(2.3)
p = po[1 - o~(T - To)]
(linear equation of state)
(2.4)
274
MARK TREW AND ROBERT McKIBBIN
where v is the Darcy velocity, or mean specific discharge, K is the permeability tensor (assumed to be diagonal, with components Kx and L~), # is the dynamic fluid viscosity, p is the pressure, p is the fluid density, g is the gravitational body force per unit mass, T is the temperature, ~ is the thermal diffusivity of fluid, To and P0 are reference values, and a is the volumetric thermal expansion coefficient. For two-dimensional flow, a streamfunction ~b, satisfying (2.1), is introduced such that: u =
0r
~z
0r
(2.5)
a n d w = - - 0---x-"
In sedimentary layers it is observed that fluid velocity is small. Thus, as a first approximation, heat transfer is assumed to occur only by diffusive processes and convective heat transfer is ignored. Solving (2.3) using this approximation gives a linear temperature distribution. This may be expressed as: (2.6)
T = To+ Cx + Dx + E .
where C, D, and E are constants. By considering the isotherms and the diffusive temperature gradient, VTa, it is found that: c -- IVTdl sin(0 - 7)
and
D = IVTdl cos(0 - ~).
(2.7)
where [VTa[ is the magnitude of the diffusive temperature gradient, 0 is the angle the system makes with the horizontal, and 3' is the angle the isotherms make with the horizontal (see Figure 1). The two components of (2.2) are combined to eliminate pressure, and together with (2.4), (2.5), (2.6), and (2.7), yield an expression for layer i of an N layer system of stratified material. 02r 1 o2r -Oz - - T + (i Oz 2 -
g K~,aIVTal sinT,
(2.8)
where ( = K x / K ~ is the permeability anisotropy ratio, u is the kinematic viscosity, and K~ is the permeability in the z direction (see Figure 1). The streamfunction ~b[mZs-1 ] is non-dimensionalised by the thermal diffusivity n[m2s-1], and z[m] and z[m] by II[m], the height of the system of sloping permeable layers. In non-dimensional form, (2.8) becomes 02q~i 1 02f~i gzl ----~ OX + ~i ~OZ - Kz~ Ra'
(2.9)
where Ra = (gKzlaH2/ua)lVTdl sin 7 is the Rayleigh number for the system. This gives a ratio of destabilising to stabilising influences in the system. Nondimensionalising the variables enables magnitude comparisons of solutions to be easily made.
275
CONVECTION IN ANISOTROPIC INCLINED POROUS LAYERS
Conducting a dimensional analysis on (2.3) and (2.8) shows that for heat transfer due to the convection to be negligible, the problem Rayleigh number, Ra, is restricted by:
(2.10)
( b + 1)
where b = lt/L is the aspect ratio of the system. Since the flow is unicellular (Bories et al., 1973), there is zero normal flow on all four sides of the rectangular system. Furthermore, both specific mass flux and pressure must be continuous across the N - 1 layer interfaces at Z = Zi, i = 1, 2 , . . . , N - 1 (Dagan, 1989, pp. 140-141). The first boundary condition is written as: ~=0
onX=O,X=l/b
(2.11)
kg=O
onZ=O,
(2.12)
Z=
1.
For the flow to be continuous, the streamfunction must be continuous across the layer interfaces. 9 i = gli+l
on Z = Zi, i = 1, 2 , . . . , N - 1.
(2.13)
From pressure, temperature and density continuity across the layer interfaces, (2.2) gives
I(x~
-
U +l Kx~+~
on Z =
Zi,
i = 1, 2 , . . . ,
N-
1.
(2.14)
3. S o l u t i o n o f the G o v e r n i n g E q u a t i o n s
The solution to the homogeneous part of (2.9) is found using the method of separation of variables. Since the homogeneous equation is linear, a particular integral satisfying the non-homogeneous part of (2.9) is then added, giving an expression for 9 in layer i:
Kzi Ra t~i(X, Z ) - ~ X ( b X - 1)+ oo
{ An, cosh (nTr~
bZ) + Bn, sinh (nTr v~/b Z) } sin ( nTrbX). (3.1)
n=l
Applying the boundary conditions (2.12)-(2.14) to equation (3.1) gives Fourier series from which expressions for the coefficients An and Bn in each layer are
276
MARK TREW AND ROBERT McK1BBIN n
0
g-
iAl BI A2
0
B2 9
~
0
A3
H .
..~
9
~
0
BN_I
0
AN BN d
Fig. 2.
Assemblyof the FCEs for the nth term (*-non-zeroterms).
obtained9 These expressions or Fourier Coefficient Equations (FCEs) from a system of 2N linear equations for the 2N coefficients of the nth term in the series in equation (3.1). The system of linear equations may be solved for each term and the infinite series approximated by a finite sum of terms. The structure of the system is given in Figure 2. The convergence of the series solution is easily proved for a one-layer model (Trew, 1991)9 Numerical experimentation shows this to be true also for multilayered models. If the system parameters b and ( are chosen to be large, numerical inaccuracy causes the system of FCEs to become rank deficient after a certain number of terms. This situation is improved by scaling the parameters, however, in general the series converges before rank deficiency occurs.
4. Application of the Model The following plots all give non-dimensionalised values. A value of R a = 0.3 was chosen for all the systems considered, following the restriction given in (2.10). The orientation of the flow in the rectangular system is anticlockwise and the streamlines are presented in a square region, regardless of aspect ratio or system inclination9 The profiles shown are of the magnitude of mean specific discharge as a function of Z, for X = 0.0 ( ), 4~ (*-*-*), ~ ( . . . . )" Ludvigsen et al. (1992) considered a symmetrical three-layer system of two identical layers separated by a low permeability layer. They showed that for a sufficiently small aspect ratio, the low permeability layer had almost no effect on the flow. Using the extended model presented here, it is desirable to consider whether this observation concerning the effect of layering on the flow is altered in any way by system asymmetry, systems of more than three layers or layers with anisotropic permeabilities.
277
CONVECTION IN ANISOTROPIC INCLINED POROUS LAYERS Ca) 1.0 . . . . . . . . . . . . . . . . . 0.7500
150000 ~
X~
..................0.25 +- = --0.13 - 0.00 ' "
0,0
"" "~fic'.
Dlschargel"
0.5
(b) 1s 0.7500 0.5000 X~
..................1.00 0.2500 ++ : =0.50 -
0.0
IMeanSpecificDischargel
-
0.00
0.5
(c) I , , . I , , , I , , , I , , ,
l , , , l , , , l ,
I
0.75 Z 0.5 0.25+ 0 , , 0.25
, i . 0.5
+ ,"i . , 0.75
, i , I
, ,
i + , , i , ,".'i + , , 1.25 1.5 1.75
b
Fig. 3. Streamlines, mean specific discharge magnitude profiles and recirculation centre movement for a four layer system. K~ = 1.0, 2.0, 4.0 and 8.0 for layers 1 to 4 respectively. = 4.0, 2.0, 1.0, and 0.5 for layers 1 to 4 respectively, b = 2.0 in (a) and 0.5 in (b). (c) tracks the vertical movement of the centre of recirculation as the aspect ratio of the system is varied.
E x p e r i m e n t a t i o n w i t h d i f f e r e n t s y s t e m s i n d i c a t e s that the a s p e c t r a t i o h a s s o m e i n f l u e n c e o n t h e r e s p o n s e o f t h e f l o w to l a y e r i n g . T h i s is n o t s u r p r i s i n g , s i n c e s y s t e m s w i t h l a r g e a s p e c t r a t i o s (b = II/L large) h a v e m o s t o f the f l o w m o v i n g p a r a l l e l to t h e e n d b o u n d a r i e s a n d t h u s c r o s s i n g l a y e r i n t e r f a c e s w h e r e t h e r e m a y b e g r a d i e n t s o n t h e m e a n s p e c i f i c d i s c h a r g e a n d c o r r e s p o n d i n g c h a n g e s in the f l o w
278
MARK TREW AND ROBERTMcKIBBIN
(a)
I 11~i i i iilii!::: iiiiii Z
X: ..................0.25 : .- ,0,13 0.00
" ~ - ~ '
0.0
"
9
,
,
[MeanSpecificDtschargcl
(b)
1.0
.....
-1
zl~ ................................... 0.5000.................. ?oo 00t~i-
'~~ -
.... , , . . . . . . . . . . . . . . . . . . ... . . . . . . . . . . . . . . . . .
9
'"
IMeanSpecificDischargel
1~176176 ~5 : : ~~176 ooo 9
(c)
0.75
J z
0.5 0.25
0 ''1 0.25 0,5
'1'',1,' 0.75
',1'7'
I
b
1.9.5
1.5
1.75
Fig. 4. Streamlines, mean specific discharge magnitude profiles and recirculation centre movement for a three layer system./fz = 1.0, 0.05, and 1.0 for layers 1 to 3 respectively. = 1.0 for all layers, b = 2.0 in (a) and 0.5 in (b). (c) tracks the vertical movement of the recirculation centres, with varying system aspect ratios, in the first and second layers. Recirculation in the first layer ceases at a system aspect ratio of approximately 0.75.
279
CONVECTION IN ANISOTROPIC INCLINED POROUS LAYERS
.~-.
"i o.i
1.C
ill
!!! I II !!!
9
~.~, '.:: ~-~
I
I
Z X:
.................0.50 -- ; =0.25 - -
0.00
..........................
0.0
IMean Specific Dischargel
1).5
Fig. 5. Streamlinesand mean specificdischargemagnitudeprofilesfor a singlelayer system. K~ = 1.0, ( = lO0.Oandb = 1.0. direction. The problem is similar to the refraction of light due to the light wave moving between materials in which it travels with different velocities. Larger aspect ratio systems tend to move the stagnation points in the flow (at the centre of recirculation rolls) toward regions of lower permeability and thus split the flow, localising it in these regions. In systems with smaller aspect ratios (b = H / L small), most of the flow is parallel to the upper and lower boundaries where it is unaffected by layer interfaces. Smaller aspect ratio systems tend to move stagnation points toward the geometrical centre of the system. This behaviour is seen in Figure 3 and Figure 4. When a single layer of material is more permeable in the x direction than the z (i.e. a high permeability anisotropy ratio), flow prefers to move parallel to the upper and lower boundaries (see Figure 5). Regions of high mean specific discharge form close to these boundaries. There tends to be a region of flow at the upper and lower boundaries where the mean specific discharge is varying across this region and constant throughout the rest of the system at a given X value. This region is reduced by a decreasing aspect ratio and increased for an increasing aspect ratio. With a low permeability in its preferred direction of flow the fluid will reach a maximum discharge closer to the boundary and so form a thin region of varying discharge. A low permeability anisotropy ratio leads to the same behaviour, with a region of varying discharge forming along the end walls, however with reduced overall magnitude of the mean specific discharge. The gradients (e.g. temperature, geometry) required to drive the primarily vertical flow associated with the low permeability anisotropy ratio would be greater than the gradients required to drive primarily horizontal flow. Thus, for a given driving gradient it might be expected that in a system with a low permeability anisotropy ratio there will be a smaller mean specific discharge compared to other systems. Varying the permeability anisotropy ratio of a layer affects the response of the flow to that layer. This is demonstrated by considering the three layer system of Ludvigsen et al. (1992), shown in Figure 6(a). When the middle layer is set at
280
M A R K T R E W A N D ROBERT M e K I B B I N
(a)
Io.6ooo
i~-- .......................................
.......................................
0.4000 .................. 0.Xi0 ;~
IMeanSpeci'filc Discharge'1 ' '
0,(
:Oo:~~
05
(b) ~.~
~
~ 1.0~----~
Z
iiiiii 0.6000 0.4000
X:
..................0.50 :0.25 - 0.00 =
0.0
IMeanSpecificDischargel
).5
(c)
..............If......................10'6000 :i~
.................. 0.50
0 D .
: = :0"25 0.00 t
i
n
,~,
i
5 1 SpeciM fic eDlschargel a n
i
.
CONVECTION IN ANISOTROPIC INCLINED POROUS LAYERS
281
(d)
0 - 7 5 1 ~ Z0 . 5 ~
.......
...........................
............... iiiiiiiiiii
0.25 0.5 0.75 ~ b 1.25 1.5 1.75 2 Fig. 6. Streamlines, mean specific discharge magnitude profiles and recirculation centre movement for a three layer system with varying permeability anisotropy in layer 2. Layers 1 and 3 are identical. In (a),/(~ = 0.05 and ff = 1.0. In (b), Kz = 1.0 and ff = 0.05. In (c), /(z = 0.05 and r = 20.0. The system has b = 1.0. (d) tracks the vertical movement of the centre of recirculation as the aspect ratio is varied for each of the three configurations considered, for ff = 0.05 ( . . . . . ); r = 1.0 ( ); ff = 20.0 ( . . . . ).
the same permeability in the z direction as the bottom and top layers and the permeability in the x direction is low, the flow is almost completely independent of the layer. Most flow crossing this layer does so while moving up or down the end walls. Since this is the preferred direction of flow, the relatively small effect of the layering is explained. However, a layer with a large permeability anisotropy ratio has a slightly increased effect on the flow at the larger aspect ratios. In Figure 6(c), flow is seen to "pinch" toward the centre of the system in the middle layer. The reasoning behind this follows the discussion concerning the single layer of anisotropic material. The fluid attempts to flow parallel to the interface for as far as possible and forms the boundary region of changing discharge within the layer itself, i.e. the discharge gradients tend to smear out into the middle layer itself, rather than remaining across the interface. This smearing-out of the discharge gradient leads to an increased tendency to form the recirculation rolls observed and reduces the tendency of a lower aspect ratio to cause the flow to ignore layering. When the centre of recirculation is tracked (Figure 6(d)), an increasing permeability anisotropy ratio seems to have a limiting effect. A twenty fold increase from ( = 0.05 to ( = 1.0 has a large effect, whereas a further twenty fold increase to ( = 20.0 has a much smaller effect. As a symmetric system is modified so that it becomes non-symmetric (see Figure 7), layering is found to have an increased effect. Asymmetry leads to increased discharge magnitudes which result in increased discharge gradients at layer interfaces. The larger the changes in the magnitude of the mean specific discharge are at the interfaces, the greater the influence of the layer on the flow. The track of the recirculation centre shows that very small system aspect ratios are required for layering to be ignored by the flow.
282
M A R K T R E W AND ROBERT McKIBBIN (a) I.C
..,....,.:..' '
0.6~0
0.4000
0.C
i .....ii.... .....iii .....i..... IMean Specific Dischargel
X: ..................1.00 : ,0.50 0.00
3.5
(b)
1.0
'/
'
X:
o
,ii i iiiii iii
0.3000 ..................1.00 .---.----0.50 3.1000 - 0.00 ).5
(c)
X: ..................50.00 ,--,--,25.00 0.00
21211121ZZIII ~176 0.0
iMean'Specihc Dischargel . . . . .
0.0500 0.5
CONVECTIONIN ANISOTROPICINCLINEDPOROUSLAYERS
283
(d) I
,,,I,,,
.,i
,i.
i
,
,
i
I
,
i
,
I
0.75"
Z 0.5"
0.25'
~
'0!#'~75"I' b'i.~ ?.# L~k"
Fig. 7. Mean specific discharge magnitude profiles and recirculation centre movement for a three layer system with varying symmetry. The middle layer has Kz = 0.05 and ( = 1.0. The top and bottom layers are identical. In (a) and (b), b = 0.5 and in (c), b = 0.01. (d) tracks the vertical movement of the centre of recirculation as the aspect ratio is varied for the configuration in (b). Recirculation in layer one is found to cease when the system aspect ratio is approximately 1.75.
5. Conclusions The model presented in this paper may be applied to a wide variety of physical configurations o f convective flow in sloping systems o f saturated porous material. The applications o f the model considered here are an extension o f the observations made by Ludvigsen et al. (1992) that low permeability layers have a reduced influence on the flow for systems with small aspect ratios. Layers with high permeability anisotropy ratios (>>1) increase the impact o f that layer on the flow only slightly. In contrast, layers with a low permeability anisotropy ratio (<< 1) cause that layer to have little effect on the flow, even at high aspect ratios. It was found that system asymmetry caused the greatest effect on the flow, even at smaller aspect ratios. A general observation, however, is that smaller system aspect ratios reduce the influence o f layering on the flow (regardless o f system composition) and larger system aspect ratios cause layering to have greater effects.
References Bodes, S. A. and Combamous, M. A., 1973, Natural convection in a sloping porous layer, J. Fluid Mech. 57, 63-79. Dagan, G., 1989, Flow and Transport in Porous Formations, Springer-Verlag, New York. Ludvigsen, A., Palm, E. and McKibbin, R., 1992, Convective momentum and mass transport in porous sloping layers, JournalofGeophys&alResearch, 97, 12315-12325. McKibbin, R. and O'Sullivan, M. J., 1980, Onset of convection in a layered porous medium heated from below, J. Fluid Mech. 96, 375-393. Trew, M., 1991, Modelling convective fluid flow in a system of sloping layered porous media, Engineering Science Project, University of Auckland.