Transport in Porous Media 12: 19-29, 1993. 9 1993 Kluwer Academic Publishers.. Printed in the Netherlands.

19

Natural Convection in a Horizontal Porous Layer with Anisotropic Thermal Diffusivity LEIV STORESLETTEN Department of Mathematics, Agder College, Kristiansand, Norway (Received: 9 August 1991)

Abstract. The present paper is concerned with free convection in a horizontal porous layer with anisotropic thermal diffusivity. It is assumed that the diffusivity has rotational symmetry, with a symmetry

axis making an arbitrary angle against the vertical. The critical Rayleigh number and wave number at marginal stability are calculated and the steady motion occurring at convection onset is examined. It is found that there are two different types of convection cells, depending on whether the longitudinal diffusivity is larger than the transverse diffusivity or not. In the former case, the convection cells are rectangular with vertical lateral walls. In the latter case, however, the lateral cell walls are tilted as well as curved. Key words. Anisotropy, thermal diffusivity, marginal stability, flow pattern, porous layer, convection.

1. I n t r o d u c t i o n Convective flows in p o r o u s media are of interest in m a n y varied problems in geophysics and energy-related systems, like geothermal energy systems, oil reservoir modelling, open-pore insulating systems, and diagenetic processes in sedimentary basins, to name but a few. N a t u r a l convection in anisotropic porous media has attracted the interest of several researchers over the last 15 years. In the middle of the seventies, it was shown that anisotropy in the mechanical and thermal properties effects the marginal stability condition as well as the preferred width of the convection cells (Castinel and C o m b a r n o u s , 1974; Epherre 1975). O n the other hand, Kvernvold and T y v a n d (1979) showed that even a three-dimensional anisotropy does not lead to any new mathematical difficulties or essential new flow patterns at convention onset compared with isotropy. This is true only as long as one of the principal axes of the anisotropic m e d i u m is vertical. This requirement has been maintained in almost all former work in the field, see the review article by M c K i b b i n (1984) as well as M c K i b b i n (1986) and Nilsen and Storesletten (1990). T y v a n d and Storesletten (1991) seem to have been the first to have studied natural convection in an anisotropic p o r o u s m e d i u m where none of the principal axes is vertical. They considered a horizontal p o r o u s layer with anisotropy in the permeability, whereas the thermal diffusivity was isotropic. This was sufficient to achieve qualitatively new flow patterns with tilted plane of m o t i o n or tilted lateral cell walls. In the present work,

20

LEIV STORESLETTEN

we study the analogous problem for a layer with anisotropy in the thermal properties. We consider natural convection in a horizontal fluid-saturated porous layer with anisotropic thermal diffusivity. For simplicity, it is assumed that the diffusivity has rotational symmetry with a symmetry axis making an arbitrary angle against the vertical. The direction of the symmetry axis is denoted longitudinal, which means that the diffusivity is transversely isotropic. Moreover, our analysis is restricted to full isotropy in the permeability. The critical Rayleigh number and wave number at marginal stability are calculated and the steady motion occurring at the onset of convection is examined. It is found that there exist two different types of convection cells (rolls), depending on whether the longitudinal diffusivity is larger than the transverse diffusivity or not. In the former case, the convection ceils are rectangular with vertical lateral walls, whereas the lateral cell walls are tilted as well as curved in the latter case. These results correspond mainly to those found in the analogous problem with anisotropic permeability, studied by Tyvand and Storesletten (1991). On the other hand, there are also a few essential differences, which are discussed in Section 3.

2. Mathematical Formulation We consider free convection in a fluid-saturated porous layer with anisotropic thermal diffusivity. The layer is bounded above and below by two infinite and impermeable heat-conducting horizontal planes. The upper and lower boundaries are separated by a distance h and are at constant temperatures To and To + AT, respectively. Here the characteristic temperature difference AT is positive, which means that the layer is heated from below. It is assumed that the thermal diffusivity has rotational symmetry with a symmetry axis making an arbitrary angle against the vertical. The direction of the symmetry axis is denoted longitudinal, which means that the diffusivity is transversely isotropic. Let ~ci.and xx be the longitudinal and transverse components of the diffusivity tensor D*, i.e. D* = / r

