- Email: [email protected]
Natural convection in a thin horizontal porous annulus
In!. J Hear Mess Transfer. Printed in Great Britain
Vol. 30. No. 4. pp. 729 739, 1987
Natural convection
Co17 9310187 f3.00+0.00 ,p 1987 Pergamon Journals Ltd.
in a thin horizontal annulus
MIHIR SENt and KENNETH Sibley School
of Mechanical
and Aerospace
Engineering,
porous
E. TORRANCE Cornell
University,
Ithaca,
NY 14853, U.S.A.
(Received 7 April 1986 and in final form 6 August 1986) Abstract-The possible modes of time-dependent natural convection in a horizontal annulus of finite length are considered. The annulus is filled with porous material and the annular thickness is assumed small in comparison with the mean radius. All boundaries are impermeable and adiabatic; heating is through a circumferentially distributed volumetric heat source. The governing equations reduce to a set of two non-linear ordinary differential equations. Steady non-linear oscillations exist for the special case of infinite Rayleigh number and symmetric heating about the vertical. For lower Rayleigh numbers, damped oscillations are obtained, the degree of damping increasing with the inclination of the line of symmetry and with decreasing Rayleigh number. Multiple stable steady states are obtained for small inclinations. Chaotic motions do not develop for non-inertial Darcy flows.
INTRODUCTION
sional numerical study by Facas and Farouk [l l] showed the convergence of transient convection to a steady-state pattern with no significant intermediate oscillatory states. In a recent extension of earlier work [7], Robillard et al. [12,13], examined a horizontal porous annulus with a constant temperature inner surface and with a circumferentially-varying, sinusoidal temperature on the outer surface. Flow symmetry was not assumed. For low Rayleigh numbers a perturbation analysis about the conduction solution was used; at higher Rayleigh numbers numerical methods were employed. At Rayleigh numbers above a critical, oscillatory flows were calculated. When a time-averaging process was applied to these flows, they displayed characteristics similar to previous steady-state flows. For much higher Rayleigh numbers, no periodicity in the timedependent flow was detected and the system exhibited qualitative ‘chaotic’ behavior. Comparison of refs. [ 10,l l] with ref. [ 121 suggests that the two-dimensional oscillatory flows for relatively low Rayleigh numbers in ref. [12] might be a consequence of the assumed circumferential variation of temperature on the outer surface of the annulus. On the other hand, flow symmetry was not assumed in ref. [ 123 but was assumed in refs. [ 10,111. Due to the complexity of the governing partial differential equations a full reconciliation of the foregoing differences may not be possible without extensive numerical calculations. Nevertheless, it is possible under certain conditions to study the onset of oscillations and possible chaotic behavior in a simpler manner. To do this, it is necessary to simplify the governing equations to a set of ordinary differential equations. This has been done for other natural convection systems (e.g. ref. [14]). Such a simplified analysis would indicate the possible modes of steadystate and time-dependent behavior.
convection in annular porous materials is important in many applications. In particular, when porous materials are used for thermal insulation on pipes, the heat losses are directly related to any convective flows that may develop within the insulation. Furthermore, such convective flows strongly depend on heating conditions and system geometry, and may be time dependent. This paper examines the possible modes of time-dependent convection that can develop in a thin-walled porous annulus. Recent reviews of natural convection in porous materials are available [1,2]. Several studies have considered horizontal annular regions, under steadystate conditions, with the two cylindrical surfaces kept at different temperatures. The effects of cylinder eccentricity [3,4], inclination of the axis of the annulus [S], thermal radiation [6], and internal heat generation [7] have been examined, along with other effects [8,9]. Most of these studies have used numerical methods to solve the governing partial differential equations. A considerable amount of information on steady-state heat transfer has been accumulated for a range of Rayleigh numbers and radius ratios. Some work has also been done on time-dependent convection and stability when the cylindrical surfaces of an annulus are kept at different temperatures. Caltagirone [lo] examined stability by the Galerkin method and found a Rayleigh number-dependent convection pattern. On increasing the Rayleigh number an initial conduction regime gave way to two-dimensional convection. With further increases of the Rayleigh number, three-dimensional flows and associated fluctuations appeared. The two-dimen-
NATURAL
ton leave from National Autcinoma
Facultad de Ingenieria, Universidad de Mkxico, M6xic0, D.F., Mexico. 729
730
M. SENand K. E. TORRANCE
NOMENCLATURE a, b
c C g K L Pp P-, P+
:I r, 8, z R Ra
t T AT u, 0, w v, w x> Y
non-dimensional heat source constants defined in equations (l&f and Wd) ratio of heat capacity of porous matrix to that of fluid integration constant in equation (17) acceleration due to gravity permeability of porous medium length of annulus pressure period of oscillation critical points volumetric heat source characteristic magnitude of heat source defined after equation (11) ~ylind~~a1 coordinates defined in Fig. 1 mean radius of annulus Darcy-Rayleigh number defined in equation (log) time temperature characteristic temperature difference defined in equation (IOf) velocity components in cylindrical coordinates Fourier coefficients defined in equations (37) and (38) non-dimensional Fourier coefficients defined in equations (10a) and (lob)
9 ;i
/I... ,::>
Z ‘.
_’)R
I+-----LB FIG. 1. Geometry of porous annular material of mean radius R and length L. The present study introduces a thin annulus approximation to make the time-dependent problem analytically tractable. By contrast, most prior studies have assumed an annulus radius ratio of about two. Non-uniform heating is simulated through a given volumetric heat source which is a function of the @ coordinate in Fig. 1. This heat source distribution is quite general, and can be thought of as arising from external heat addition or removal at the inner or outer circumferential boundaries of the thin annulus.
Ymin,y,,,
minimum and maximum values of Y during a periodic oscillation.
Greek symbols a thermal diffusivity coefficient of volumetric expansion B discriminant defined in equation (22) Y 6 magnitude of inertia term in equation (28) dynamic viscosity of fluid P V Liapunov function in equation (26) fluid density P u eigenvalue 7 non-dimensional time defined in equation (1Oe). Subscripts n, m 0
indices of Fourier coefficients index of first term in Fourier expansion.
Superscripts c cosine coefficient in Fourier expansion S sine coefficient in Fourier expansion perturbation from critical point * reference value. Symbols 0
critical point or steady-state value average over a period.
The only constraint is that the heat source distribution must integrate to zero around the circumference of the annulus in order to allow steady-state solutions. The present study simplifies the governing equations and examines the temporal and steady-state behavior of the solutions. The objective is to classify the non-Iinear modes of time-de~ndent motion in a thin porous annulus. The following sections outline the governing equations under the thin annulus approximation, and the solution. Results are presented for high and intermediate Rayleigh numbers, and for symmetric and asymmetric heating. The paper closes with a discussion of inertial effects.
THIN
ANNULUS
GOVERNING
EQUATIONS
We can begin our analysis by considering the finitelength thr~-dimensional annulus shown in Fig. 1. The radial, circumferential and axial coordinates are r, S, and z; the corresponding velocity components are u, v, and w. The annulus is thin in comparison with its mean radius and consists of a fluid-saturated
Natural convection in a thin horizontal porous annulus
porous medium. All boundary surfaces are adiabatic and impermeable. Heating is through a prescribed rate of volumetric heat addition. As the thickness of the annulus gets smaller, so will the radial fluid velocity in the porous medium because of the impermeability of the boundaries. On neglecting the radial velocity it can be shown that there is no axial flow in the annulus. In the interests of brevity the details of this proof are relegated to the appendix. Thus, under the thin annulus approximation, there is a considerable simplification of the governing equations and the problem reduces to a circumferential flow in the annulus. With the convective motion confined to the circumferential direction, the temperature is a function only of 0. An obvious victim of this approximation will be any recirculation zones within the radial thickness of the porous material. Only if higher order expansions are made in terms of the thickness of the annulus can such cellular activity be detected. The mode of heating is also an important factor, since, for certain conditions in a thin annulus, only the trivial radial conductive solution will be obtained. The internal volumetric heat source distribution used here will avoid this difficulty. It is important to point out that the present approach differs in detail though not in basic analytical technique from solutions obtained from series expansions in terms of the Rayleigh number [3, 4, 121. The advantage here is that the results are valid for large Rayleigh numbers and are thus susceptible to diverse time-dependent behaviors. The derivation of the governing equations parallels the procedure for toroidal thermosyphons [lS]. The thin annulus porous medium equations are given in the appendix and only the planar versions will be used here. If u is the circumferential velocity component, by mass conservation (29a) it must be a function of time alone. The Darcy-Oberbeck-Boussinesq equations, equations (29b) and (29d), are u=-
K
convection in porous media [l, 12, 163. Integrating equation (1) from 0 = 0 to 21r, the pressure term is eliminated T(0, t) cos 0 d0.
-PI{1
- /?(T-
T*)}gcose - ;$j
T=
f [T:(t)sin(nB) “=I
q”’ = f
+
7yt)COS(ne)l
[q: sin (no) + q;
00s fne)].
“=I
Substituting
(4)
(5)
in equation (3) we obtain ”
-
Kf lgpy 2/l
”
The energy equation, equation (2), gives for the sine and cosine modes, respectively
Equation (6) can be substituted in equations (7) and (8) to obtain an infinite set of autonomous ordinary differential equations. Of these the first two with n = 1 decouple from the rest. The nondimensional form of the first two equations is dx -=Y2-b-& dr dY ds=a-xY-Ka
Y
where
(1)
x = T”,/AT
(loa)
y = T’,IAT
(lob)
a = @i/Q
b = -G/Q 7 =
where K is the permeability of the porous medium, p the dynamic viscosity, p the density, /i’the coefficient of thermal expansion of the fluid, p the pressure, T the temperature, R the mean radius of the annulus, c the ratio of the heat capacity of the porous matrix to that of the fluid, q”’ the volumetric rate of heat generation divided by the heat capacity of the fluid, and a the effective thermal diffusivity. The starred variables are at reference conditions. It should be pointed out that the non-inertial form of Darcy’s law has been used, which is suitable for fine grained materials as indicated in ref. [lo]. This form has been used often in analyses of time-dependent natural
(3)
We expand T and q”’ in terms of Fourier series in the following form
1 c~+~g=4”‘(8)+-g (2)
P[
731
AT=
Kp+g/lATt/2cpR
&PQIKp*dO
Ra = Kp*gjlATR/2pa.
(104 ww uw
uw Wf3)
The dependent variables x and y represent the nondimensional sine and cosine amplitudes of the first Fourier mode. In addition, since the azimuthal velocity component, u is proportional to the amplitude of the first cosine mode (see equation (6)), y is related to u and dimensional parameters by y = 2pv/Kp*gpAT.
(11)
Only the first Fourier coefficients of the heat source
132
M. SENand K. E. TORRANCE
distribution are of importance in the determination of the velocity. A characteristic heat input is defined by Q = ,/((&)* + (qcl)*). Thus, a and b are of order unity and satisfy the relation a* + bZ = 1, so that only one of them need be specified. The relative magnitudes of a and b depend on the inclination of the line of symmetry of the heating, with a = 0 corresponding to heating which is symmetrical about a vertical line. The tangent of the angle of inclination of the line of symmetry is a/b. In equation (log) Ra is the Darcy-Rayleigh number. Note that the governing equations, equations (9a) and (9b), have been nondimensionalized such that the geometrical effect of the volumetric heating is through the parameters a and b. These parameters do not vanish even if the heating tends to zero. The Rayleigh number employs the mean radius of the annulus as the length scale. Thus, in the limit of an infinite Rayleigh number, conduction becomes negligible with respect to convection, and Ra disappears from equations (9a) and (9b). However, the volumetric heating parameters a and b remain and control the structure of the convective motion. Thus, the infinite Rayleigh number limit can provide important information on the preferred structure of time-dependent convective flows in the annulus. In the present study, the flow characteristics at Ra = 100 or higher are very similar to those at an infinite Ra. Thus, to illustrate the structure of the flows, we will consider Ra values of infinity and 10. The problem has been reduced to the solution of two ordinary differential equations, equations (9a) and (9b), under the thin annulus approximation. These equations will be studied using a variety of techniques including linear stability analysis, exact integration under special conditions, Liapunov method, and numerical analysis. Whenever numerical solutions are reported they are obtained with a fourth order Runge-Kutta scheme with a time step of 0.01.
below and cooling from above and will be assumed in all that follows. A linear stability analysis around the critical points of the non-linear solution gives the eigenvalues (T= ;[-2
f &’
with x and y proportional for P+ and P- are thus
- Sj*)]
(14)
to exp(ar). The eigenvalues
(15a)
n=$-$~/(~-8b)]
and o=f[$fJ(;-8b)]
(W
respectively. Instability results from a positive real part of the eigenvalue. Thus P+ (Pm) is unstable if a < 0 (a z=-0) and stable otherwise. The special case of symmetrical heating about a vertical line is the present analog to the problem studied by Robillard et al. [123 and is considered next. Symmetrical
heating, a = 0
From the paragraph following equation (1 l), when a = 0, it follows that b = 1. By equations (15a) and (15b), neutral stability is implied for a = 0, with linear oscillations around the critical points P*:(x,y)
= (0, * 1)
of radian frequency 42. This special non-linear problem has an exact solution in phase space, given by x2 + y2 - lny* = C.
(17)
Constant C is determined from initial conditions and is always larger than unity. The value C = 1 corresponds to critical point initial conditions. The solution is invariant under a y + -y transformation and is thus symmetrical about the x-axis. INFINITE RAYLEIGH NUMBER LIMIT Some solutions of equation (17), in the form of closed orbits in the x, y phase plane, are shown in In the limit of an infinitely high Rayleigh number the governing equations, equations (9a) and (9b), Fig. 2(a) for C = 2, 5 and 10. These correspond to initial values of (1, l), (2,l) and (3,l) in the upper half simplify to of the plane and to (1, - l), (2, - 1) and (3, - 1) in the dx lower half. The time dependence of x is the same for dz=y2-b (124 both positive and negative initial conditions and is shown in Fig. 2(b) from numerical calculations. The dy (12b) corresponding y variation is shown in Fig. 2(c) for z = a - xy. positive initial conditions. For negative initial conThe steady-state solutions are the critical points of ditions, the y variation is found by replacing the y the set of equations. These are denoted by P+ and ordinate in Fig. 2(c) with -y. The larger amplitude oscillations tend to form peaks indicating high freP-, with values quency content. The period P is generally a function P+:(x,y) = (a/,,‘b,Jb) (134 of amplitude for non-linear oscillations. This is shown in Fig. 3 where P is found to increase as the amplitude P - : (2, j) = - (u/,/b, ,,/‘b) (W gets larger. Very small oscillations about the critical point exist as C + 1, for which the period is found to and exist only for b > 0. This represents heating from
Natural convection in a thin horizontal porous annulus
733
FIG. 4. Maximum (y,,,) and minimum (y,J values of y(r) during an oscillation for infinite Rayleigh number symmetric heating. Broken lines are approximations.
-21 -3
x
the present case all orbits are entirely in the upper or lower halves of the phase plane. In other words the fluid velocity does not reverse while going through a series of oscillations. This differs from ref. [ 123, where even the steady-state oscillations had zero mean circumferential velocity. The maximum and minimum of the absolute values of Y during the oscillations are given by the solutions to the transcendental equation
4 3 2 I x
0 -I -2
Y2 - lnY* = C.
-3 -d
-4o-
I5
IO
5
20
r
-41
0
5
IO
I5
20
T
FIG. 2. (a) Closed trajectories in x-y phase plane for infinite Rayleigh number symmetrical heating, C= 2, 5, and 10; (b) X(T)transients; (C)J(T)transients for positive initial values (transients for negative initial values are mirror symmetric).
01 j
I
c
1
’
5
0
’
1
’C IO
FIG. 3. Period of oscillation, P, for infinite Rayleigh number symmetric heating and different values of C. A value of C = 1 corresponds to critical point initial conditions.
be 71,/2 as predicted by eigenvalues (15a) and (15b) of the linear theory. These oscillations have a superficial resemblance to the numerical results of ref. [12], although there are significant differences in detail. For example, in
(18)
These values are denoted by Y,,,.. and Y,,,~,and are shown in Fig. 4. For large C, they can be approximated by Y,,, = JC
(19a)
Ymin= exp ( - C)
(lob)
as shown by broken lines in the figure. As C increases, the lowest velocity during each oscillation (as indicated by Ymi,,)decreases so that the flow almost comes to a complete stop before accelerating again. The information presented in ref. [12] also suggests that for that problem the time-averaged characteristics were similar to those for the steady-state solution. In the present case we can easily show from the governing equations that (x) = 0
(2Oa)
(Y2> = 1
(2Ob)
(XY> = 0
(2W
where averages over a period are represented by ( ). On comparing equation (20b) with equations (13a) and (13b) we find that the r.m.s. velocity for periodic oscillations is exactly equal to what would occur under steady state. Asymmetrical
heating, a # 0 On either side of a = 0, one or the other of the critical points becomes asymptotically stable (from equations (Isa) and (15b)). From numerical integrations, almost the entire x-y plane seems to be the basin of attraction of these critical points, ruling out finite amplitude oscillations. For very small a, slowly damped oscillations around the stable critical point can be observed which for large a are highly damped. This behavior is shown in Fig. 5 for an initial condition
134
M. SEN and K. E.
4,
-1
3
TORRANCE 4
I
r
3
2 I
0
I
(b)
-I
-2
-3 -4
-4 1 0
4
3
IU T
I
15
20
FIG. 5. (a) Trajectory in x-y phase plane for infinite Rayleigh number asymmetric heating with n = 0.1 and initial condition (0, -0.5). The solution migrates from an unstable initial state to a stable steadystate. (b) y (5) transient.
(0, -0.5) and for a = 0.1, which represents an inclination of 5.7” in the line of symmetry of the heating with respect to the vertical. Figure 5(a) presents the x-y phase plane in which the solution curve is seen to first spiral out from the unstable critical point Pand then spiral in towards the stable critical point P+. According to the eigenvalues, equations (15a) and (15b), the rates of divergence and convergence are equal. The change in direction of the velocity is also observed in Fig. S(b) which shows the time variation of y. FINITE
RAYLEIGH
NUMBER
(
-&b _a
>
y-&=0
FIG. 6. Parameter values for finite Rayleigh number convection which lead to three real solutions (hatched region) and to one real solution (unhatched region). Symmetric heating corresponds to a = 0, and asymmetric heating to a # 0.
Symmetrical heating, a = 0
CONVECTION
For finite Rayleigh numbers, the complete equations, equations (9a) and (9b), must be considered. However, the results should asymptotically match those for infinite Ra of the preceding section. The critical points for finite Ra are given by solutions of the cubic equation
y3+
l/ Ro2
Once again we consider, the special case of a = 0 and b = 1 for which the roots of equation (21a) can be directly determined by factorization. The solutions are 1
(21a)
x== -G’ 1
“=Y-G.
@lb)
Equation (21a) has either one real root or three depending on the value of the discriminant ,=(-&--b)l/27+a’/4Ra2.
(22)
The hatched and unhatched areas in Fig. 6 correspond to parameter values which give three solutions and one solution, respectively. A linear stability analysis can be carried out in a similar fashion to the preceding section to obtain the eigenvalues
.=~[-(i+~)‘J(il-8~z)].
1
i=
p=
(24a)
-J(I
-Ra,
--&)
j = 0.
(24b) (24~)
These are shown in Figs. 7(a) and (b), and are denoted by I, II and III, respectively. It is seen that solutions I and II are only valid for Ra 2 1 while the trivial solution III is valid for all Ra. Substituting these in the eigenvalue expression, equation (23), the linear stability theory indicates that I and II are always stable while III is only stable for Ra < 1. For small enough Ra there is only the conducting solution with zero velocity. That this is asymptotically stable for all perturbations can be shown in the following manner. The governing equations can be written in local form as
(23)
As compared to equation (14), the finite Rayleigh number analysis is seen to lead to a term which stabilizes small deviations from the critical point.
ji=
x=-Ra’
z =(Y')~d - y’ = d7
x’/Ra
Ray’ - x’y’ - y’/Ra
(25b)
Natural convection in a thin horizontal porous annulus
-.2
II
1
(01 2-
I
TV
I,
*
I
,
-I
-
-2
-
5
10
II
, Ra
x
III
FIG. 7. Steady-state solutions for finite Rayleigh number convection with symmetric heating. (a) f and (b)j vs Ra. The
three solution branches are denoted by I, II, and III. Dotted lines represent unstable solutions.
where the primes denote perturbations from the steady solution. We consider the positive definite function v = (x32 + (y’)2.
(26)
From this we have dv avdx avdy dt=axds+ay;iT = -2(~‘)~/Ra - 2Ra
(27)
which is negative definite for Ra < 1, so that under this condition v is a Liapunov function. This means that solutions starting from any initial condition tend to the conducting one. On increasing Ra from zero, the initially stable conduction solution gives way to two convective solutions at Ra = 1. By comparison with the infinite Ra approximation, the extra steady-state solution III for finite Ra is unstable and is thus not important. Furthermore, on the asymptotically stable branches I and II, the convergence is of the order of exp ( - l/Ra). Thus, initial conditions away from the critical point will lead to damped oscillations which could be mistaken for periodic motion in numerical calculations at large Ra. Figures 8(a) and (b) show an example of such a slowly damped motion for Ra = 10 and for positive and negative initial values of y. Once again the velocity does not change direction during the oscillations. The constant amplitude oscillations of Fig. 2 can be compared to these. Asymmetric heating, a # 0
This is the most general case of equations (9a) and (9b). Figures 9(a) and (b) show the evolution from the initial condition (0, -0.5) for Ra = 10, a = 0.1. In the phase plane graph the solution is seen to spiral away very slowly from the unstable critical point before it
735
changes direction and spirals in fairly rapidly to the stable critical point. Figure 9 can be compared to Fig. 5 where Ra was considered infinite. Though the qualitative trends are the same, the speed of divergence from and convergence to critical points is not. Two stable steady-state velocities are sometimes possible. This is illustrated in Figs. lo(a) and (b) for u = 0.1, which show X and jj for different values of Ra. The three solution branches are indicated by I, II and III with the unstable parts shown as dotted. Branch I exists for all Ra and is always stable. Branch II exists above a certain critical value Ramin but is entirely unstable. Branch III also exists above this critical value and is stable only in an interval Ramin < Ra < Ramax between the crosses. Figure 11 shows Ramin and Ramax for different values of a. As a -+ 0, Ramaxbecomes infinite. It can also be seen that multiple stable steady states are possible only for a < 0.335, which corresponds to an inclination of the symmetry line of the heating of less than 20” from the vertical. Whenever two stable steady states are possible, the long time tendency of the time-dependent solution towards either of these will depend on the initial condition. We consider as an example Ra = 2 and points a = 0.1, for which the two critical (- 0.62, -0.83) and (- 0.39,0.89) are linearly stable. In Fig. 12 all initial conditions marked with a filled circle lead to the first critical point indicated with a larger filled circle. This defines the finite basin of attraction of this critical point attractor. Initial conditions corresponding to the crosses tend to the second critical point, marked with a bold cross, and this basin of attraction is found to be infinite. According to the eigenvalues of the linear stability theory, equation (23), the convergence is 60% faster for the second critical point than for the first. Furthermore it is not evidently the initial velocity alone that determines the final steady-state condition. The two basins of attraction are separated by an unstable limit cycle and it is not easy to determine what an initial condition on this limit cycle would do. In terms of the foregoing terminology, for the symmetric heating case a = 0, both basins of attraction are infinite, being the upper half of the phase plane for the upper critical point and the lower half for the lower critical point. In such a case the initial velocity alone determines the final steady state to which the system would converge.
INERTIAL
EFFECTS AND POSSIBLE SOLUTIONS
CHAOTIC
The non-linear ordinary differential equations, equations (9a) and (9b), which govern this problem have solutions which can be represented as trajectories in a two-dimensional phase plane. The kind of attractors possible for bounded long time solutions under such circumstances are only steady states and closed
M.
736
and K. E. TORRANCE
SEN
4
4 3
I
3
2
I
2
I .--__------__
Y 0 . I
_
1
__
‘I
-1
j
_\
’
;
-2
-__-___.-. -I 0 X
’
-2 -3
._ _I_.__ -3
,I
‘,
I,,’ /’
-3
\,
-I
)
-2
4-
Lb)
I---
_j
‘\\
-4
_
0
(01
.I
I
2
.._ ___,_ 4 3
i
-4
20
FIG. 8. (a)Trajectories in x-y phase plane for finite Rayleigh number symmetric heating, Ra = 10. Initial conditions are (3,1) for the upper trajectory and (3, - 1) for the lower trajectory. (b)y(s) transients.
I -I
7
-I
-2
-2
-3 -4 -4L-
--.Ibl
-3 -3
-2
-I
0 X
I
2
3
4
-4
0
IO
20
30
40 T
50
60
70
80
FIG. 9. (a) Trajectory in x-y phase plane for finite Rayleigh number asymmetric heating, Ra = 10 and a = 0.1. The solution migrates from an unstable initial state to a stable steady state. (b)y(r) transient.
-a 04
FIG. 11. Maximum (Ram.3 and minimum (Ra,i,,) values of the Rayleigh number for the conditionally stable solution branch III in Fig. 10.
Cb)
FIG. 10. Steady-state solutions for finite Rayleigh number convection with asymmetric heating, a = 0.1. (a)x and (b)j vs Ra. The solution branches are denoted by I, II, and III. Dotted lines represent unstable solutions.
orbits. Strange attractors with chaotic solutions are not possible as would be the case for dynamical systems in three-dimensional phase spaces, which is what occurs in natural convection in closed tubes with toroidal geometry [14,15]. The reason that the present problem has a two-dimensional and not a three-dimensional representation in phase space is that no time derivatives of the velocity appear in Darcy’s law. Consideration of fluid inertia leads to
the inclusion of other terms [7, lo] also, but as they do not increase the dimension of the phase space they will be ignored in the present discussion. Given the importance of the foregoing point in relation to possible chaos in the system, it is relevant to ask what effect the inclusion of a time derivative in Darcy’s law would have, even if it were very small in magnitude. The governing equations would then be of the form dx -=zy-b dt
dy
(28)
iG=a-rx
sg=y-z where for simplicity
we consider
the non-conducting
Natural convection in a thin horizontal porous annulus T
.......
:::::l.......::::: :::::::::::::::::;::::::::f:::;;;;;;;;; ..............................
:::.::::::::::::::::::::::::::::::::::: ...+*...*.....*.+*. .................... ....................................... ......................................... :::::::::::.........~....., :::::::::::::I::::::::::::;:::::::::::: _....................... ......-_. ._~~...-~...-....__~_.~-.-.~~~.~~..__...! ............. ................. ............ ... .+..*.....................* ..* ................... ..‘.......................~........* ................................................ ............ . .... ..*., ..P ...................... ................... *... ..l...........* .................... ................................... ............ ...................... ..:::: .. .. . ................................... i:............ .. *. .. - . ..*.+................* ..* :: :: ............. ............ ............. ..* ............ ............ . . .. . “.............R :: :: :::;:::;::::t:::;::::::::::::~ :: ............*..f . . ............................ ............ -2 -I -3 0 I 2 I I
t
.
I
.
.+..
*
.
::
FIG. 12. Basins of attraction of the two stable critical points for Ra = 2 and a = 0.1 are indicated by filled circles and crosses in the x-y phase plane. The corresponding critical points are shown with bold symbols.
case. Inertial effects are included through 6 and it is easy to see that for 6 = 0, the non-inertial equations, equations (12a) and (12b), are recovered. The critical points of equations (12a), (12b) and (28) are identical and are given by equations (13a) and (13b). The critical point P+ is stable for a > ~56~” and P- for a < -8~~“. There is a narrow range of values of a near zero of width proportional to 6, where both critical points are unstable. In this region either closed orbits or much more complicated chaotic behaviors are possible. A description of the different possibilities is given in ref. [ 151 for 6 = 1. However, 6 is very small here and the range of values of a for which chaotic solutions could occur is also very small and of little practical importance in the present context. Moreover, it must be remembered that the non-inertial form of Darcy’s law was used in ref. [ 123 and thus the chaotic motion found there must be due to some other reason, either the boundary conditions used or two dimensionality. This merits further study.
CONCLUSIONS
The governing equations for time-dependent natural convection arising from an arbitrary heat source distribution in a thin annulus have been reduced to the solution of two non-linear ordinary differential equations. The approximation does not involve any restriction on the Rayleigh number though effects due to possible recirculation zones within the annular thickness have been suppressed. For infinite Rayleigh numbers and symmetrical heating an exact solution exists. This solution, combined with a linear stability analysis, provides information on the general behavior of the system. For asymmetric heating or finite Rayleigh numbers, damped velocity oscillations can be observed leading to a steady-state motion. The damping factors can be qualitatively predicted from the eigenvalues of linear stability theory. Under cer-
737
tain conditions two stable steady states are possible, and either can be reached depending on the initial conditions. It can also be demonstrated that chaotic behavior is not possible under the thin annulus approximation unless fluid inertia effects are included. Two other observations can be made with respect to possible numerical solutions of the governing partial differential equations under these and similar circumstances. Any results implying steady oscillations should be carefully verified in order to distinguish between strictly periodic phenomena and slowly damped oscillations. Also, in time-dependent calculations, different initial conditions can lead to different steady states.
Acknowledfement-We acknowledge the oartial suooort of the Na&nal Science Foundation _ under _-Grant MEA 8401489.
REFERENCES 1. P. Cheng, Heat transfer in geothermal systems, Adu. Heat Transfer 14, l-105 (1978). 2. M. A. Combamous and S. A. Bories, Hydrothermal convection in saturated porous media, Ado. Hydrosci.
10,231-307 (1975). 3. H. H. Bau, Low Rayleigh number thermal convection in a saturated porous medium bounded by two horizontal eccentric cylinders, ASME J. Heat Transfer 106, 166175 (1984). 4. H. H. Bau, Thermal convection in a horizontal, eccentric
annulus containing a saturated porous medium-an extended perturbation analysis, Int. J. Heat Mass Transfer 27, 2277-2287 (1984). 5. Y. Takata, K. Fukuda,
S. Hasegawa, K. Iwashige, H. Shimomura and K. Sanokawa, Three-dimensional natural convection in a porous medium enclosed with concentric inclined cylinders, Proc. 7th Int. Heat Transfer Conf., Munich, Vol. 2, pp. 351-356 (1982). 6. R. Echigo, S. Hasegawa, S. Tottori, H. Shimomura and Y. Okamoto, An analysis on the radiative and free convective heat transfer in a horizontal annulus with permeable insulator, Proc. 6th Int. Heat Transfer Conf., Toronto, Vol. 3, pp. 385-390 (1978). 7. P. Vasseur, T. H. Nguyen, L. Robillard and V. K. T. Thi, Natural convection between horizontal concentric cylinders filled with a porous layer with internal heat generation, Int. J. Heat Mass Transfer 27, 337-349 (1984). 8. V. A. Brailovskaya,
G. B. Petrazhitskii and V. I. Polezhaer, Natural convection and heat transfer in porous interlayers between horizontal coaxial cylinders, J. Appl. Mech. Tech. Phys. 19, 781-786 (1978). 9. P. J. Burns and C. L. Tien, Natural convection in porous media bounded by concentric spheres and horizontal cylinders, Int. J. Heat Mass Transfer 22,929-939 (1979). 10. J.-P. Caltagirone, Thermoconvective instabilities in a porous medium bounded by two concentric horizontal cylinders, J. Fluid Mech. 76, 337-362 (1976). 11. G. N. Facas and B. Farouk, Transient and steady-state natural convection in a porous medium between two concentric cylinders, ASME J. Hear Transfer 105, 660663 (1983). 12. L. Robillard, T. H. Nguyen, M. G. Satish and P. Vasseur,
Oscillating flow in an annular porous layer heated from
M. SEN and K. E. TORRANCE
738
a++-, -_=I)
below, Heat Transfer in Porous Media and Particulate Flows (edited by L. S. Yao, M. M. Chen, C. E. Hickox,
P. G. Simpkins, L. C. Chow, M. Kaviany, P. Cheng and L. R. Davis), HTD-Vol. 46, pp. 41-47. The American Society of Mechanical Engineers, New York (1985). 13. T. H. Nguyen, M. G. Satish, L. Robillard and P. Vasseur, Free convection in a non-uniformly heated annular porous layer, ASME Paper No. 85-WA/HT-8. The American Society of Mechanical Engineers, New York (1985). 14. J. A. Yorke and E. D. Yorke, Chaotic behavior and fluid dynamics, Hydrodynamic Instabilities and the Transition to Turbulence (edited by H. L. Swinney and J. L. Gollub), 2nd edn, pp. 77-95. Springer, Berlin (1985). 15. M. Sen, E. Ramos and C. Treviiio, The toroidal thermosyphon with known heat flux, Int. J. Heat Mass Transfer 28, 219-233
az
-
’
n;+2
= 0,
2+as=o (31a-c)
uo = Kp*gBT;/&, II: = -nKpyRp
vf, = nKp’JRp, w.
=
---
K ah p az'
w:,=_--2, K aps
wc” =
KapC,
F aZ ~2
+ -!-x(-nv:TE, 2R n
+ nuC,T:) + w,,$
=a2
ar nv c._!-__?Te,fw at R
(33a-c)
P az
(1985).
16. T. Masuoka, M. Tashiro and T. Katsuhara, Transient natural convection in a porous medium within a closed sphere, Proc. JSME-ASME Thermal Engng Joint Conf., Honolulu, pp. 243-249 (1983).
(32a-c)
ar ---‘+w:o O a2
a*T,
(34a)
aT
aZ
=q:-!$~+a!$
(34b)
APPENDIX
Here we show that the present study is applicable not only for a planar geometry but also for a three-dimensional annulus with a horizontal axis, as shown in Fig. 1. This is due to the fact that under the thin annulus approximation, the assumption of impermeable and adiabatic boundaries leads to zero flow and heat conduction in the axial direction. For a thin annulus with negligible radial velocity, we have the following governing equations for flow in porous media. Mass conservation gives
f!+g=O
- B(T- T*)JgcosO
!!%T’+w ’
az+wcs
O az
” az
(34c) Equations (31a)-(31c) can be integrated in the z-direction to give wg = 0
(294
where w is the axial component of the fluid velocity. The Darcy-Oberbeck-Boussinesq equations in the circumferential and axial coordinate directions are respectively
u = -$p’(l
aTE,
‘dt+R
+k$]
(29b)
dz = 0
{2} = ${z}sin(Znmz/L)
k ap
+ {g} The Fourier components expanded as
The energy equation is [I]
{z} The four dependent variables II,w, T and p are Zn-periodic in 0 so that they can be expanded in Fourier series of the form
(36)
where impermeability at z = 0 and L has been used. Equation (36) indicates that Fourier series expansions of vi and vi in z may also be used. These can be written as
w=-lz’
I$ = & + C [& sin (no) + 4’. cos (no)] n
(35)
COS(2nnu,Q].
(37)
of the axial velocity can be
= ~{~‘“}sin(2nrnz/~).
(38)
Expansions (37) and (38) can be substituted in equations (31b) and (31~) to obtain yes nn = y’s nm= 0
(39)
(30)
where 4 represents each one of the variables. The new dependent variables & and +‘, are functions of z and t alone. The summations are from n = 1 to infinity, as will be all the other summations. The heat source can also be expanded as in equation (5). Three component equations can be obtained for each one of equations (30), one for the integral itself, and then two for the integral of sin(m0) and cos(m0) times the equation. These equations are
In addition equations (31b), (31c), (32b), (32c), (33b) and (33~) can be used to obtain w’ I
R2azw: 2 a2
0
’
wc - $2
”
= 0.
(41a,b)
Substituting equation (38) in equations (41a) and (41b) and ,.^~ using equation (4U),we get
Natural
convection
in a thin horizontal
porous annulus
739
Thus, under the thin annulus approximation, ‘j$&
+ @$){~}sin(2nnz/L) “In
= 0
component
(43)
It can be shown in a similar manner from the energy equations (34a)-(34c) that, as a consequence of the adiabatic condition we have assumed at z = 0 and L, the temperature is independent ofz. Thus, only the circumferential component of the velocity field exists and the temperature varies only
(44)
with 9.
from which
This gives
CONVECTION
NATURELLE
the velocity
(42)
DANS UN ANNEAU
in the z-direction
is zero.
POREUX HORIZONTAL
considtre les modes possibles de coinvection naturelle, dependant du temps, dans un espace finie. L’anneau est rempli de matiere poreuse et son epaisseur est supposee petite en comparaison du rayon moyen. Toutes les frontieres sont impermeables et adiabatiques; le chauffage est fourni par une source volumetrique de chaleur distribuie circonferentiellement. Les equations se reduisent a un systcme de deux equations differentielles non lintaires. Des oscillations stables non lineaires existent pour le cas special d’un nombre de Rayleigh intini et un chauffage symttrique par rapport a la verticale. Pour des nombres de Rayleigh plus faibles, on obtient des oscillations amorties dont le degre d’amortissement augmente avec I’inclinaison de la ligne de symttrie et quand le nombre de Rayleigh diminue. Des etats stables permanents sont obtenus pour des petites inclinaisons. Des mouvements chaotiques ne se developpent pas dans les Ccoulements de Darcy non inertiels. R&urn&On
annulaire de longueur
NATtiRLICHE
KONVEKTION IN EINEM DtiNNEN, PORBSEN RINGSPALT
WAAGERECHTEN
Zusammenfassung-Die
moglichen Zustlnde der zeitabhangigen natiirlichen Konvektion in einem horizontalen Ringspalt endlicher Lange werden betrachtet. Der Ringspalt ist mit porosem Material gefiillt, und die Ringspaltbreite wird als klein im Vergleich zum mittleren Radius angenommen. Alle Grenzflachen sind dicht und adiabat ; die Beheizung erfolgt durch eine iiber den Umfang verteilte volumetrische Warmequelle. Die Bilanzgleichungen reduzieren sich auf zwei nichtlineare gewohnliche Differentialgleichungen. Im Spezialfall bei unendlicher Rayleigh-Zahl und symmetrischer Beheizung iiber die Vertikale treten stetige nichtlineare Schwingungen auf. Bei niedrigen Rayleigh-Zahlen treten gediimpfte Schwingungen auf, wobei der Grad der Dampfung mit der Neigung der Symmetrieachse und mit abnehmender Rayleigh-Zahl ansteigt. Bei kleineren Neigungen ergeben sich mehrere stabile Zustande. Bei nicht tragheitsbedingten Darcy-Stromungen entwickeln sich keine chaotischen Bewegungen.
ECTECTBEHHAg
KOHBEKHHR B TOHKOM TOPM30HTAJIbHOM KOJIbHEBOM KAHAJIE
AnHoTaluln-PaccMaTpesalorcn
IIOPHCTOM
B03MOEHbIe ~WiMbl eCTeCTBeHHOi? KOHBeKUWli B TOPR30HTanbHOM KonbueBoM Kaaane KoHeqHoii nnuHbI. KOJIbIIeBOfi KaHan 3anonHeH nopHcTbIM MaTepHanoM, ero TonmliHa cqHTaeTca Manoii no cpaeseaew, co cpenHHM paneycoM. Bee rpaHHubI HenpomiUaeMble H anHa6aTHrecKHe; HarpeB OCyIUeCTBJTaeTCa06aeMHbtMn BCTOYHWKaMH TelTJIa, paCllOnO)l(eHHbIMWII0 OKPy~HOCT1(.O~Pe~eJDI‘OUIHe ypaBHeHHK CBeneHbl KCACTeMenByX HenHHetiHbIX06blKHOBeHHbIX nH+j,ePeHUHanbHbIX ypaBHeHk& B Cny'iae 6eCKOHevHOrO WiCJIa kWIe!4 N CHMMeTpHHHOrO HarpeBa HaBnIonamrca ycToiiHHBbIe HenHUeiiHbIe KOne6aHHa. J&m 6onee HW~KWXqricen Psnea nonyqeno 3aryxanne KOne6aHHti,npWieM CTeneHb 3aTy3aHHa BO3paCTaeT C HaKJIOHOMnHHUA CHMMeTpHH H C yMeHb"IeHHeM wfcna Psnen. ,&a ManbIx yrnos HaKnoHa Ha6nmnanHCb MHoromcneHHbIe ycToi+isBbIeCOCTOIIHHR.B CJIy’Iae 6e3bIHepUHOHHbIXTe’IeHHfi flapC!H XaOTH'teCKHeLIBH)KeHAII He pa3BHaaEOTCR.
Vol. 30. No. 4. pp. 729 739, 1987
Natural convection
Co17 9310187 f3.00+0.00 ,p 1987 Pergamon Journals Ltd.
in a thin horizontal annulus
MIHIR SENt and KENNETH Sibley School
of Mechanical
and Aerospace
Engineering,
porous
E. TORRANCE Cornell
University,
Ithaca,
NY 14853, U.S.A.
(Received 7 April 1986 and in final form 6 August 1986) Abstract-The possible modes of time-dependent natural convection in a horizontal annulus of finite length are considered. The annulus is filled with porous material and the annular thickness is assumed small in comparison with the mean radius. All boundaries are impermeable and adiabatic; heating is through a circumferentially distributed volumetric heat source. The governing equations reduce to a set of two non-linear ordinary differential equations. Steady non-linear oscillations exist for the special case of infinite Rayleigh number and symmetric heating about the vertical. For lower Rayleigh numbers, damped oscillations are obtained, the degree of damping increasing with the inclination of the line of symmetry and with decreasing Rayleigh number. Multiple stable steady states are obtained for small inclinations. Chaotic motions do not develop for non-inertial Darcy flows.
INTRODUCTION
sional numerical study by Facas and Farouk [l l] showed the convergence of transient convection to a steady-state pattern with no significant intermediate oscillatory states. In a recent extension of earlier work [7], Robillard et al. [12,13], examined a horizontal porous annulus with a constant temperature inner surface and with a circumferentially-varying, sinusoidal temperature on the outer surface. Flow symmetry was not assumed. For low Rayleigh numbers a perturbation analysis about the conduction solution was used; at higher Rayleigh numbers numerical methods were employed. At Rayleigh numbers above a critical, oscillatory flows were calculated. When a time-averaging process was applied to these flows, they displayed characteristics similar to previous steady-state flows. For much higher Rayleigh numbers, no periodicity in the timedependent flow was detected and the system exhibited qualitative ‘chaotic’ behavior. Comparison of refs. [ 10,l l] with ref. [ 121 suggests that the two-dimensional oscillatory flows for relatively low Rayleigh numbers in ref. [12] might be a consequence of the assumed circumferential variation of temperature on the outer surface of the annulus. On the other hand, flow symmetry was not assumed in ref. [ 123 but was assumed in refs. [ 10,111. Due to the complexity of the governing partial differential equations a full reconciliation of the foregoing differences may not be possible without extensive numerical calculations. Nevertheless, it is possible under certain conditions to study the onset of oscillations and possible chaotic behavior in a simpler manner. To do this, it is necessary to simplify the governing equations to a set of ordinary differential equations. This has been done for other natural convection systems (e.g. ref. [14]). Such a simplified analysis would indicate the possible modes of steadystate and time-dependent behavior.
convection in annular porous materials is important in many applications. In particular, when porous materials are used for thermal insulation on pipes, the heat losses are directly related to any convective flows that may develop within the insulation. Furthermore, such convective flows strongly depend on heating conditions and system geometry, and may be time dependent. This paper examines the possible modes of time-dependent convection that can develop in a thin-walled porous annulus. Recent reviews of natural convection in porous materials are available [1,2]. Several studies have considered horizontal annular regions, under steadystate conditions, with the two cylindrical surfaces kept at different temperatures. The effects of cylinder eccentricity [3,4], inclination of the axis of the annulus [S], thermal radiation [6], and internal heat generation [7] have been examined, along with other effects [8,9]. Most of these studies have used numerical methods to solve the governing partial differential equations. A considerable amount of information on steady-state heat transfer has been accumulated for a range of Rayleigh numbers and radius ratios. Some work has also been done on time-dependent convection and stability when the cylindrical surfaces of an annulus are kept at different temperatures. Caltagirone [lo] examined stability by the Galerkin method and found a Rayleigh number-dependent convection pattern. On increasing the Rayleigh number an initial conduction regime gave way to two-dimensional convection. With further increases of the Rayleigh number, three-dimensional flows and associated fluctuations appeared. The two-dimen-
NATURAL
ton leave from National Autcinoma
Facultad de Ingenieria, Universidad de Mkxico, M6xic0, D.F., Mexico. 729
730
M. SENand K. E. TORRANCE
NOMENCLATURE a, b
c C g K L Pp P-, P+
:I r, 8, z R Ra
t T AT u, 0, w v, w x> Y
non-dimensional heat source constants defined in equations (l&f and Wd) ratio of heat capacity of porous matrix to that of fluid integration constant in equation (17) acceleration due to gravity permeability of porous medium length of annulus pressure period of oscillation critical points volumetric heat source characteristic magnitude of heat source defined after equation (11) ~ylind~~a1 coordinates defined in Fig. 1 mean radius of annulus Darcy-Rayleigh number defined in equation (log) time temperature characteristic temperature difference defined in equation (IOf) velocity components in cylindrical coordinates Fourier coefficients defined in equations (37) and (38) non-dimensional Fourier coefficients defined in equations (10a) and (lob)
9 ;i
/I... ,::>
Z ‘.
_’)R
I+-----LB FIG. 1. Geometry of porous annular material of mean radius R and length L. The present study introduces a thin annulus approximation to make the time-dependent problem analytically tractable. By contrast, most prior studies have assumed an annulus radius ratio of about two. Non-uniform heating is simulated through a given volumetric heat source which is a function of the @ coordinate in Fig. 1. This heat source distribution is quite general, and can be thought of as arising from external heat addition or removal at the inner or outer circumferential boundaries of the thin annulus.
Ymin,y,,,
minimum and maximum values of Y during a periodic oscillation.
Greek symbols a thermal diffusivity coefficient of volumetric expansion B discriminant defined in equation (22) Y 6 magnitude of inertia term in equation (28) dynamic viscosity of fluid P V Liapunov function in equation (26) fluid density P u eigenvalue 7 non-dimensional time defined in equation (1Oe). Subscripts n, m 0
indices of Fourier coefficients index of first term in Fourier expansion.
Superscripts c cosine coefficient in Fourier expansion S sine coefficient in Fourier expansion perturbation from critical point * reference value. Symbols 0
critical point or steady-state value average over a period.
The only constraint is that the heat source distribution must integrate to zero around the circumference of the annulus in order to allow steady-state solutions. The present study simplifies the governing equations and examines the temporal and steady-state behavior of the solutions. The objective is to classify the non-Iinear modes of time-de~ndent motion in a thin porous annulus. The following sections outline the governing equations under the thin annulus approximation, and the solution. Results are presented for high and intermediate Rayleigh numbers, and for symmetric and asymmetric heating. The paper closes with a discussion of inertial effects.
THIN
ANNULUS
GOVERNING
EQUATIONS
We can begin our analysis by considering the finitelength thr~-dimensional annulus shown in Fig. 1. The radial, circumferential and axial coordinates are r, S, and z; the corresponding velocity components are u, v, and w. The annulus is thin in comparison with its mean radius and consists of a fluid-saturated
Natural convection in a thin horizontal porous annulus
porous medium. All boundary surfaces are adiabatic and impermeable. Heating is through a prescribed rate of volumetric heat addition. As the thickness of the annulus gets smaller, so will the radial fluid velocity in the porous medium because of the impermeability of the boundaries. On neglecting the radial velocity it can be shown that there is no axial flow in the annulus. In the interests of brevity the details of this proof are relegated to the appendix. Thus, under the thin annulus approximation, there is a considerable simplification of the governing equations and the problem reduces to a circumferential flow in the annulus. With the convective motion confined to the circumferential direction, the temperature is a function only of 0. An obvious victim of this approximation will be any recirculation zones within the radial thickness of the porous material. Only if higher order expansions are made in terms of the thickness of the annulus can such cellular activity be detected. The mode of heating is also an important factor, since, for certain conditions in a thin annulus, only the trivial radial conductive solution will be obtained. The internal volumetric heat source distribution used here will avoid this difficulty. It is important to point out that the present approach differs in detail though not in basic analytical technique from solutions obtained from series expansions in terms of the Rayleigh number [3, 4, 121. The advantage here is that the results are valid for large Rayleigh numbers and are thus susceptible to diverse time-dependent behaviors. The derivation of the governing equations parallels the procedure for toroidal thermosyphons [lS]. The thin annulus porous medium equations are given in the appendix and only the planar versions will be used here. If u is the circumferential velocity component, by mass conservation (29a) it must be a function of time alone. The Darcy-Oberbeck-Boussinesq equations, equations (29b) and (29d), are u=-
K
convection in porous media [l, 12, 163. Integrating equation (1) from 0 = 0 to 21r, the pressure term is eliminated T(0, t) cos 0 d0.
-PI{1
- /?(T-
T*)}gcose - ;$j
T=
f [T:(t)sin(nB) “=I
q”’ = f
+
7yt)COS(ne)l
[q: sin (no) + q;
00s fne)].
“=I
Substituting
(4)
(5)
in equation (3) we obtain ”
-
Kf lgpy 2/l
”
The energy equation, equation (2), gives for the sine and cosine modes, respectively
Equation (6) can be substituted in equations (7) and (8) to obtain an infinite set of autonomous ordinary differential equations. Of these the first two with n = 1 decouple from the rest. The nondimensional form of the first two equations is dx -=Y2-b-& dr dY ds=a-xY-Ka
Y
where
(1)
x = T”,/AT
(loa)
y = T’,IAT
(lob)
a = @i/Q
b = -G/Q 7 =
where K is the permeability of the porous medium, p the dynamic viscosity, p the density, /i’the coefficient of thermal expansion of the fluid, p the pressure, T the temperature, R the mean radius of the annulus, c the ratio of the heat capacity of the porous matrix to that of the fluid, q”’ the volumetric rate of heat generation divided by the heat capacity of the fluid, and a the effective thermal diffusivity. The starred variables are at reference conditions. It should be pointed out that the non-inertial form of Darcy’s law has been used, which is suitable for fine grained materials as indicated in ref. [lo]. This form has been used often in analyses of time-dependent natural
(3)
We expand T and q”’ in terms of Fourier series in the following form
1 c~+~g=4”‘(8)+-g (2)
P[
731
AT=
Kp+g/lATt/2cpR
&PQIKp*dO
Ra = Kp*gjlATR/2pa.
(104 ww uw
uw Wf3)
The dependent variables x and y represent the nondimensional sine and cosine amplitudes of the first Fourier mode. In addition, since the azimuthal velocity component, u is proportional to the amplitude of the first cosine mode (see equation (6)), y is related to u and dimensional parameters by y = 2pv/Kp*gpAT.
(11)
Only the first Fourier coefficients of the heat source
132
M. SENand K. E. TORRANCE
distribution are of importance in the determination of the velocity. A characteristic heat input is defined by Q = ,/((&)* + (qcl)*). Thus, a and b are of order unity and satisfy the relation a* + bZ = 1, so that only one of them need be specified. The relative magnitudes of a and b depend on the inclination of the line of symmetry of the heating, with a = 0 corresponding to heating which is symmetrical about a vertical line. The tangent of the angle of inclination of the line of symmetry is a/b. In equation (log) Ra is the Darcy-Rayleigh number. Note that the governing equations, equations (9a) and (9b), have been nondimensionalized such that the geometrical effect of the volumetric heating is through the parameters a and b. These parameters do not vanish even if the heating tends to zero. The Rayleigh number employs the mean radius of the annulus as the length scale. Thus, in the limit of an infinite Rayleigh number, conduction becomes negligible with respect to convection, and Ra disappears from equations (9a) and (9b). However, the volumetric heating parameters a and b remain and control the structure of the convective motion. Thus, the infinite Rayleigh number limit can provide important information on the preferred structure of time-dependent convective flows in the annulus. In the present study, the flow characteristics at Ra = 100 or higher are very similar to those at an infinite Ra. Thus, to illustrate the structure of the flows, we will consider Ra values of infinity and 10. The problem has been reduced to the solution of two ordinary differential equations, equations (9a) and (9b), under the thin annulus approximation. These equations will be studied using a variety of techniques including linear stability analysis, exact integration under special conditions, Liapunov method, and numerical analysis. Whenever numerical solutions are reported they are obtained with a fourth order Runge-Kutta scheme with a time step of 0.01.
below and cooling from above and will be assumed in all that follows. A linear stability analysis around the critical points of the non-linear solution gives the eigenvalues (T= ;[-2
f &’
with x and y proportional for P+ and P- are thus
- Sj*)]
(14)
to exp(ar). The eigenvalues
(15a)
n=$-$~/(~-8b)]
and o=f[$fJ(;-8b)]
(W
respectively. Instability results from a positive real part of the eigenvalue. Thus P+ (Pm) is unstable if a < 0 (a z=-0) and stable otherwise. The special case of symmetrical heating about a vertical line is the present analog to the problem studied by Robillard et al. [123 and is considered next. Symmetrical
heating, a = 0
From the paragraph following equation (1 l), when a = 0, it follows that b = 1. By equations (15a) and (15b), neutral stability is implied for a = 0, with linear oscillations around the critical points P*:(x,y)
= (0, * 1)
of radian frequency 42. This special non-linear problem has an exact solution in phase space, given by x2 + y2 - lny* = C.
(17)
Constant C is determined from initial conditions and is always larger than unity. The value C = 1 corresponds to critical point initial conditions. The solution is invariant under a y + -y transformation and is thus symmetrical about the x-axis. INFINITE RAYLEIGH NUMBER LIMIT Some solutions of equation (17), in the form of closed orbits in the x, y phase plane, are shown in In the limit of an infinitely high Rayleigh number the governing equations, equations (9a) and (9b), Fig. 2(a) for C = 2, 5 and 10. These correspond to initial values of (1, l), (2,l) and (3,l) in the upper half simplify to of the plane and to (1, - l), (2, - 1) and (3, - 1) in the dx lower half. The time dependence of x is the same for dz=y2-b (124 both positive and negative initial conditions and is shown in Fig. 2(b) from numerical calculations. The dy (12b) corresponding y variation is shown in Fig. 2(c) for z = a - xy. positive initial conditions. For negative initial conThe steady-state solutions are the critical points of ditions, the y variation is found by replacing the y the set of equations. These are denoted by P+ and ordinate in Fig. 2(c) with -y. The larger amplitude oscillations tend to form peaks indicating high freP-, with values quency content. The period P is generally a function P+:(x,y) = (a/,,‘b,Jb) (134 of amplitude for non-linear oscillations. This is shown in Fig. 3 where P is found to increase as the amplitude P - : (2, j) = - (u/,/b, ,,/‘b) (W gets larger. Very small oscillations about the critical point exist as C + 1, for which the period is found to and exist only for b > 0. This represents heating from
Natural convection in a thin horizontal porous annulus
733
FIG. 4. Maximum (y,,,) and minimum (y,J values of y(r) during an oscillation for infinite Rayleigh number symmetric heating. Broken lines are approximations.
-21 -3
x
the present case all orbits are entirely in the upper or lower halves of the phase plane. In other words the fluid velocity does not reverse while going through a series of oscillations. This differs from ref. [ 123, where even the steady-state oscillations had zero mean circumferential velocity. The maximum and minimum of the absolute values of Y during the oscillations are given by the solutions to the transcendental equation
4 3 2 I x
0 -I -2
Y2 - lnY* = C.
-3 -d
-4o-
I5
IO
5
20
r
-41
0
5
IO
I5
20
T
FIG. 2. (a) Closed trajectories in x-y phase plane for infinite Rayleigh number symmetrical heating, C= 2, 5, and 10; (b) X(T)transients; (C)J(T)transients for positive initial values (transients for negative initial values are mirror symmetric).
01 j
I
c
1
’
5
0
’
1
’C IO
FIG. 3. Period of oscillation, P, for infinite Rayleigh number symmetric heating and different values of C. A value of C = 1 corresponds to critical point initial conditions.
be 71,/2 as predicted by eigenvalues (15a) and (15b) of the linear theory. These oscillations have a superficial resemblance to the numerical results of ref. [12], although there are significant differences in detail. For example, in
(18)
These values are denoted by Y,,,.. and Y,,,~,and are shown in Fig. 4. For large C, they can be approximated by Y,,, = JC
(19a)
Ymin= exp ( - C)
(lob)
as shown by broken lines in the figure. As C increases, the lowest velocity during each oscillation (as indicated by Ymi,,)decreases so that the flow almost comes to a complete stop before accelerating again. The information presented in ref. [12] also suggests that for that problem the time-averaged characteristics were similar to those for the steady-state solution. In the present case we can easily show from the governing equations that (x) = 0
(2Oa)
(Y2> = 1
(2Ob)
(XY> = 0
(2W
where averages over a period are represented by ( ). On comparing equation (20b) with equations (13a) and (13b) we find that the r.m.s. velocity for periodic oscillations is exactly equal to what would occur under steady state. Asymmetrical
heating, a # 0 On either side of a = 0, one or the other of the critical points becomes asymptotically stable (from equations (Isa) and (15b)). From numerical integrations, almost the entire x-y plane seems to be the basin of attraction of these critical points, ruling out finite amplitude oscillations. For very small a, slowly damped oscillations around the stable critical point can be observed which for large a are highly damped. This behavior is shown in Fig. 5 for an initial condition
134
M. SEN and K. E.
4,
-1
3
TORRANCE 4
I
r
3
2 I
0
I
(b)
-I
-2
-3 -4
-4 1 0
4
3
IU T
I
15
20
FIG. 5. (a) Trajectory in x-y phase plane for infinite Rayleigh number asymmetric heating with n = 0.1 and initial condition (0, -0.5). The solution migrates from an unstable initial state to a stable steadystate. (b) y (5) transient.
(0, -0.5) and for a = 0.1, which represents an inclination of 5.7” in the line of symmetry of the heating with respect to the vertical. Figure 5(a) presents the x-y phase plane in which the solution curve is seen to first spiral out from the unstable critical point Pand then spiral in towards the stable critical point P+. According to the eigenvalues, equations (15a) and (15b), the rates of divergence and convergence are equal. The change in direction of the velocity is also observed in Fig. S(b) which shows the time variation of y. FINITE
RAYLEIGH
NUMBER
(
-&b _a
>
y-&=0
FIG. 6. Parameter values for finite Rayleigh number convection which lead to three real solutions (hatched region) and to one real solution (unhatched region). Symmetric heating corresponds to a = 0, and asymmetric heating to a # 0.
Symmetrical heating, a = 0
CONVECTION
For finite Rayleigh numbers, the complete equations, equations (9a) and (9b), must be considered. However, the results should asymptotically match those for infinite Ra of the preceding section. The critical points for finite Ra are given by solutions of the cubic equation
y3+
l/ Ro2
Once again we consider, the special case of a = 0 and b = 1 for which the roots of equation (21a) can be directly determined by factorization. The solutions are 1
(21a)
x== -G’ 1
“=Y-G.
@lb)
Equation (21a) has either one real root or three depending on the value of the discriminant ,=(-&--b)l/27+a’/4Ra2.
(22)
The hatched and unhatched areas in Fig. 6 correspond to parameter values which give three solutions and one solution, respectively. A linear stability analysis can be carried out in a similar fashion to the preceding section to obtain the eigenvalues
.=~[-(i+~)‘J(il-8~z)].
1
i=
p=
(24a)
-J(I
-Ra,
--&)
j = 0.
(24b) (24~)
These are shown in Figs. 7(a) and (b), and are denoted by I, II and III, respectively. It is seen that solutions I and II are only valid for Ra 2 1 while the trivial solution III is valid for all Ra. Substituting these in the eigenvalue expression, equation (23), the linear stability theory indicates that I and II are always stable while III is only stable for Ra < 1. For small enough Ra there is only the conducting solution with zero velocity. That this is asymptotically stable for all perturbations can be shown in the following manner. The governing equations can be written in local form as
(23)
As compared to equation (14), the finite Rayleigh number analysis is seen to lead to a term which stabilizes small deviations from the critical point.
ji=
x=-Ra’
z =(Y')~d - y’ = d7
x’/Ra
Ray’ - x’y’ - y’/Ra
(25b)
Natural convection in a thin horizontal porous annulus
-.2
II
1
(01 2-
I
TV
I,
*
I
,
-I
-
-2
-
5
10
II
, Ra
x
III
FIG. 7. Steady-state solutions for finite Rayleigh number convection with symmetric heating. (a) f and (b)j vs Ra. The
three solution branches are denoted by I, II, and III. Dotted lines represent unstable solutions.
where the primes denote perturbations from the steady solution. We consider the positive definite function v = (x32 + (y’)2.
(26)
From this we have dv avdx avdy dt=axds+ay;iT = -2(~‘)~/Ra - 2Ra
(27)
which is negative definite for Ra < 1, so that under this condition v is a Liapunov function. This means that solutions starting from any initial condition tend to the conducting one. On increasing Ra from zero, the initially stable conduction solution gives way to two convective solutions at Ra = 1. By comparison with the infinite Ra approximation, the extra steady-state solution III for finite Ra is unstable and is thus not important. Furthermore, on the asymptotically stable branches I and II, the convergence is of the order of exp ( - l/Ra). Thus, initial conditions away from the critical point will lead to damped oscillations which could be mistaken for periodic motion in numerical calculations at large Ra. Figures 8(a) and (b) show an example of such a slowly damped motion for Ra = 10 and for positive and negative initial values of y. Once again the velocity does not change direction during the oscillations. The constant amplitude oscillations of Fig. 2 can be compared to these. Asymmetric heating, a # 0
This is the most general case of equations (9a) and (9b). Figures 9(a) and (b) show the evolution from the initial condition (0, -0.5) for Ra = 10, a = 0.1. In the phase plane graph the solution is seen to spiral away very slowly from the unstable critical point before it
735
changes direction and spirals in fairly rapidly to the stable critical point. Figure 9 can be compared to Fig. 5 where Ra was considered infinite. Though the qualitative trends are the same, the speed of divergence from and convergence to critical points is not. Two stable steady-state velocities are sometimes possible. This is illustrated in Figs. lo(a) and (b) for u = 0.1, which show X and jj for different values of Ra. The three solution branches are indicated by I, II and III with the unstable parts shown as dotted. Branch I exists for all Ra and is always stable. Branch II exists above a certain critical value Ramin but is entirely unstable. Branch III also exists above this critical value and is stable only in an interval Ramin < Ra < Ramax between the crosses. Figure 11 shows Ramin and Ramax for different values of a. As a -+ 0, Ramaxbecomes infinite. It can also be seen that multiple stable steady states are possible only for a < 0.335, which corresponds to an inclination of the symmetry line of the heating of less than 20” from the vertical. Whenever two stable steady states are possible, the long time tendency of the time-dependent solution towards either of these will depend on the initial condition. We consider as an example Ra = 2 and points a = 0.1, for which the two critical (- 0.62, -0.83) and (- 0.39,0.89) are linearly stable. In Fig. 12 all initial conditions marked with a filled circle lead to the first critical point indicated with a larger filled circle. This defines the finite basin of attraction of this critical point attractor. Initial conditions corresponding to the crosses tend to the second critical point, marked with a bold cross, and this basin of attraction is found to be infinite. According to the eigenvalues of the linear stability theory, equation (23), the convergence is 60% faster for the second critical point than for the first. Furthermore it is not evidently the initial velocity alone that determines the final steady-state condition. The two basins of attraction are separated by an unstable limit cycle and it is not easy to determine what an initial condition on this limit cycle would do. In terms of the foregoing terminology, for the symmetric heating case a = 0, both basins of attraction are infinite, being the upper half of the phase plane for the upper critical point and the lower half for the lower critical point. In such a case the initial velocity alone determines the final steady state to which the system would converge.
INERTIAL
EFFECTS AND POSSIBLE SOLUTIONS
CHAOTIC
The non-linear ordinary differential equations, equations (9a) and (9b), which govern this problem have solutions which can be represented as trajectories in a two-dimensional phase plane. The kind of attractors possible for bounded long time solutions under such circumstances are only steady states and closed
M.
736
and K. E. TORRANCE
SEN
4
4 3
I
3
2
I
2
I .--__------__
Y 0 . I
_
1
__
‘I
-1
j
_\
’
;
-2
-__-___.-. -I 0 X
’
-2 -3
._ _I_.__ -3
,I
‘,
I,,’ /’
-3
\,
-I
)
-2
4-
Lb)
I---
_j
‘\\
-4
_
0
(01
.I
I
2
.._ ___,_ 4 3
i
-4
20
FIG. 8. (a)Trajectories in x-y phase plane for finite Rayleigh number symmetric heating, Ra = 10. Initial conditions are (3,1) for the upper trajectory and (3, - 1) for the lower trajectory. (b)y(s) transients.
I -I
7
-I
-2
-2
-3 -4 -4L-
--.Ibl
-3 -3
-2
-I
0 X
I
2
3
4
-4
0
IO
20
30
40 T
50
60
70
80
FIG. 9. (a) Trajectory in x-y phase plane for finite Rayleigh number asymmetric heating, Ra = 10 and a = 0.1. The solution migrates from an unstable initial state to a stable steady state. (b)y(r) transient.
-a 04
FIG. 11. Maximum (Ram.3 and minimum (Ra,i,,) values of the Rayleigh number for the conditionally stable solution branch III in Fig. 10.
Cb)
FIG. 10. Steady-state solutions for finite Rayleigh number convection with asymmetric heating, a = 0.1. (a)x and (b)j vs Ra. The solution branches are denoted by I, II, and III. Dotted lines represent unstable solutions.
orbits. Strange attractors with chaotic solutions are not possible as would be the case for dynamical systems in three-dimensional phase spaces, which is what occurs in natural convection in closed tubes with toroidal geometry [14,15]. The reason that the present problem has a two-dimensional and not a three-dimensional representation in phase space is that no time derivatives of the velocity appear in Darcy’s law. Consideration of fluid inertia leads to
the inclusion of other terms [7, lo] also, but as they do not increase the dimension of the phase space they will be ignored in the present discussion. Given the importance of the foregoing point in relation to possible chaos in the system, it is relevant to ask what effect the inclusion of a time derivative in Darcy’s law would have, even if it were very small in magnitude. The governing equations would then be of the form dx -=zy-b dt
dy
(28)
iG=a-rx
sg=y-z where for simplicity
we consider
the non-conducting
Natural convection in a thin horizontal porous annulus T
.......
:::::l.......::::: :::::::::::::::::;::::::::f:::;;;;;;;;; ..............................
:::.::::::::::::::::::::::::::::::::::: ...+*...*.....*.+*. .................... ....................................... ......................................... :::::::::::.........~....., :::::::::::::I::::::::::::;:::::::::::: _....................... ......-_. ._~~...-~...-....__~_.~-.-.~~~.~~..__...! ............. ................. ............ ... .+..*.....................* ..* ................... ..‘.......................~........* ................................................ ............ . .... ..*., ..P ...................... ................... *... ..l...........* .................... ................................... ............ ...................... ..:::: .. .. . ................................... i:............ .. *. .. - . ..*.+................* ..* :: :: ............. ............ ............. ..* ............ ............ . . .. . “.............R :: :: :::;:::;::::t:::;::::::::::::~ :: ............*..f . . ............................ ............ -2 -I -3 0 I 2 I I
t
.
I
.
.+..
*
.
::
FIG. 12. Basins of attraction of the two stable critical points for Ra = 2 and a = 0.1 are indicated by filled circles and crosses in the x-y phase plane. The corresponding critical points are shown with bold symbols.
case. Inertial effects are included through 6 and it is easy to see that for 6 = 0, the non-inertial equations, equations (12a) and (12b), are recovered. The critical points of equations (12a), (12b) and (28) are identical and are given by equations (13a) and (13b). The critical point P+ is stable for a > ~56~” and P- for a < -8~~“. There is a narrow range of values of a near zero of width proportional to 6, where both critical points are unstable. In this region either closed orbits or much more complicated chaotic behaviors are possible. A description of the different possibilities is given in ref. [ 151 for 6 = 1. However, 6 is very small here and the range of values of a for which chaotic solutions could occur is also very small and of little practical importance in the present context. Moreover, it must be remembered that the non-inertial form of Darcy’s law was used in ref. [ 123 and thus the chaotic motion found there must be due to some other reason, either the boundary conditions used or two dimensionality. This merits further study.
CONCLUSIONS
The governing equations for time-dependent natural convection arising from an arbitrary heat source distribution in a thin annulus have been reduced to the solution of two non-linear ordinary differential equations. The approximation does not involve any restriction on the Rayleigh number though effects due to possible recirculation zones within the annular thickness have been suppressed. For infinite Rayleigh numbers and symmetrical heating an exact solution exists. This solution, combined with a linear stability analysis, provides information on the general behavior of the system. For asymmetric heating or finite Rayleigh numbers, damped velocity oscillations can be observed leading to a steady-state motion. The damping factors can be qualitatively predicted from the eigenvalues of linear stability theory. Under cer-
737
tain conditions two stable steady states are possible, and either can be reached depending on the initial conditions. It can also be demonstrated that chaotic behavior is not possible under the thin annulus approximation unless fluid inertia effects are included. Two other observations can be made with respect to possible numerical solutions of the governing partial differential equations under these and similar circumstances. Any results implying steady oscillations should be carefully verified in order to distinguish between strictly periodic phenomena and slowly damped oscillations. Also, in time-dependent calculations, different initial conditions can lead to different steady states.
Acknowledfement-We acknowledge the oartial suooort of the Na&nal Science Foundation _ under _-Grant MEA 8401489.
REFERENCES 1. P. Cheng, Heat transfer in geothermal systems, Adu. Heat Transfer 14, l-105 (1978). 2. M. A. Combamous and S. A. Bories, Hydrothermal convection in saturated porous media, Ado. Hydrosci.
10,231-307 (1975). 3. H. H. Bau, Low Rayleigh number thermal convection in a saturated porous medium bounded by two horizontal eccentric cylinders, ASME J. Heat Transfer 106, 166175 (1984). 4. H. H. Bau, Thermal convection in a horizontal, eccentric
annulus containing a saturated porous medium-an extended perturbation analysis, Int. J. Heat Mass Transfer 27, 2277-2287 (1984). 5. Y. Takata, K. Fukuda,
S. Hasegawa, K. Iwashige, H. Shimomura and K. Sanokawa, Three-dimensional natural convection in a porous medium enclosed with concentric inclined cylinders, Proc. 7th Int. Heat Transfer Conf., Munich, Vol. 2, pp. 351-356 (1982). 6. R. Echigo, S. Hasegawa, S. Tottori, H. Shimomura and Y. Okamoto, An analysis on the radiative and free convective heat transfer in a horizontal annulus with permeable insulator, Proc. 6th Int. Heat Transfer Conf., Toronto, Vol. 3, pp. 385-390 (1978). 7. P. Vasseur, T. H. Nguyen, L. Robillard and V. K. T. Thi, Natural convection between horizontal concentric cylinders filled with a porous layer with internal heat generation, Int. J. Heat Mass Transfer 27, 337-349 (1984). 8. V. A. Brailovskaya,
G. B. Petrazhitskii and V. I. Polezhaer, Natural convection and heat transfer in porous interlayers between horizontal coaxial cylinders, J. Appl. Mech. Tech. Phys. 19, 781-786 (1978). 9. P. J. Burns and C. L. Tien, Natural convection in porous media bounded by concentric spheres and horizontal cylinders, Int. J. Heat Mass Transfer 22,929-939 (1979). 10. J.-P. Caltagirone, Thermoconvective instabilities in a porous medium bounded by two concentric horizontal cylinders, J. Fluid Mech. 76, 337-362 (1976). 11. G. N. Facas and B. Farouk, Transient and steady-state natural convection in a porous medium between two concentric cylinders, ASME J. Hear Transfer 105, 660663 (1983). 12. L. Robillard, T. H. Nguyen, M. G. Satish and P. Vasseur,
Oscillating flow in an annular porous layer heated from
M. SEN and K. E. TORRANCE
738
a++-, -_=I)
below, Heat Transfer in Porous Media and Particulate Flows (edited by L. S. Yao, M. M. Chen, C. E. Hickox,
P. G. Simpkins, L. C. Chow, M. Kaviany, P. Cheng and L. R. Davis), HTD-Vol. 46, pp. 41-47. The American Society of Mechanical Engineers, New York (1985). 13. T. H. Nguyen, M. G. Satish, L. Robillard and P. Vasseur, Free convection in a non-uniformly heated annular porous layer, ASME Paper No. 85-WA/HT-8. The American Society of Mechanical Engineers, New York (1985). 14. J. A. Yorke and E. D. Yorke, Chaotic behavior and fluid dynamics, Hydrodynamic Instabilities and the Transition to Turbulence (edited by H. L. Swinney and J. L. Gollub), 2nd edn, pp. 77-95. Springer, Berlin (1985). 15. M. Sen, E. Ramos and C. Treviiio, The toroidal thermosyphon with known heat flux, Int. J. Heat Mass Transfer 28, 219-233
az
-
’
n;+2
= 0,
2+as=o (31a-c)
uo = Kp*gBT;/&, II: = -nKpyRp
vf, = nKp’JRp, w.
=
---
K ah p az'
w:,=_--2, K aps
wc” =
KapC,
F aZ ~2
+ -!-x(-nv:TE, 2R n
+ nuC,T:) + w,,$
=a2
ar nv c._!-__?Te,fw at R
(33a-c)
P az
(1985).
16. T. Masuoka, M. Tashiro and T. Katsuhara, Transient natural convection in a porous medium within a closed sphere, Proc. JSME-ASME Thermal Engng Joint Conf., Honolulu, pp. 243-249 (1983).
(32a-c)
ar ---‘+w:o O a2
a*T,
(34a)
aT
aZ
=q:-!$~+a!$
(34b)
APPENDIX
Here we show that the present study is applicable not only for a planar geometry but also for a three-dimensional annulus with a horizontal axis, as shown in Fig. 1. This is due to the fact that under the thin annulus approximation, the assumption of impermeable and adiabatic boundaries leads to zero flow and heat conduction in the axial direction. For a thin annulus with negligible radial velocity, we have the following governing equations for flow in porous media. Mass conservation gives
f!+g=O
- B(T- T*)JgcosO
!!%T’+w ’
az+wcs
O az
” az
(34c) Equations (31a)-(31c) can be integrated in the z-direction to give wg = 0
(294
where w is the axial component of the fluid velocity. The Darcy-Oberbeck-Boussinesq equations in the circumferential and axial coordinate directions are respectively
u = -$p’(l
aTE,
‘dt+R
+k$]
(29b)
dz = 0
{2} = ${z}sin(Znmz/L)
k ap
+ {g} The Fourier components expanded as
The energy equation is [I]
{z} The four dependent variables II,w, T and p are Zn-periodic in 0 so that they can be expanded in Fourier series of the form
(36)
where impermeability at z = 0 and L has been used. Equation (36) indicates that Fourier series expansions of vi and vi in z may also be used. These can be written as
w=-lz’
I$ = & + C [& sin (no) + 4’. cos (no)] n
(35)
COS(2nnu,Q].
(37)
of the axial velocity can be
= ~{~‘“}sin(2nrnz/~).
(38)
Expansions (37) and (38) can be substituted in equations (31b) and (31~) to obtain yes nn = y’s nm= 0
(39)
(30)
where 4 represents each one of the variables. The new dependent variables & and +‘, are functions of z and t alone. The summations are from n = 1 to infinity, as will be all the other summations. The heat source can also be expanded as in equation (5). Three component equations can be obtained for each one of equations (30), one for the integral itself, and then two for the integral of sin(m0) and cos(m0) times the equation. These equations are
In addition equations (31b), (31c), (32b), (32c), (33b) and (33~) can be used to obtain w’ I
R2azw: 2 a2
0
’
wc - $2
”
= 0.
(41a,b)
Substituting equation (38) in equations (41a) and (41b) and ,.^~ using equation (4U),we get
Natural
convection
in a thin horizontal
porous annulus
739
Thus, under the thin annulus approximation, ‘j$&
+ @$){~}sin(2nnz/L) “In
= 0
component
(43)
It can be shown in a similar manner from the energy equations (34a)-(34c) that, as a consequence of the adiabatic condition we have assumed at z = 0 and L, the temperature is independent ofz. Thus, only the circumferential component of the velocity field exists and the temperature varies only
(44)
with 9.
from which
This gives
CONVECTION
NATURELLE
the velocity
(42)
DANS UN ANNEAU
in the z-direction
is zero.
POREUX HORIZONTAL
considtre les modes possibles de coinvection naturelle, dependant du temps, dans un espace finie. L’anneau est rempli de matiere poreuse et son epaisseur est supposee petite en comparaison du rayon moyen. Toutes les frontieres sont impermeables et adiabatiques; le chauffage est fourni par une source volumetrique de chaleur distribuie circonferentiellement. Les equations se reduisent a un systcme de deux equations differentielles non lintaires. Des oscillations stables non lineaires existent pour le cas special d’un nombre de Rayleigh intini et un chauffage symttrique par rapport a la verticale. Pour des nombres de Rayleigh plus faibles, on obtient des oscillations amorties dont le degre d’amortissement augmente avec I’inclinaison de la ligne de symttrie et quand le nombre de Rayleigh diminue. Des etats stables permanents sont obtenus pour des petites inclinaisons. Des mouvements chaotiques ne se developpent pas dans les Ccoulements de Darcy non inertiels. R&urn&On
annulaire de longueur
NATtiRLICHE
KONVEKTION IN EINEM DtiNNEN, PORBSEN RINGSPALT
WAAGERECHTEN
Zusammenfassung-Die
moglichen Zustlnde der zeitabhangigen natiirlichen Konvektion in einem horizontalen Ringspalt endlicher Lange werden betrachtet. Der Ringspalt ist mit porosem Material gefiillt, und die Ringspaltbreite wird als klein im Vergleich zum mittleren Radius angenommen. Alle Grenzflachen sind dicht und adiabat ; die Beheizung erfolgt durch eine iiber den Umfang verteilte volumetrische Warmequelle. Die Bilanzgleichungen reduzieren sich auf zwei nichtlineare gewohnliche Differentialgleichungen. Im Spezialfall bei unendlicher Rayleigh-Zahl und symmetrischer Beheizung iiber die Vertikale treten stetige nichtlineare Schwingungen auf. Bei niedrigen Rayleigh-Zahlen treten gediimpfte Schwingungen auf, wobei der Grad der Dampfung mit der Neigung der Symmetrieachse und mit abnehmender Rayleigh-Zahl ansteigt. Bei kleineren Neigungen ergeben sich mehrere stabile Zustande. Bei nicht tragheitsbedingten Darcy-Stromungen entwickeln sich keine chaotischen Bewegungen.
ECTECTBEHHAg
KOHBEKHHR B TOHKOM TOPM30HTAJIbHOM KOJIbHEBOM KAHAJIE
AnHoTaluln-PaccMaTpesalorcn
IIOPHCTOM
B03MOEHbIe ~WiMbl eCTeCTBeHHOi? KOHBeKUWli B TOPR30HTanbHOM KonbueBoM Kaaane KoHeqHoii nnuHbI. KOJIbIIeBOfi KaHan 3anonHeH nopHcTbIM MaTepHanoM, ero TonmliHa cqHTaeTca Manoii no cpaeseaew, co cpenHHM paneycoM. Bee rpaHHubI HenpomiUaeMble H anHa6aTHrecKHe; HarpeB OCyIUeCTBJTaeTCa06aeMHbtMn BCTOYHWKaMH TelTJIa, paCllOnO)l(eHHbIMWII0 OKPy~HOCT1(.O~Pe~eJDI‘OUIHe ypaBHeHHK CBeneHbl KCACTeMenByX HenHHetiHbIX06blKHOBeHHbIX nH+j,ePeHUHanbHbIX ypaBHeHk& B Cny'iae 6eCKOHevHOrO WiCJIa kWIe!4 N CHMMeTpHHHOrO HarpeBa HaBnIonamrca ycToiiHHBbIe HenHUeiiHbIe KOne6aHHa. J&m 6onee HW~KWXqricen Psnea nonyqeno 3aryxanne KOne6aHHti,npWieM CTeneHb 3aTy3aHHa BO3paCTaeT C HaKJIOHOMnHHUA CHMMeTpHH H C yMeHb"IeHHeM wfcna Psnen. ,&a ManbIx yrnos HaKnoHa Ha6nmnanHCb MHoromcneHHbIe ycToi+isBbIeCOCTOIIHHR.B CJIy’Iae 6e3bIHepUHOHHbIXTe’IeHHfi flapC!H XaOTH'teCKHeLIBH)KeHAII He pa3BHaaEOTCR.