내부열원을 가지는 다공성매체 내에서의 자연대류 현상

최응수

포스코 기술연구소 경상북도, 포항시, 남구, 괴동동 1, 790-785

A BIFURCATION STUDY ON NATURAL CONVECTION

IN POROUS MEDIA WITH INTERNAL HEAT SOURCES

Eung-Soo Choi

POSCO Technical Res. Labs. 1 Koedong-dong, Nam-ku, Pohang-shi, 790-785

INTRODUCTION

As an analogue of the Rayleigh-Ben□rd problem [Ben□rd, 1901] Natural convection in fluid saturated porous media has been extensively studied over the last decades since it has numerous practical applications such as utilization of geothermal energy, fiber and granular insulations, design of packed bed reactors and underground disposal of nuclear waste materials. Earlier works of natural convection in porous media focused on the marginal stability of basic conduction state and heat transfer rate on the temperature difference across porous media [Lapwood, 1948; Elder, 1967]. Based on the steady state natural convection studies the dynamic behavior of convection flow, in which regularly oscillating modes develope into a chaotic mode as the instability of the fluid flow is increased, was explored [Caltagirone, 1975; Kimura et al., 1986]. Recently, solution structures which consist of various solution branches were built by employing the bifurcation theory [Riley and Winters, 1989; Weinitschke et al., 1990]. The bifurcation studies revealed a fundamental picture that how natural convection developes as the buoyancy force is increased. Most of the studies on natural convection in porous media used the Darcy equation to describe flow behavior. The internal geometry in porous media is too complex to apply the Navier-Stokes equation. However, the Darcy equation, which is a simple empirical relation between the discharge rate of fluid and the pressure gradient in a fictitious infinite porous layers, does not satisfy the no-slip condition at the bounding surfaces of porous media. This limitation may result overestimation of the heat transfer rate by convective flow, especially when the hydrodynamic layer is larger than the thermal boundary layer, since the damping effect by the viscous force near the solid boundary of porous media is neglected. To overcome this problem the Brinkman's extension of the Darcy equation is used in the present study. The solution structures for the problem of natural convection in fluid saturated porous media having internal heat sources are constructed using the bifurcation algorithm. The effects of the viscous force term on the natural convection is of primary interest.

GOVERNING EQUATIONS