+ lCT(j'j' + k'k')

(1)

where i', j' and k' are unit vectors along the principal axes. A Cartesian frame of reference is chosen, with x- and y-axis at the lower boundary, where the x-axis is aligned along the horizontal projection of i'. The z-axis is directed opposite to gravity. The unit vectors in the x, y, and z directions are denoted i, j, and k. By introducing the anisotropy parameter ,~ =

1~T

--,

KL

(2)

NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER

21

the dimensionless diffusivity tensor D = D*/tcT may be written as D = q-~i'i' + j'j' + k'k' = D l l i i + jj + D13(ik + ki) + D33kk

(3)

Here Dll = t/-1 cos/fl + sin2 fl, D3 3 =

D13

=(t/- 1 __ 1) cos fl sin fl,

COS2 fl ..~ /I - 1 sin / fl

(4)

and fl is the angle between the longitudinal and horizontal direction, i.e. the angle determined by i' and i. The linearized version of the governing equations in dimensionless form may be written (Kvernvold and Tyvand, 1979) v + Vp - R a 0 k = 0,

(5)

V.v = 0,

(6)

00 0 t - w = V-(D-V0)

(7)

Here Darcy law and the Boussinesq approximation have been used and the density is assumed to be a linear function of the temperature. Moreover, v = ui + v] + wk is the velocity, p the pressure, t the time, and 0 is the deviation from the linear temperature distribution corresponding to the motionless conduction state. Ra is the Rayleigh number defined by Ra - K97 ATh - - , /~TV

(8)

where K is the permeability, g the acceleration due to gravity, 7 the thermal expansion coefficient, and v the kinematic viscosity of the saturating fluid. Equation (7) may be written ~0 ~t

~20 ~20 ~20 ~20 w -----D ~ l ~ x2 + -3y - 2 + D33 ~22 4- 2 D l a -C3X - OZ"

(9)

Applying the operator

t?y & ' to Equation (5), adding the component equations, and using the continuity Equation (6), we get the equation V2w = Ra VlZ0,

(10)

22

LEIV STORESLETTEN

where 02

(~2

V2 -- ~ x2 + --@2"

(11)

The operator Ra V1z applied to Equation (9), substituting Equation (10), finally gives the equation V2

Faw a2w uEat - O ~ ax 2

a2w ay 2

a2w O3~-~z 2 -

a2w 7

2O~3aTb~zJ

Fa2w a wl

(12)

= aak~x2 + t~y2]" Impermeable and perfectly heat-conducting boundaries require that w=O=O

at

z=O

and

z=l.

(13)

By using Equation (10), these boundary conditions are expressed by vertical velocity alone: 82w W-

~Z2 - 0

at

z=O

and

z=l.

(14)

3. Marginal Stability and Steady Convection Cells At the onset of convection, the preferred flow cells tend to arrange themselves such that the tangential diffusivity along the streamlines is as small as possible. Let us first consider the case fl = 0. When t1 < 1, we have minimum diffusivity in directions orthogonal to the x-axis. This indicates that the preferred motion is independent of x. When q > 1, it is clear that motion in the y-direction should be avoided, as the diffusivity is maximum in that direction. Thus, it is reasonable to expect convective motion in the (x, y)-plane, independent of y. In the case fi -- 0, these physical arguments suggest that the preferred motion at the onset of convection is independent of x or y depending, respectively, on whether t / < 1 or q > 1. This situation turns out to be valid for arbitrary values of fl, i.e. 0 ~< fl < n. In Appendix A, these hypothesis are numerically confirmed by showing that the Rayleigh number in both cases is a local minimum. Case I: t / < 1. In this case, fl denotes the angle between the x-axis and the direction with maximum thermal diffusivity. Since the solutions are independent of x, it is easy to solve the problem analytically. The curl of Equation (5) gives au

au

ay

az

o,

(15)

NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER

23

which indicates that u is a constant, here given the value zero. It follows that the motion at convection onset is two-dimensional. The convection cells are rectangular with vertical planes of motion. At this point, there is an essential difference between the present case and the analogous case for a porous layer with anisotropic permeability, where the planes of motion were tilted (Tyvand and Storesletten, 1991). This tilt was purely mechanical determined, decoupled from any thermal effects. In the absence of x-dependence the governing Equation (12) reduces to