Consider a horizontal, square, two-dimensional porous medium of width [eqn: 2a], height [eqn: 2b] and permeability [eqn: K]. The porous medium is saturated with a fluid of density [eqn: rho] and viscosity [eqn: mu].  heat generation [eqn: Q sub g] is assumed and the boundary of the porous medium is maintained at a uniform temperature [eqn: T' sub w]. The dimensionless conservation equations of continuity, momentum and energy for a two-dimensional flow field are:

[eqn: Ra { PARTIAL theta } over { PARTIAL x } ~+~ OMEGA ~-~ Da DEL sup 2 OMEGA ~=~ 0 ~~~~~~~~~~~~~~~~~~~~~(1)]

[eqn: DEL sup 2 PSI ~+~ OMEGA ~=~ 0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(2)]

[eqn: DEL sup 2 theta ~-~ left ( { PARTIAL PSI } over { PARTIAL y} { PARTIAL theta } over { PARTIAL x} ~-~ { PARTIAL PSI } over { PARTIAL x} { PARTIAL theta } over { PARTIAL y} right ) ~+~ 1/4 ~=~ {PARTIAL theta} over {PARTIAL tau} ~~(3)]

with computation domain, [eqn: R ~ := ~ left [ left ( x,y right ) vert -1 <= x <= 1, ~-1 <= y <= 1 right ] .] The Boussinesq approximation [eqn: rho sub f ~=~ rho sub w left [ 1 ~-~ beta left ( T' - T' sub w right ) right ]] is invoked to allow for density variation with temperature in the gravity term in the momentum equation. The dimensionless variables are defined as follows.

[eqn: x = x'/a, ~~~ y=y'/b , ~~~ t=t'/tau , ~~~u = u' over {alpha/a} , ~~~ v = v' over {alpha/b}]

[eqn: PSI = PSI' over alpha , ~~~ OMEGA = OMEGA' over {alpha/a sup 2}, ~~~ theta = {left ( T' sub w -T' right ) } over {Q sub g A / k}]

The control parameters appearing in Equations (1)∼(3) are the Rayleigh number [eqn: Ra = K bold g beta Q sub g A a /k alpha nu], the Darcy number [eqn: Da = K / a sup 2]. On the solid boundary, the Dirichlet type boundary conditions are used for the temperature and the stream function [eqn: theta = 0, PSI = 0] and the vorticity can be obtained by taking the Taylor series expansion for the stream function [Roache, 1972].

NUMERICAL METHODS

The non-linear boundary value problem defined by Equations (1)∼(3) can be represented in the following implicit form

[eqn: F left ( bold u, Ra, Da right ) ~=~ 0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(4)]

where [eqn: bold u ~: = left ( theta left ( x,y right ) , OMEGA left ( x,y right ) , PSI left ( x,y right ) right ) .] The non-linear Equation (4) can be solved by using the Newton-Raphson method in which the successive correctors, [eqn: U sub n = left ( T left ( x,y right ) , O left ( x,y right ) , P left ( x,y right ) right )] are sought by solving the linear boundary value problem

[eqn: F sub bold u left ( bold u sub n , Ra, Da right ) U sub n ~=~ bold {r} sub n ~~~~~~~~~~~~~~~~~~~~~~~~~~(5)]

where [eqn: F sub bold u] is the Fr□ch□t derivative and[eqn: bold {r sub n}] is the residuals of Equations (4). The predictor which is needed for the next step [eqn: Ra + TRIANGLE Ra] can be generated by employing the Euler-Newton method

[eqn: bold u left ( Ra + TRIANGLE Ra right ) ~=~ bold u left ( Ra right ) ~+~ { PARTIAL bold u } over { PARTIAL Ra } TRIANGLE Ra ~~~~~~~~~~~~(6)]

where [eqn: PARTIAL bold u / PARTIAL Ra] is obtained from

[eqn: F sub bold u left ( bold u , Ra , Da right ) { PARTIAL bold u } over { PARTIAL Ra } ~=~ -{ PARTIAL F } over { PARTIAL Ra } . ~~~~~~~~~~~~~~(7)]

When there is a limit point nearby, the regular continuation scheme fails. This can be overcome by parametrizing the problem with the pseudo-arclength parameter[eqn: s], measured along the solution path. Then, the dynamic parameter Ra is treated as an element of the unknown vector[eqn: left ( theta left ( s right ) , OMEGA left ( s right ) , PSI left ( s right ) , Ra left ( s right ) right ).] Various parametrizing equations have been proposed in the literature. We use the one similar to that proposed by Keller[20]. viz.

[eqn: w INT INT left [ left ( theta left ( s right ) - theta left ( s sub o right ) right ) sup 2 ~+~ left ( OMEGA left ( s right ) - OMEGA left ( s sub o right ) right ) sup 2 ~+~ left ( PSI left ( s right ) - PSI left ( s sub o right ) right ) sup 2 right ] dxdy]

[eqn: + left ( 1-w right ) left ( Ra left ( s right ) - Ra left ( s sub o right ) right ) sup 2 ~=~ 0 ~~~~~~~~~~~~(8)]

where [eqn: w] is a weight factor, used to give equal weight to all of the variables that go to define the arclength.

COMPUTATION RESULTS

Figure 1 shows the solution structure which consists of various solution branches for [eqn: Da = 10 sup -3] and up to [eqn: Ra=15000]. The dotted lines are non-symmetric solutions and they occur in pairs, with one being the mirror image of the other. Since [eqn: N u] is an averaged quantity both parts of an non-symmetric branch have the same value and hence they appear as a single curve in these bifurcation diagrams. A unique solution on the primary branch, PM1 was easily obtained for [eqn: Ra=0] and then the entire branch of PM1 was generated using the arclength continuation scheme over the half domain [eqn: R sup +]. The fist non-symmetric solution branch, AS1, which bifurcates at the two symmetry breaking bifurcation points, SB1, SB2, was found by introducing non-symmetric perturbation to the Euler-Newton method over the full domain [eqn: R]. The isolated symmetric solution branches, IS1 and IS2, and the non-symmetric solution branches, AS2 and AS3 were located by solving the extended system. It is seen that all the symmetry breaking bifurcations happen near the limit points. However, the limit points and the symmetry breaking bifurcation points do not coincide since the symmetry breaking bifurcation occurs as a pitchfork. The multiple solutions start to appear[eqn: Ra=2586]. Each of the different solution branches in the bifurcation diagram correspond to different flow patterns. Across any singular point, there is a qualitative change in the flow pattern. In Figure 1 a total of five different solutions are found at [eqn: Ra= 7000]. Three symmetric solution branches are on PM1 and IS2 while two non-symmetric solutions are on AS1. The contours of the streamlines and isotherms for the symmetric four cell solution and non-symmetric two cell solution among the five different solutions are shown in Figure 2.

Figure 1

Figure 2

REFERENCES

1.	Ben□rd, H., Les tourbillons cellulaires dans une nappe liquide transportant       chaleur par convection en regime permanant. Ann. Chem. Phys, 7(Ser.        23), 1901.

2.	Caltagirone, J. P., Thermoconvective instabilities in a horizontal porous	layer, J. of Fluid Mech., 72, 169, 1975.

3.	Elder, J. W., Steady free convection in a porous medium heated from	below, J of Fluid Mech., 27, 29, 1967.

4.	Keller, H. B., Applications of bifurcation theory(ed. P. H. Rabinowitz).	Academic, New York, 1977.

5..	Kimura, S., Shubert, G. and Straus, J. M., Route to chaos in porous-	medium thermal convection, J of Fluid Mech., 166, 305, 1986.

6.	Lapwood, E. R., Convection of a fluid in a porous medium, Proc. Camb.	Phil. Soc., 44, 508, 1948.

7.	Riley, D. S. and Winters, K. H., Modal exchange mechanism in Lapwood	convection, J. Fluid Mech., 204, 325, 1989.

8.	Roache, P. J., Computational fluid dynamics. Hermosa, 1972.

9.	Weinitschke H. J., Nandakumar K. and Ravi Sankar S., A bifurcation	study of convective heat transfer in porous media. Phys. Fluids A2, 912,	1990.