(02 \

+ 0:VLet

a w7

v33Vj

= Ra j o y2.

(16)

The preferred mode of disturbance which satisfies the boundary conditions is given by w = sin(fez)exp(imy + ~t),

(17)

where m is a wave number, a is the growth rate, and i is the imaginary unit. It is shown in Appendix B that the principle of exchange of stabilities is valid, i.e. a is real and marginal stability is defined by ~r = 0. Substituting solution (17) into Equation (16), we find the Rayleigh number at onset of convection to be Ra = (1 + rc2/mS)(m 2 + D337~s).

(18)

Minimizing Ra with respect to m, we obtain the critical Rayleigh number Rac = ~2(1 + {cos2/~ + q-1 sin s/3}1/2)2.

(19)

The corresponding critical wave number is me = ~ ( c o s 2 j~ + ~ - 1 sin s/3)1/4.

(20)

From Equation (19), it follows that Rao ~ 4res as r / ~ 1. For t/fixed, Rac obtains its minimum 4res for/~ = 0 and its maximum ~s(1 + 1/~/r/) s for /3 = ~r/2. For these two values of/3, the results coincide with those found by Kvernvold and Tyvand (1979). Case II: t / > 1. In this case/~ denotes the angle between the x-axis and the direction with minimal thermal diffusivity. When t / > 1, it turns out that the solutions are independent of y. In order to demonstrate this numerically (see Appendix A), we also include solutions dependent on y, i.e. we consider solutions of the general form w = Z(z)exp(i(kx + my) + ~t),

(21)

where k and m are wave numbers and the amplitude Z has to satisfy the boundary conditions (14), which imply Z(0) = Z"(0) = Z(1) = Z'(1) = 0.

(22)

The solution (21) substituted into Equation (12) leads to a fourth-order linear differential equation for Z with constants coefficients. Its general solution is Z(z) = A1 exp(rlz) + A2 exp(r2z) + As exp(r3z) + A4 exp(r4z),

(23)

24

LEIV STORESLETTEN

where rl, r2, r3, and r4 are roots of the polynomial equation (k 2 + m E - rZ)(a + Dxxk 2 + m 2 - D33 ra - 2D13ikr) = Ra(k 2 + mZ).

(24)

The constants A1, A2, A3, and A4 have to satisfy the boundary conditions (22), which leads to a linear homogeneous system o f algebraic equations. Nontrivial solutions imply that the determinant of the coefficient matrix is zero:

1

1

1

1

rl exp(rl)

exp(rz)

exp(r3)

exp(r4)

rlz exp(rx)

rzexp(r2)

r~exp(r3)

rzexp(r4)

= O.

(25)

As in Case I (q < 1), it follows from Appendix B that the growth rate a = 0 at marginal stability. Thus, in order to find the critical Rayleigh number Rac, we put a = 0 in Equation (24). Given the parameters q, fl and the wave numbers k, m, Equations (24) and (25) represent an eigenvalue problem, the eigenvalues being the Rayleigh numbers Ro < R1 < Rz < ... where the critical Rayleigh number is Rac = min Ro (t/, fl, k, m),

k ~> 0, m >~ 0.

(26)

The eigenvalue problem (24) and (25) (with a = 0) is solved numerically. It turns out that Rac is obtained at m = 0 which means that the steady solutions at convection onset are independent of y. This fact is demonstrated numerically in the Appendix A. In Figure 1 we present marginal stability curves, displaying the Rayleigh number Ra as a function of the wave number k. The curves are computed for the value t/-- 2, for the cases fl = 17.5 ~ 35 ~ 52.5 ~ and 70 ~ respectively. Table I shows the computed values of Ra~ for various values of the anisotropy ratio 11and the angle ft. For moderate values of t/, less than 2 for instance, the critical Rayleigh number depends very weakly on the angle ft. For larger values of ~/there is a stronger dependence, see Table I. For given t/, the critical Rayleigh numbers are equal for the cases fl = 0 ~ and fl = 90 ~ which is known from Kvernvold and Tyvand (1979). Moreover, Rar depends on the angle fl, and for each q there exists an angle tim giving a minimal critical Rayleigh number Rm. Table II shows tim, R,, and the corresponding wave number kin. We observe that tim ~ 45 ~ as r / ~ 1 (r/< 1). It turns out that the computed Tables I and II in the present problem have a close relationship to the corresponding tables in the analogous problem with anisotropic permeability (Tyvand and Storesletten, 1991). In Case II, the tables are in fact identical if t / a n d / / a r e replaced by 1/~ and (~/2 - fl) in the analogous problem. This fact is easily deduced from the governing equations. Putting o- = m = 0 and applying the above-mentioned transformation, the polynomial Equation (24) becomes identical with the corresponding Equation (26) in Tyvand and Storesletten (1991). Concerning the convection cells there exists no such close relationship between the tfeo problems. The plane of motion is vertical and the lateral cell walls are tilted as

25

NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER

Ra

~oo.

:

70 ~

17.5

20.~ o.oo

I

I

I

1

I

(

T

I

i.oo

2.00

3.00

4. oo

s.oo

6.oo

7.00

8.oo

Fig. 1. Marginal stability curves when q = 2 for the cases fl = 17.5~ 35 ~ 52.5 ~ and 70 ~ well as c u r v e d a l s o i n t h e p r e s e n t p r o b l e m . H o w e v e r , t h e tilt is o p p o s i t e d i r e c t e d a n d c a u s e d b y p u r e l y t h e r m a l effects. The computed

s t r e a m l i n e s a r e d i s p l a y e d i n F i g u r e 2 a t fl = 40.1 ~ f o r t h e c a s e s

q = 1.00, 2.00, 4.00, a n d

8.00. A s t r e a m

function has been defined in order to

c o n s t r u c t t h e s e c u r v e s . T h e r e is a c o n s t a n t i n c r e m e n t i n t h e s t r e a m f u n c t i o n b e t w e e n two neighbouring streamlines. T h e l a t e r a l cell w a l l s h a v e n o tilt a t t h e s t a g n a t i o n p o i n t s a t t h e t o p a n d b o t t o m of the layer, whereas the maximum z = 1/2.

tilt is f o u n d i n t h e m i d d l e o f t h e layer, i.e. a t

Table I. The Computed values of Rac for various values of q and fl

fl/q

1.143

1.333

1.60

2.00

2,667

4.00

8.00

0~ 10~ 20~ 30~ 40 ~ 50~ 60~ 70~ 80 ~ 90 ~

36.970 36.968 36.963 36,957 36.952 36.952 36.956 36.962 36.968 36.970

34.367 34.358 34.337 34.312 34,292 34.289 34.303 34.330 34.356 34.367

31.643 31.625 31.577 31.516 31,466 31.452 31.483 31.548 31.615 31.643

28.762 28,730 28,643 28.529 28.427 28.385 28.434 28.561 28.701 28.762

25.658 25.609 25.472 25.282 25.100 25.000 25.053 25.268 25,535 25.658

22.207 22.173 21.940 21.653 21.346 21.131 21.143 21.463 21.955 22.207

18.082 17.994 17.737 17.335 16.848 16.401 16.215 16.571 17.489 18,082

26 Table II.

LEIV STORESLETTEN The c o m p u t e d values of tim, k,,, a n d

R m

for various values of ~/.

r/

1.010

1.143

1.333

1.600

2.000

2.667

4.000

8.000

fl~, k,, R,,

45.1 ~ 3.1416 39.2826

46.0 ~ 3.13 36.952

46.6 ~ 3.14 34.288

48.6 ~ 3.12 31.451

50.1 ~ 3.10 28.385

51.9 ~ 3.07 24.977

54.7 ~ 3.00 21.I01

59.3 ~ 2.83 16.214

@@

% M_

-

-

m

m

J

q = 2.00

r / = 1.00

J = 4.00 Fig. 2.

r / = 8.00

C o m p u t e d streamlines at fl = 40.1 ~ for the cases r / = 1.00, 2.00, 4.00, a n d 8.00.

4. Summary We have considered natural convection in a horizontal porous layer with anisotropy in the thermal diffusivity. It is assumed that the diffusivity has rotational symmetry with a symmetry axis making an angle (90 ~ fl) against the vertical direction. The direction of the symmetry axis is denoted as longitudinal, which means that the diffusivity is transversely isotropic. We have examined the linear stability and the steady flow patterns at the onset of convection. Two different types of convection cells (rolls) were found, both twodimensional: If the longitudinal diffusivity is larger than the transverse (0 < q < 1), the convection cells are rectangular with vertical lateral cell walls like the isotropic case. For the converse case (0 > 1), the plane of motion is vertical whereas the lateral cell walls are tilted as well as curved. The preference for these different flow patterns is explained as a preference for flow directions with as small a tangential diffusivity along the streamlines as possible.

NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER

27

In Case I, 0 < t / < 1, the problem is solved analytically. The critical Rayleigh number is found to be Rac = ~2(1 + {cos2 fl + q-1 sin 2 fl}1/2)2. For given r/, Rac attains its minimum 4~z2 for fl = 0 and its maximum ~z2(1 + 1/X/~) 2 for fl = ~/2. Here fl denotes the angle between the horizontal and longitudinal direction. In Case II, t / > 1, the critical Rayleigh number has to be determined numerically, see Table I. For given t/, the number Rac attains its maximum for fl = 0 and fl -- ~/2. The minimum values Rm depending on fl are calculated in Table II. The computed streamlines are shown in Figure 2.

Appendix A For physical reasons, given in Section 3, it is expected that the solutions at the onset of convection are independent of x (k = 0) or y (m = 0) depending on whether t / < 1 or t / > 1. These conjectures are confirmed numerically by showing that the Rayleigh number for marginal stability attains a local minimum for k -- 0 when t / < 1, and for m = 0 when t / > 1. We have tested the hypothesis for various values of r/and ft. The results are given in the tables on page 28, for q = 0.25, 0.50, 2.00, and 4.00 for the angle f l - - 3 0 ~ From these tables we note that the Rayleigh number at marginal stability is not always a local minimum when the wave number is sufficiently far from its preferred value.

Appendix B: A Proof of the Principle of Exchange of Stabilities (im (a) = 0) We shall now show that the principle of exchange of stabilities is valid for the boundary value problem (5), (6), (7), and (13), i.e. the growth rate a is real and the marginal stability is defined by o- = 0. Multiplying Equations (5) and (7) by v* and Ra 0* (* denotes the complex conjugate), respectively, and adding the equations, we get o Ra 00* = Ra(w0* + w*O) - v*-v - v*. Vp +

+

R / 620 620 620 a~Dll~x2 + ~ y 2 + O 3 3 ~ z 2 +

620 \ 2D13~z)

(B1)

Here ?O/& is replaced by aO, where a is a constant which can be complex. Assuming solutions periodic in x and y, we introduce the following notation ( ( ) ) = 1/Ra

( ) dx dy dz, ,dO

o

o

where (xl - Xo) and (Yl - Yo) are the periods in the x- and y-directions, respectively.

28

LEIV STORESLETTEN

T h e case r / = 0.25 a n d fl = 30~ k=0.00 m=2.7 m=3.00 m=3.25 m=3.50 m=3.75 m=4.00 m=4.25 m=4.50 m=4.75

k=0.01

Ra=57.244859 Ra=55.082069 Ra=53.842696 Ra=53.306997 Ra=53.325932 Ra=53.795531 Ra=54.641471 Ra=55.809482 Ra=57.259188

Ra=57.245007 Ra=55.082274 Ra=53.842936 Ra=53.307257 Ra=53.326205 Ra=53.795811 Ra=54.641755 Ra=55.809768 Ra=57.259475

k=0.10 Ra=57.259714 Ra=55.102637 Ra=53.866668 Ra=53.333021 Ra=53.353197 Ra=53.823536 Ra=54.669902 Ra=55.838140 Ra=57.287948

T h e case r / = 0.50 a n d fl = 30~ k=0.00 m=2.50 m=2.75 m=3.00 m=3.25 m=3.50 m=3.75 m=4.00 m=4.25 m=4.50

k=0.01

Ra=47.938428 Ra=45.869786 Ra=44.735650 Ra=44.296813 Ra=44.396313 Ra=44.927696 Ra=45.816695 Ra=47.010223 Ra=48.469517

Ra=47.938380 Ra=45.869820 Ra=44.735734 Ra=44.296927 Ra=44.396448 Ra=44.927843 Ra=45.816851 Ra=47.010386 Ra=48.469683

k=0.10 Ra=47.933689 Ra=45.873194 Ra=44.744004 Ra=44.308267 Ra=44.409762 Ra=44.942454 Ra=45.832326 Ra=47.026442 Ra=48.486133

T h e case q = 2 a n d fl = 30~ m=0.00 k=2.50 k=2.75 k=3.00 k=3.25 k=3.50 k=3.75 k=4.00 k=4.25 k=4.50

m=0.01

Ra=31.681117 Ra=30.085635 Ra=29.130747 Ra=28.649077 Ra=28.530817 Ra=28.701502 Ra=29.109190 Ra=29.716757 Ra=30.497097

T h e case r / = 4 a n d fl

=

Ra=31.675987 Ra=30.086267 Ra=29.134885 Ra=28.655420 Ra=28.538583 Ra=28.710205 Ra=29.118520 Ra=29.726511 Ra=30.507141

30~

m=0.00 k=3.00 k=3.25 k=3.50 k=3.75 k=4.00 k=4.25 k=4.50 k=4.75 k=5.00

Ra=31.681065 Ra=30.085641 Ra=29.130788 Ra=28.649141 Ra=28.530895 Ra=28.701589 Ra=29.109283 Ra=29.716855 Ra=30.497198

m=0.10

Ra=23.135779 Ra=22.364496 Ra=21.904812 Ra=21.687596 Ra=21.664629 Ra=21.801453 Ra=22.072908 Ra=22.460281 Ra=22.949416

m=0.01 Ra=23.135858 Ra=22.364590 Ra=21.904916 Ra=21.687706 Ra=21.664743 Ra=21.801569 Ra=22.073025 Ra=22.460398 Ra=22.949533

m=0 .i0 Ra=23.143594 Ra=22.373917 Ra=21.915205 Ra=21.698580 Ra=21.675970 Ra=21.813001 Ra=22.084570 Ra=22.471996 Ra=22.961145

NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER

29

Integration of Equation (B1) over the fluid volume (Xo, x l)x(yo, Yl )x(0, 1), application of Equations (6) and (13), and partial integration, finally produces the equation ~r(O0*) = (wO* + w * O ) -

/

(1/Ra)(v*-v) -

00 00*

/ 0 0 00" - D13

00 00*

00 0 0 * \

00" 0 0 \

V; +

(B2)

Since all terms at the right-hand side are real, it follows that im(a) = O. Acknowledgements The author is grateful to Dr Peder A. Tyvand, Agricultural University of Norway, for having directed his attention to the problem investigated in this paper. He also wishes to thank his colleague Dr Asvald Lima for helpful instruction in connection with the computing and the graphics, and Professor Enok Palm for helpful discussions. References Castinel, G. and Combarnous, M., 1974, crit~re d'apparition de la convection naturelle dans une couche poreuse anisotrope horizontale, C.R. Acad. Sci. Ser. B 287, 701-704. Epherre, J. F., 1975, Critbre d'apparition de la convection naturelle dans une couche poreuse anisotrope, Rev. Therrnique 168, 949-950. Kvernvold, O. and Tyvand, P. A., 1979, Nonlinear thermal convection in anisotropic porous media, J. Fluid Mech. 90, 609-624. McKibbin, R., 1984, Thermal convection in layered and anisotropic porous media: A review, in R. A. Wooding and I. White (eds), Convective Flows in Porous Media, DSIR, Wellington, New Zealand, pp. t13-127. McKibbin, R., 1986, Thermal convection in a porous layer: Effects of anisotropy and surface boundary conditions, Transport in Porous Media 1, 271-292. Nilsen, T. and Storesletten, L., 1990, An analytical study on natural convection in isotropic and anisotropic porous channels, Trans. ASME C: J. Heat Transfer 112, 396~401. Tyvand, P. A. and Storesletten, L., 1991, Onset of convection in an anisotropic porous medium with oblique principal axes, J. Fluid Mech. 226, 371-382.