- Email: [email protected]
Rayleigh–Bénard stability of a solidifying porous medium
International Journal of Heat and Mass Transfer 42 (1999) 3337±3350
Rayleigh±BeÂnard stability of a solidifying porous medium C. Mackie a, P. Desai b,*, C. Meyers c a
Tulane University, School of Mechanical Engineering, New Orleans, LA 70118, USA Georgia Institute of Technology, School of Mechanical Engineering, Atlanta, GA 30332, USA c North Carolina A&T, College of Engineering, Greensboro, NC 27411, USA
b
Received 17 February 1998; received in revised form 8 October 1998
Abstract In many solid/liquid phase transformation processes, natural convection controls the freezing or melting rate of the material. Freezing or melting in the presence of nonmelting components is important in numerous industrial processes and is also present in a wide range of systems in nature. The kinematics of the solid±liquid interface, coupled with bulk convection in the melt, play an important role in determining the microstructure of a solidi®ed material. The fundamental problem of thermoconvective instability of a single-component ¯uid in a horizontal porous layer has been studied extensively; additional complexities arise during solidi®cation due to the presence of a non-melting component, namely the porous medium. This study addresses the problem of Rayleigh±BeÂnard instability of a liquid layer undergoing a phase transformation within a porous medium. A linear stability analysis is performed to determine the eects of the medium and phase-change on the conditions for incipient instability under a range of thermal boundary conditions. The analysis reveals that the onset of convection, or the stability of the system, is signi®cantly aected by the presence of the porous medium, state of solidi®cation and the thermal boundary conditions. # 1999 Elsevier Science Ltd. All rights reserved.
1. Introduction The problem of thermoconvective instability in a horizontal ¯uid layer driven by buoyancy eects studied extensively in the literature [1±3] encompasses unresolved complexities associated with the freezing or melting of a liquid in the presence of a nonmelting component, such as a solid matrix or soil. This study addresses hydrodynamic stability issues during solidi®cation of such mixtures by modeling the matrix±liquid slurry as a liquid saturated horizontal porous medium in a gravitational ®eld going through a phase change under a full range of thermal conditions at the lower bounding surface. Freezing or melting in
* Corresponding author. Tel: +1-404-894-3233; fax: +1 404-894-7790. E-mail address: [email protected] (P. Desai)
the presence of nonmelting components is important in industrial applications as well as in a wide range of systems in nature. Examples include, but are not limited to, seasonal freezing and melting of soil, lakes and rivers, arti®cial freezing of the ground as a construction technique for supporting poor soils, insulation of underground buildings, the melting of the upper permafrost in the Arctic due to buried pipelines, thermal energy storage in porous media and storage of frozen foods. Metallurgical applications include manufacturing of composite materials and puri®cation of metals. The study of freezing and melting in porous media has gained a great interest over the past few years. A fundamental understanding of convective stability in such a range of applications is thus very important. It is recognized that ¯uid ¯ow, in general, and buoyancy-driven ¯ow, in particular, plays a critical role in determining interfacial traits of a solidifying material;
0017-9310/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 7 - 9 3 1 0 ( 9 8 ) 0 0 3 6 1 - 5
3338
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
Nomenclature a ac A BL ca ci cp D Da f g g hSL h1 kÃ, zà kL,S K L n, tÃi Pr Ra Rac (Rac) Rc (Rc) St T T1 T10 TM v w wà z
separation constant or wave number critical wave number solid thickness ratio Biot number, h1Z0/kL acceleration coecient general constants for characteristic equation speci®c heat at constant pressure dierential operator, D=d/dz, D 2=d2/dz 2 Darcy number, K/Z 20 temperature perturbation variable (liquid), see Eq. (33) temperature perturbation variable (solid), see Eq. (33) gravitational acceleration vector speci®c latent heat of solid heat transfer coecient of the bulk environment unit vector in the z-direction thermal conductivity for liquid and solid, respectively permeability of the porous medium length of the solidi®cation system unit normal and tangent vectors Prandtl number, nL/aL Rayleigh number, gb(TLÿTM)Z 30/nLaL critical Rayleigh±Darcy number, Rc Da critical Rayleigh number Stefan number, hSL/cpDTL temperature upper plate temperature lower plate temperature interface melting temperature velocity vector for liquid, v=v(u, v, w ) velocity component in the z-direction for the liquid vertical velocity perturbation variable spatial coordinate in vertical direction for liquid.
Greek symbols a ratio of the thermal diusivities of the liquid and the solid, aL/aS b volumetric coecient of thermal expansion acceleration coecient, ca/sH Pr ga DTL temperature dierence across the porous layer Z interface position coordinate Z normalized perturbation variable for the interface position Z0 length of the liquid saturated porous region time rate of change of the interface Zt m me/mL mL dynamic viscosity nL kinematic viscosity r ratio of the densities of the liquid and the solid, rL/rS s real part of the complex growth rate [s] sH heat capacity ratio, see Eq. (3) f medium porosity F pattern planform, where F=F(x, y ) H2 Laplacian operator H2H 2-D Laplacian operator, see Eq. (35).
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
Subscripts e L M S 0 1 Superscripts '
3339
eective property for porous medium liquid property or region melting temperature solid property or region reference state bulk environment. bar, basic state quantity perturbation variable normalized.
¯uid motion adjacent to a solidifying interface aects the local thermal and solutal ®elds which control the geometric, dynamic and thermodynamic characteristics of the interface. For instance, convection during fabrication or solidi®cation processing of metal matrix composites aects the distribution of particulates in the composites, in turn in¯uencing the properties of the material. In the pure metal analog, independent of the solute diusion, it is also known that the convective stability is aected by the density and presence of particulates as well as by the conductivity of the boundary. The porous medium model simulates the in¯uences of a non-interacting matrix, coupled with a range of thermal conditions at the lower boundary, on the initiation and growth of bulk convection in a solidifying nonhomogeneous system. This study focuses on the eect of the permeability of a porous medium and phase change on the incipient convection in the interstitial liquid; by parametrically changing the permeability from that corresponding to a Darcy model of the ¯ow through it to a vanishing one, an understanding is gained for the in¯uence of volume fraction and its distribution on the critical conditions for convective stability of the interstitial ¯uid in porous media in the presence of solidi®cation. This work couples solidi®cation dynamics with the thermal and hydrodynamic ®elds in a nonhomogeneous solidifying system via a Rayleigh±BeÂnard stability analysis. The porous medium is saturated with the liquid in all available pores, or the interparticle space of the medium, leaving no space for a vapor, gas, or other liquid within the medium Since the pioneering work on the porous-medium analog of the Rayleigh±BeÂnard stability problem [4], the subject of phase-change in a porous medium has received increasing attention over the past two decades. While liquid/gas phase change in a saturated porous medium has received considerable attention, solid/ liquid phase change in liquid-saturated porous medium
has received considerably less attention. Such problems and their manifestation represent a relatively new area in the ®eld of convection in porous media [5]. This study considers a horizontal, liquid saturated, porous medium undergoing solidi®cation while being cooled uniformly from above. A recent study [6] examined the coupled eects on stability arising from thermal convection and solidi®cation of a single-component liquid in a porous medium. The stability of the basic state of heat conduction and the stability of ®nite-amplitude convection in the Darcy limit of vanishing permeability for isothermal boundaries were examined via a twoparameter perturbation analysis. The present work investigates the eects of the solid matrix in a pure melt on the thermal convection for a range of permeabilities and the full range of lower wall thermal boundary conditions. The cited references are based on the Darcy ¯ow model, and hence are relevant to well packed low permeability porous media. This study diers by including the Brinkman extended ¯ow model that allows for the investigation of sparsely packed media. Beyond the context of convective stability, there exists substantial research in porous media transport phenomena germane to the present work. NonDarcian eects have been modeled illustrating the transition between the porous and clear ¯uid limits by modifying the Darcy equation with the Brinkman term [10]. Also experiments have been conducted using a modi®ed Darcy equation incorporating the Forchheimer term which attempts to account for microscopic drag forces due to acceleration [11]. It was demonstrated experimentally [7] that an eective porous medium thermal conductivity model is adequate for heat transfer estimates only when the solid matrix and pore-material have suciently similar thermal conductivities to allow for the requisite local thermal equilibrium. A continuum model of the porous medium was developed [8] for analyzing solid/liquid phase
3340
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
change in a porous medium that accounted for inertia as well as boundary friction eects. In particular, these include the Brinkman viscous and the Forchheimer pressure drag terms. Such results are relevant to the present work where the liquid and the porous material have similar conductivities, and where friction due to microscopic shear is also important.
2. The mathematical model A quiescent horizontal layer of a pure ¯uid in a gravitational ®eld heated from below does not always remain so in the presence of an adverse density gradient beyond a threshold value which is capable of giving rise to ¯uid motion. This observation may be applied to examine the solidi®cation (or melting) of a porous layer saturated with a pure liquid enclosed between two rigid, parallel, thermally dissimilar, horizontal plates of in®nite extent, a distance L apart, in a gravitational ®eld, as shown in Fig. 1. The analytical model is formulated under the assumptions that the liquid in the upper part of the layer is frozen and the eective thermophysical and transport properties of the porous medium are homogeneous and isotropic. The porous matrix and the solidi®ed layer are both assumed to be rigid [8]. Fluid motion within the por-
ous medium is described via a Newtonian, Boussinesq model of the Navier±Stokes equations, written for a liquid saturated porous medium. Motion due to density change upon phase change is not considered to be signi®cant. The liquid velocity through the medium is deemed suciently small, conforming to the Darcy regime, or, at most, to the Brinkman-modi®ed Darcy regime. The melt±freeze front is assumed to be a thin surface of negligible thickness that remains at the melting point of the phase-change material. The ¯uid layer is separated at the bottom from the surroundings by a thin rigid lamina of negligible heat capacity. The upper plate is maintained at a constant temperature, T1. 2.1. Governing equations In the porous layer, the heat transfer at the lower surface of the system is via a convective±conductive mechanism. The ¯ux at the bottom plate is maintained such that the melting temperature, TM, is bounded between the lower plate temperature, T10, and the upper plate temperature, T1 (T10>TM>T1). Therefore, there is a solid±liquid interface at z=Z0, 0 < Z0 < L. The governing equations are written in a Cartesian frame ®xed with respect to the composite system on the lower plate with the positive z-axis pointing
Fig. 1. Solidi®cation model with saturated porous medium.
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
upward (see Fig. 1). The temperature ®eld within the solidi®ed layer is described by the thermal energy equation @ TS aS r2 TS : @t
1
For the liquid, the temperature ®eld may be obtained from a simultaneous solution of the thermal energy and continuity equations, together with the Boussinesq model of the Darcy-extended momentum equation. These are expressed, respectively, as sH
@ TL v rTL aL r2 TL @t
2
where
2.2. Boundary conditions The boundary conditions at the phase boundary are the continuity of temperature TL TS TM
6
where the interface curvature due to undercooling of the melt is neglected, and the conservation of energy rS hSL
@Z j kS rTS ÿ kL rTL j n: @t
3
rv0
4
@Z kn @t
ca rL
@v ÿrP ÿ @t
mL v meff r2 v K
5
rL b TL ÿ Tref gk^ where ca is a constant which depends on the geometry of the porous medium and is determined by the nature of the assumed model pore tubes [5]. In most cases, the term becomes negligible when suitably scaled. This is due to the small value of permeability which multiplies the scaled inertia term. In the proceeding governing equations, the Brinkman-extended Darcy momentum model is employed. The frictional eects represented by Darcy law in the ¯uid momentum equation in a porous medium, are two orders of magnitude lower than those represented in the Navier±Stokes equations. Thus, only the impermeability condition at the surface is satis®ed, leaving the no-slip boundary condition unresolved. The Brinkman-extended Darcy model satis®es not only the impermeability and no-slip conditions but accounts also for the friction caused by microscopic shear by introducing an additional viscous force term via an eective dynamic viscosity, meH2v. The Brinkman-extended model reduces to a form of the Navier±Stokes model as the medium permeability approach in®nity (K 4 1) and to the conventional Darcy model as the medium permeability approach zero (K 4 0). Such a formulation allows a comparison of the ¯ows in porous media with those in clear ¯uids.
8
and the no-slip condition v t1 v t2 :
and
7
Also de®ned on this boundary (interface) are the conservation of mass (the kinematic condition) rL v n rL ÿ rS
f rcr L 1 ÿ f rcr S sH heat capacity ratio rcr L
3341
9
The boundary conditions on the two bounding surfaces are dTL at z 0, ÿ kL h1 T1 ÿ T10 dz 0 10 and at z L,
TS T1 :
3. Basic state and linear stability analysis Initially, the stationary system is in a static equilibrium state with no ¯uid motion and a planar interface at z=Z0. The purely conductive temperature pro®les in the liquid and the solid regions are, respectively, given by
TL
ÿh1 T1 ÿ TM z T1 ÿ TM h1 Z0 TM h1 Z0 kL h1 Z0 kL
11
TM ÿ T1 Z TM ÿ T1 zÿ 0 TM : Z0 ÿ L Z0 ÿ L
12
and
TS
The static equilibrium temperature pro®les are used to normalize the solid and liquid temperatures, respectively, as T L
TL ÿ TM ; T10 j ÿTM
T s
TS ÿ T1 : TM ÿ T1
13
3342
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
The scales for, v, x (and y, z ), t and p are chosen, respectively, as
distribution. In normalized form the basic state becomes (denoted by variables with bars on top)
aL sH Z20 ,Z0 , Z0 aL
v 0;
and
rL nL aL : Z20
The normalized equations become, after dropping the asterisks, for the solid, 1 @ TS aS 14 r 2 TS sH @ t aL and for the liquid, @ TL v rTL r2 TL , @t
rv0
15a,b
and @v ga rP @t
1 v ÿ m~ r2 v Ra TL k 0: Da
16
The normalized boundary conditions on the interface are v t1
1 ÿ r Zt k n, sH
r St Pr
@Z ArTL ÿ rTS n, @t
v1 t1 v2 t2 0
17a,b
and at z L,TS 1:
19a,b
It may be noted that the dimensionless group, BL, acts like a Biot number for heat transfer from the bulk environment to the liquid. As BL approaches zero, the boundary condition at the surface approaches that of an insulated surface; as BL approaches in®nity, the boundary condition approaches that of a constant surface temperature; also appearing is another dimensionless parameter: A
L ÿ Z0 kS DTS Z0 kL DTL
T L 1 ÿ z;
zÿ1 : T S 1 ÿ Lÿ1
22
It is noted that TL and TS are scaled by dierent reference lengths. Thus at z=1, the temperature for the solid and liquid regions (TL and TS) are equal to the melting temperature, TM [Eq. (12)]. Next, a linear perturbation model of the system is developed by considering a slight change in temperature due to a slight melting of the solid. Such a problem is formulated by introducing perturbations of the basic state de®ned via a conventional barred-primed scheme of notations as TS T S T 0 S ,
TL T L T 0 L
p p p 0 ,
23
Z 1 Z 0:
24
18a,b
The surface boundary conditions become dTL BL TL , dz
21
and
v v 0,
and TL TS 0:
at z 0,
dp Ra T L k dz
20
where z=Z0 is the position of the planar interface and z=L is the thickness of the combined solid±liquid system. For a small A, A ÿ1 may be considered the equivalent Biot number for heat transfer from the liquid to the solid; also, A measures the amount of solid present in the system [9]. The basic state under consideration is that of a stagnant layer of liquid, with a hydrostatic pressure distribution, and a purely conductive temperature
Introduction of the perturbation expressions, Eqs. (23) and (24), into the governing Eqs. (14)±(16), followed by linearization in the disturbance quantities, yields the perturbation model for the solid, 0 1 @T S aS r2 T 0 S 25 sH @t aL and for the liquid @ T 0L w 0 r2 T 0 L , @t
r v0 0
26
and ga
@v0 rp 0 @t
1 v 0 ÿ m~ r2 v 0 Ra T 0 L k 0 Da
27
where w' is the vertical component of the perturbation velocity v '. The linearized boundary conditions for the perturbation model are at z 1,
T 0 S ÿZ 0
@ T S , @z
T 0 L ÿZ 0
0 @Z @ T 0S @ T 0L r St Pr A ÿ @t @z @z
@ T L @z
28
29
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
and u 0 1 ÿ r
@Z @x
0
,
@Z 0 v 0 1 ÿ r , @y
@Z w 0 1 ÿ r : @t
Substitution of the normal mode into the linear perturbation model [Eq. (32)] results in the normal-mode disturbance equations for the solid 30
at z 0, T
0
S
and at z L,
31
0:
@ 2 1 2 r w ÿ r w r2 r2 w Ra r2H TL 0 @t Da
@ TL ÿ w ÿ r 2 TL 0 @t r2 TS 0
aS 2 D ÿ a2 g aL
36
sf w^ D2 ÿ a2 f
The solution to the preceding set of equations results in a sucient condition to ascertain the stability boundary of the liquid layer, as would be expected from the linear analysis. The linear perturbation system, [Eqs. (25)±(31)], is examined to determine the critical and necessary conditions for incipient nonlinear convection. After eliminating the pressure from the liquid momentum equation [Eq. (27)] by taking its curl, the curl is taken again, yielding
ÿga
sg
and, for the liquid
The surface boundary conditions become dT 0 L BL T 0 L , dz
3343
(32)
37
and 2
2
2
2 2
ÿga s D ÿ a w^
1 D2 ÿ a2 w^ Da
38
2
ÿ D ÿ a w^ Ra a f where an ordinary dierential operator, D, has been de®ned as D
d dz
and D2
d2 : dz2
39
The boundary conditions transform as g L 0,
g 1
Z^ A
40
for the solid region, and ^ Df 1 ADg 1 w 1 1 ÿ rsZ
^ Dw 1 1 ÿ rsZ
^ ^ w 0 0 Dw 0 0 Df 0 BL f 0 where the primes have been dropped. Next, to ®nd normal modes solution to the linear stability problem it is noted that the variables become separable under this assumption, yielding solutions of the form ^ Zg est F x,y fTS ,TL ,w,Zg f g z,f z,w,^
33
where the time constant, s=sr+isi, contains sr, the growth rate, and si, the frequency of the disturbance, and F is the plan form function which determines the cellular structure of the ¯uid motion and satis®es the two-dimensional wave or membrane equation Fxx Fyy ÿa2 F
34
or r2H F
2
a F 0,
where
r2H
@ @ : @ x 2 @ y2
35
It may be noted that `a' is the separation constant, mathematically, in addition to being the nondimensional overall wave number, physically.
41 in the liquid region.
4. Solutions 4.1. Solid region solution It is noticed that the general solution for the solidi®ed layer may be obtained independent of the liquid equations, resulting in g z C1 sinh az C2 cosh az:
42
Use of the boundary conditions, [Eq. (41)], yields the temperature pro®le in the solid as g z
ÿZ sinh a 1 A ÿ z : A sinh aA
43
It was necessary ®rst to solve for the solid region, independent of the liquid region, in order to express all of the governing equations in the liquid region in terms of the temperature perturbation variable, f.
3344
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
4.2. Liquid region solution The ®rst two (liquid momentum and energy) equations in Eq. (32) are combined, thus eliminating the vertical velocity component, forming a single-variable perturbation equation in terms of the normal mode perturbation temperature, f. The combined expression becomes 1 1 2 2 2 ÿ ga s ÿ s D ÿ a f ga s m~ s Da Da 44 D2 ÿ a2 2 f ÿ m~ D2 ÿ a2 3 f Ra a2 f 0: Next, it is assumed that the principle of exchange of stability is valid, or that all possible values of the growth rate allowable are real. The implication of s being real is that the neutral or marginal stability curve is characterized by a growth rate, s, of zero. Accordingly, the single-variable temperature perturbation equation for the liquid [Eq. (44)] reduces to a sixth-order, linear, homogeneous equation 1 2 2 3 D ÿ a f ÿ D2 ÿ a2 2 f Ra a2 f 0 45 Da with six linear, homogeneous boundary conditions (determined by applying the solid solution) and evaluating the continuity and thermal energy equations at z=(0, 1), given as D2 ÿ a2 f 0 0, D D2 ÿ a2 f 0 0, Df 0 BL f 0
D2 ÿ a2 f 1 0 D D2 ÿ a2 f 1 0
and Df 1 a Coth aAf 1 0: 46
The linear, homogeneous, sixth-order system [Eqs. (45) and (46)] constitutes an eigenvalue problem. The system obviously has the trivial solution, f=0, for any values of a, BL, Ra, and Da, but, in general, also possesses non-trivial solutions depending on an eigenvalue, say l, and has a general solution obtained as a linear combination of six linearly independent solutions. In this case, a, BL, and Da become ®xed parameters, with Ra serving as the eigenvalue, producing non-trivial continuous solutions for the perturbation temperature eigenfunction, f(z ). Use of the six boundary conditions and the general solution, f(z ), yields six simultaneous linear equations. This de®nes the eigenvalue problem in which Ra serves as the eigenvalue for a particular set of values for a, BL, Da and A. The critical Rayleigh number for the onset of convection has a corresponding critical wave number, both of which being
the minimum values satisfying each value as an eigenvalue. The solution for the stability problem is obtained for given values of BL, Da and A, by determining iteratively the lowest characteristic values of a and Ra. 5. Results and discussion The linear stability analysis for the system in this work has two experimentally veri®able results: the critical temperature dierence for the onset of convection and the corresponding critical wavelength. For a range of boundary conditions varying from a constant temperature to a constant heat ¯ux, the incipient conditions for convection for a pure liquid solidifying in a porous medium are determined. 5.1. Eect of varying solid thickness and boundary conditions Since the solid portion of the porous medium is typically held at temperatures below the melt temperature, initial freezing begins with the contact between the medium and the superheated liquid. Physically, the removal of the latent heat takes place rapidly, creating a thin coating on the porous medium, the degree to which the permeability is aected due to such a phenomenon being quite small. Therefore, any change in permeability will be due to the solidi®cation of the interstitial ¯uid within the slightly coated matrix. In the model, cooled from above, and/or heated from below, solidi®cation occurs only at the interface; interstitial solidi®cation does not occur. Therefore, Figs. 2±8, the critical conditions for incipient convection at ®xed permeabilities are noted. As the solidi®ed portion of the system grows, or as, the aspect ratio of the layer, A, increases for a set of ®xed values of the Darcy number, Da, and the Biot number, BL, the critical Rayleigh number of the system decreases, rendering the system less stable. It may be noted, for a system with no solidi®cation (A=0), constant temperature boundary, BL=1000, and a dense medium, DaE10ÿ6, that the classical problem of convection in a porous medium between two in®nite walls, the so-called `Lapwood problem', and the solution of Rac=4p 2=39.45 are recovered. Departure from the classical problem of rigid boundaries was determined as a function of solid thickness and the results are designated on Figs. 2±7 as the Lapwood problem recovering the results of Karcher et al. Results for other conditions may be compared with those resulting from the modi®ed classical Lapwood case. Figs. 2±4 illustrate the eect of the solid thickness ratio on the critical Rayleigh number for ®xed permeability for dierent boundary conditions, compared to the modi-
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
3345
Fig. 2. Critical Rayleigh±Darcy number vs. solid thickness ratio for Da=10ÿ6.
®ed Lapwood solution, which is denoted by a thin dashed line. The ®gures also show that as the boundary condition at the bottom wall is modi®ed from that of a ®xed temperature (BL 4 1), to that of an insulated boundary (BL 4 0), the system becomes less stable. For all values of the Darcy number considered, simultaneously increasing both the solid thickness ratio and the Biot number has a strictly stabilizing in¯uence. Figs. 2±7 indicate further that the critical Rayleigh±
Darcy and wave numbers approach their lower limits asymptotically as A 4 1. Such a destabilizing eect was previously documented for solidi®cation in a bulk liquid-layer [9] as well as in a dense liquid saturated porous medium [6]. Considering the limiting case of an isothermal lower wall (BL=1000) and a densely packed medium (Da=10ÿ6), there is a 16.44% decrease in the critical Rayleigh±Darcy number as the (normalized) solid thickness, A, increases monotonically from 0±1;
Fig. 3. Critical Rayleigh±Darcy number vs. solid thickness ratio for Da=10ÿ4.
3346
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
Fig. 4. Critical Rayleigh±Darcy number vs. solid thickness ratio for Da=10ÿ3.
Fig. 5. Critical wave number vs. solid thickness ratio for Da=10ÿ6.
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
3347
Fig. 6. Critical wave number vs. solid thickness ratio for Da=10ÿ4.
for the less dense case of Da=10ÿ3 (and an isothermal bottom wall), as the solid thickness ratio increases monotonically, the critical Rayleigh number decreases by 6.4%. Such a comparison at two dierent permeabilities suggests that as the system becomes more dense it is more susceptible to instability, a conclusion
that may appear to be counterintuitive at ®rst. However, it should be noted that Fig. 8 presents a plot of the critical Rayleigh±Darcy number as a function of the Darcy number. It shows that the medium will be less stable as its permeability approaches that of a liquid layer. Physically, the liquid layer presents less
Fig. 7. Critical wave number vs. solid thickness ratio for Da=10ÿ3.
3348
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
Fig. 8. Critical Rayleigh number vs. Darcy number a variety of boundary conditions (A=1).
resistance to the growth and development of an initial disturbance; whereas the porous medium acts to constrict such ¯ow initiation. Thus, as the Darcy number is increased towards the clear liquid limit, non-Darcian eects become more pervasive and the critical conditions deviate more from the predicted limiting values for the porous case. The Brinkman-extended Darcy model employed in this work allows for the comparison between the critical phenomena in a homogeneous liquid vis-aÁ-vis that in a saturated porous layer. Thus, an increase in the amount of solid in the system as well as in the Darcy number of the medium has a destabilizing in¯uence. The results in Figs. 5±7 indicate that the critical wave numbers tend to become decreased fractions of the critical wave number for the corresponding classical cases as the boundary condition approaches that of an insulated boundary, Theoretically, it is known that in the case of upper and lower insulated boundaries, a vanishing wave number may be assumed. The present results indicate that as the lower boundary becomes less conductive, it acts to constrain the critical wave number as expected from theory. However, in contrast, as the Darcy number increases, the critical wave number tends towards a value that represents a higher fraction of the classical result. Therefore, the largest critical wave number corresponds to a case with maximized Biot (constant temperature boundary) and Darcy numbers. In the BeÂnard case (pure liquid layer), which may be viewed as the limiting case of a porous
medium with an in®nite permeability (Daw0), the largest critical wave number corresponds to the limiting case of a constant temperature boundary and a minimum solid thickness (A=0). In other words, the presence of the porous medium and the solid act to constrain the size of the cell at the onset of convection. 5.2. Eect of ®xed boundary conditions In a typical physical situation the boundary during the process has ®xed heat transfer characteristics, usually that of an arbitrarily conducting boundary. By comparing Figs. 2±5, it is noted that decreasing the Biot number or moving towards a constant heat ¯ux boundary has a stabilizing eect on the system (the critical Rayleigh±Darcy number increases). Figs. 2±4 show the critical Rayleigh±Darcy number, Rac, dependence on the solid thickness for a variety of medium permeabilities for a ®xed thermal boundary condition. Fig. 8 also indicates that for all boundary conditions the system is more stable as the medium becomes more dense (Da 4 0), with the constant heat ¯ux boundary among these being the least stable. For Da 4 1, the permeability of the medium may be considered in®nite as the dierence in critical Rayleigh numbers for this case and the classical BeÂnard model is negligible. On the other hand, a signi®cant increase in Rc is observed for Da>10ÿ2. This introduces a transition region for values of Da where the medium is too dense to be treated as a bulk liquid and too porous to
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
model as a porous medium thus introducing nonDarcian behavior [10]. 5.3. Solidi®cation of an aluminum preform The linear stability analysis for the system in this work has two experimentally veri®able physical results: the critical temperature dierence for the onset of convection and the corresponding critical wavelength. This analysis has demonstrated, for a range of boundary conditions varying from a constant temperature to a constant heat ¯ux, how the incipient conditions for convection for a pure liquid solidifying in a homogeneous system may be determined. Speci®c calculations may be performed for a speci®c system of pure aluminum solidifying in an alumina preform in a gravitational ®eld, possessing the typical parameter values for such aluminum melt systems as [12], rL 2:39 103
kg , m3
mL 1:3 10ÿ3 Pa s, b 23:1 10ÿ6 Kÿ1 ,
kL 93
W , mK
DTL 208C,
aL 37 10ÿ6
m2 s
Tf 933:6 K
L 0:01 m: 47
3349
Fig. 9 indicates that it is possible to establish a relation between the critical Rayleigh numbers based on solid thickness ratio and the medium permeability as determined from the physical properties of the system to the modi®ed Lapwood solution from the theoretical analysis. It is shown that when the temperature dierence across the liquid layer is greater than 100 K, the aluminum system is inde®nitely unstable compared to the modi®ed Lapwood solution. If the system has a temperature dierence of 50 K, the system is unstable depending on the solid thickness ratio. To maintain a stable system other parameters would have to be adjusted to compensate for the eects of the temperature gradient. For the given properties, the system is always stable for a temperature dierence of less than 10 K.
6. Closing remarks The stability analysis performed for a solidifying liquid±saturated porous medium allows calculations of the critical Rayleigh and wave numbers for the incipient convective motion. In the asymptotic limits of A 4 0, Da 4 0, and BL 4 1, the critical Rayleigh±
Fig. 9. Critical conditions for the solidi®cation of aluminum for various temperature dierences across the liquid layer compared to critical conditions for the modi®ed Lapwood solutions.
3350
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
Darcy number (4p 2) and wave number (p ) of the classical Lapwood porous medium analysis are recovered from the model examined in this work. The presence of a solidi®cation boundary acts to destabilize the system. As the solid thickness increases the system becomes less stable as ascertained by a decrease in the Rayleigh number marking the onset of convection. Also, as the solid thickness is increased, the nondimensional wave number decreases, thus rendering smaller convection cells. The number of cells is directly related to the solidi®cation rate and hence the heat transfer rate. As the cell size decreases, a greater heat transfer rate should be expected. The largest wave number corresponds to the case with the highest Biot (constant temperature lower boundary) and Darcy numbers. Furthermore, as the lower boundary condition is modi®ed from an isothermal surface to an insulated boundary, the system becomes less stable. In comparison to the classical BeÂnard problem or that of solidi®cation of a pure liquid layer, the presence of the porous medium renders the porous Lapwood model more stable than the bulk liquid analog. As the system becomes more dense or as the permeability decreases, the system becomes more stable. Also it was found that the presence of the porous medium and the solidi®ed layer thickness act to constrain the size of the cell at the onset of convection.
References [1] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Dover, New York, 1961. [2] P.G. Drazin, W.H. Reid, Hydrodynamic Stability, Cambridge University Press, New York, 1981. [3] E.L. Koschmieder, in: BeÂnard Cells and Taylor Vortices, Cambridge University Press, New York, 1993, pp. 1±35. [4] E.R. Lapwood, Convection of a ¯uid in a porous medium, Proc. Camb. Phil. Soc. 44 (1948) 508±521. [5] A. Bejan, D.A. Nield, in: Convection in Porous Media, Springer-Verlag, New York, 1992, pp. 8±9, 305. [6] C. Karcher, U. MuÈller, Convection in a porous medium with solidi®cation, Fluid Dynamic Research 15 (1995) 25±26. [7] J.A. Weaver, R. Viskanta, Freezing of liquid-saturated porous media, ASME J. Heat Transfer 108 (1986) 654± 659. [8] C. Beckermann, R. Viskanta, Natural convection solid/ liquid phase change in porous media, Int. J. Heat Mass Transfer 31 (1) (1988) 35±46. [9] S.H. Davis, U. MuÈller, C. Dietsche, Pattern selection in single-component systems coupling BeÂnard convection and solidi®cation, J. Fluid Mech. 144 (1984) 133±151. [10] K.L. Walker, G.M. Homsy, A note on convective instabilities in Boussinesq ¯uids and porous media, ASME J. Heat Transfer 99 (1977) 338±339. [11] J.C. Ward, Turbulent ¯ow in porous media 1964, ASCE J. Hydraul. Div. 90 (1964) 6±12. [12] K.M. Fisher, The eects of ¯uid ¯ow on the solidi®cation of industrial castings and ingots, Physicochemical Hydrodynamics 2 (1981) 311±326.
Rayleigh±BeÂnard stability of a solidifying porous medium C. Mackie a, P. Desai b,*, C. Meyers c a
Tulane University, School of Mechanical Engineering, New Orleans, LA 70118, USA Georgia Institute of Technology, School of Mechanical Engineering, Atlanta, GA 30332, USA c North Carolina A&T, College of Engineering, Greensboro, NC 27411, USA
b
Received 17 February 1998; received in revised form 8 October 1998
Abstract In many solid/liquid phase transformation processes, natural convection controls the freezing or melting rate of the material. Freezing or melting in the presence of nonmelting components is important in numerous industrial processes and is also present in a wide range of systems in nature. The kinematics of the solid±liquid interface, coupled with bulk convection in the melt, play an important role in determining the microstructure of a solidi®ed material. The fundamental problem of thermoconvective instability of a single-component ¯uid in a horizontal porous layer has been studied extensively; additional complexities arise during solidi®cation due to the presence of a non-melting component, namely the porous medium. This study addresses the problem of Rayleigh±BeÂnard instability of a liquid layer undergoing a phase transformation within a porous medium. A linear stability analysis is performed to determine the eects of the medium and phase-change on the conditions for incipient instability under a range of thermal boundary conditions. The analysis reveals that the onset of convection, or the stability of the system, is signi®cantly aected by the presence of the porous medium, state of solidi®cation and the thermal boundary conditions. # 1999 Elsevier Science Ltd. All rights reserved.
1. Introduction The problem of thermoconvective instability in a horizontal ¯uid layer driven by buoyancy eects studied extensively in the literature [1±3] encompasses unresolved complexities associated with the freezing or melting of a liquid in the presence of a nonmelting component, such as a solid matrix or soil. This study addresses hydrodynamic stability issues during solidi®cation of such mixtures by modeling the matrix±liquid slurry as a liquid saturated horizontal porous medium in a gravitational ®eld going through a phase change under a full range of thermal conditions at the lower bounding surface. Freezing or melting in
* Corresponding author. Tel: +1-404-894-3233; fax: +1 404-894-7790. E-mail address: [email protected] (P. Desai)
the presence of nonmelting components is important in industrial applications as well as in a wide range of systems in nature. Examples include, but are not limited to, seasonal freezing and melting of soil, lakes and rivers, arti®cial freezing of the ground as a construction technique for supporting poor soils, insulation of underground buildings, the melting of the upper permafrost in the Arctic due to buried pipelines, thermal energy storage in porous media and storage of frozen foods. Metallurgical applications include manufacturing of composite materials and puri®cation of metals. The study of freezing and melting in porous media has gained a great interest over the past few years. A fundamental understanding of convective stability in such a range of applications is thus very important. It is recognized that ¯uid ¯ow, in general, and buoyancy-driven ¯ow, in particular, plays a critical role in determining interfacial traits of a solidifying material;
0017-9310/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 7 - 9 3 1 0 ( 9 8 ) 0 0 3 6 1 - 5
3338
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
Nomenclature a ac A BL ca ci cp D Da f g g hSL h1 kÃ, zà kL,S K L n, tÃi Pr Ra Rac (Rac) Rc (Rc) St T T1 T10 TM v w wà z
separation constant or wave number critical wave number solid thickness ratio Biot number, h1Z0/kL acceleration coecient general constants for characteristic equation speci®c heat at constant pressure dierential operator, D=d/dz, D 2=d2/dz 2 Darcy number, K/Z 20 temperature perturbation variable (liquid), see Eq. (33) temperature perturbation variable (solid), see Eq. (33) gravitational acceleration vector speci®c latent heat of solid heat transfer coecient of the bulk environment unit vector in the z-direction thermal conductivity for liquid and solid, respectively permeability of the porous medium length of the solidi®cation system unit normal and tangent vectors Prandtl number, nL/aL Rayleigh number, gb(TLÿTM)Z 30/nLaL critical Rayleigh±Darcy number, Rc Da critical Rayleigh number Stefan number, hSL/cpDTL temperature upper plate temperature lower plate temperature interface melting temperature velocity vector for liquid, v=v(u, v, w ) velocity component in the z-direction for the liquid vertical velocity perturbation variable spatial coordinate in vertical direction for liquid.
Greek symbols a ratio of the thermal diusivities of the liquid and the solid, aL/aS b volumetric coecient of thermal expansion acceleration coecient, ca/sH Pr ga DTL temperature dierence across the porous layer Z interface position coordinate Z normalized perturbation variable for the interface position Z0 length of the liquid saturated porous region time rate of change of the interface Zt m me/mL mL dynamic viscosity nL kinematic viscosity r ratio of the densities of the liquid and the solid, rL/rS s real part of the complex growth rate [s] sH heat capacity ratio, see Eq. (3) f medium porosity F pattern planform, where F=F(x, y ) H2 Laplacian operator H2H 2-D Laplacian operator, see Eq. (35).
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
Subscripts e L M S 0 1 Superscripts '
3339
eective property for porous medium liquid property or region melting temperature solid property or region reference state bulk environment. bar, basic state quantity perturbation variable normalized.
¯uid motion adjacent to a solidifying interface aects the local thermal and solutal ®elds which control the geometric, dynamic and thermodynamic characteristics of the interface. For instance, convection during fabrication or solidi®cation processing of metal matrix composites aects the distribution of particulates in the composites, in turn in¯uencing the properties of the material. In the pure metal analog, independent of the solute diusion, it is also known that the convective stability is aected by the density and presence of particulates as well as by the conductivity of the boundary. The porous medium model simulates the in¯uences of a non-interacting matrix, coupled with a range of thermal conditions at the lower boundary, on the initiation and growth of bulk convection in a solidifying nonhomogeneous system. This study focuses on the eect of the permeability of a porous medium and phase change on the incipient convection in the interstitial liquid; by parametrically changing the permeability from that corresponding to a Darcy model of the ¯ow through it to a vanishing one, an understanding is gained for the in¯uence of volume fraction and its distribution on the critical conditions for convective stability of the interstitial ¯uid in porous media in the presence of solidi®cation. This work couples solidi®cation dynamics with the thermal and hydrodynamic ®elds in a nonhomogeneous solidifying system via a Rayleigh±BeÂnard stability analysis. The porous medium is saturated with the liquid in all available pores, or the interparticle space of the medium, leaving no space for a vapor, gas, or other liquid within the medium Since the pioneering work on the porous-medium analog of the Rayleigh±BeÂnard stability problem [4], the subject of phase-change in a porous medium has received increasing attention over the past two decades. While liquid/gas phase change in a saturated porous medium has received considerable attention, solid/ liquid phase change in liquid-saturated porous medium
has received considerably less attention. Such problems and their manifestation represent a relatively new area in the ®eld of convection in porous media [5]. This study considers a horizontal, liquid saturated, porous medium undergoing solidi®cation while being cooled uniformly from above. A recent study [6] examined the coupled eects on stability arising from thermal convection and solidi®cation of a single-component liquid in a porous medium. The stability of the basic state of heat conduction and the stability of ®nite-amplitude convection in the Darcy limit of vanishing permeability for isothermal boundaries were examined via a twoparameter perturbation analysis. The present work investigates the eects of the solid matrix in a pure melt on the thermal convection for a range of permeabilities and the full range of lower wall thermal boundary conditions. The cited references are based on the Darcy ¯ow model, and hence are relevant to well packed low permeability porous media. This study diers by including the Brinkman extended ¯ow model that allows for the investigation of sparsely packed media. Beyond the context of convective stability, there exists substantial research in porous media transport phenomena germane to the present work. NonDarcian eects have been modeled illustrating the transition between the porous and clear ¯uid limits by modifying the Darcy equation with the Brinkman term [10]. Also experiments have been conducted using a modi®ed Darcy equation incorporating the Forchheimer term which attempts to account for microscopic drag forces due to acceleration [11]. It was demonstrated experimentally [7] that an eective porous medium thermal conductivity model is adequate for heat transfer estimates only when the solid matrix and pore-material have suciently similar thermal conductivities to allow for the requisite local thermal equilibrium. A continuum model of the porous medium was developed [8] for analyzing solid/liquid phase
3340
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
change in a porous medium that accounted for inertia as well as boundary friction eects. In particular, these include the Brinkman viscous and the Forchheimer pressure drag terms. Such results are relevant to the present work where the liquid and the porous material have similar conductivities, and where friction due to microscopic shear is also important.
2. The mathematical model A quiescent horizontal layer of a pure ¯uid in a gravitational ®eld heated from below does not always remain so in the presence of an adverse density gradient beyond a threshold value which is capable of giving rise to ¯uid motion. This observation may be applied to examine the solidi®cation (or melting) of a porous layer saturated with a pure liquid enclosed between two rigid, parallel, thermally dissimilar, horizontal plates of in®nite extent, a distance L apart, in a gravitational ®eld, as shown in Fig. 1. The analytical model is formulated under the assumptions that the liquid in the upper part of the layer is frozen and the eective thermophysical and transport properties of the porous medium are homogeneous and isotropic. The porous matrix and the solidi®ed layer are both assumed to be rigid [8]. Fluid motion within the por-
ous medium is described via a Newtonian, Boussinesq model of the Navier±Stokes equations, written for a liquid saturated porous medium. Motion due to density change upon phase change is not considered to be signi®cant. The liquid velocity through the medium is deemed suciently small, conforming to the Darcy regime, or, at most, to the Brinkman-modi®ed Darcy regime. The melt±freeze front is assumed to be a thin surface of negligible thickness that remains at the melting point of the phase-change material. The ¯uid layer is separated at the bottom from the surroundings by a thin rigid lamina of negligible heat capacity. The upper plate is maintained at a constant temperature, T1. 2.1. Governing equations In the porous layer, the heat transfer at the lower surface of the system is via a convective±conductive mechanism. The ¯ux at the bottom plate is maintained such that the melting temperature, TM, is bounded between the lower plate temperature, T10, and the upper plate temperature, T1 (T10>TM>T1). Therefore, there is a solid±liquid interface at z=Z0, 0 < Z0 < L. The governing equations are written in a Cartesian frame ®xed with respect to the composite system on the lower plate with the positive z-axis pointing
Fig. 1. Solidi®cation model with saturated porous medium.
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
upward (see Fig. 1). The temperature ®eld within the solidi®ed layer is described by the thermal energy equation @ TS aS r2 TS : @t
1
For the liquid, the temperature ®eld may be obtained from a simultaneous solution of the thermal energy and continuity equations, together with the Boussinesq model of the Darcy-extended momentum equation. These are expressed, respectively, as sH
@ TL v rTL aL r2 TL @t
2
where
2.2. Boundary conditions The boundary conditions at the phase boundary are the continuity of temperature TL TS TM
6
where the interface curvature due to undercooling of the melt is neglected, and the conservation of energy rS hSL
@Z j kS rTS ÿ kL rTL j n: @t
3
rv0
4
@Z kn @t
ca rL
@v ÿrP ÿ @t
mL v meff r2 v K
5
rL b TL ÿ Tref gk^ where ca is a constant which depends on the geometry of the porous medium and is determined by the nature of the assumed model pore tubes [5]. In most cases, the term becomes negligible when suitably scaled. This is due to the small value of permeability which multiplies the scaled inertia term. In the proceeding governing equations, the Brinkman-extended Darcy momentum model is employed. The frictional eects represented by Darcy law in the ¯uid momentum equation in a porous medium, are two orders of magnitude lower than those represented in the Navier±Stokes equations. Thus, only the impermeability condition at the surface is satis®ed, leaving the no-slip boundary condition unresolved. The Brinkman-extended Darcy model satis®es not only the impermeability and no-slip conditions but accounts also for the friction caused by microscopic shear by introducing an additional viscous force term via an eective dynamic viscosity, meH2v. The Brinkman-extended model reduces to a form of the Navier±Stokes model as the medium permeability approach in®nity (K 4 1) and to the conventional Darcy model as the medium permeability approach zero (K 4 0). Such a formulation allows a comparison of the ¯ows in porous media with those in clear ¯uids.
8
and the no-slip condition v t1 v t2 :
and
7
Also de®ned on this boundary (interface) are the conservation of mass (the kinematic condition) rL v n rL ÿ rS
f rcr L 1 ÿ f rcr S sH heat capacity ratio rcr L
3341
9
The boundary conditions on the two bounding surfaces are dTL at z 0, ÿ kL h1 T1 ÿ T10 dz 0 10 and at z L,
TS T1 :
3. Basic state and linear stability analysis Initially, the stationary system is in a static equilibrium state with no ¯uid motion and a planar interface at z=Z0. The purely conductive temperature pro®les in the liquid and the solid regions are, respectively, given by
TL
ÿh1 T1 ÿ TM z T1 ÿ TM h1 Z0 TM h1 Z0 kL h1 Z0 kL
11
TM ÿ T1 Z TM ÿ T1 zÿ 0 TM : Z0 ÿ L Z0 ÿ L
12
and
TS
The static equilibrium temperature pro®les are used to normalize the solid and liquid temperatures, respectively, as T L
TL ÿ TM ; T10 j ÿTM
T s
TS ÿ T1 : TM ÿ T1
13
3342
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
The scales for, v, x (and y, z ), t and p are chosen, respectively, as
distribution. In normalized form the basic state becomes (denoted by variables with bars on top)
aL sH Z20 ,Z0 , Z0 aL
v 0;
and
rL nL aL : Z20
The normalized equations become, after dropping the asterisks, for the solid, 1 @ TS aS 14 r 2 TS sH @ t aL and for the liquid, @ TL v rTL r2 TL , @t
rv0
15a,b
and @v ga rP @t
1 v ÿ m~ r2 v Ra TL k 0: Da
16
The normalized boundary conditions on the interface are v t1
1 ÿ r Zt k n, sH
r St Pr
@Z ArTL ÿ rTS n, @t
v1 t1 v2 t2 0
17a,b
and at z L,TS 1:
19a,b
It may be noted that the dimensionless group, BL, acts like a Biot number for heat transfer from the bulk environment to the liquid. As BL approaches zero, the boundary condition at the surface approaches that of an insulated surface; as BL approaches in®nity, the boundary condition approaches that of a constant surface temperature; also appearing is another dimensionless parameter: A
L ÿ Z0 kS DTS Z0 kL DTL
T L 1 ÿ z;
zÿ1 : T S 1 ÿ Lÿ1
22
It is noted that TL and TS are scaled by dierent reference lengths. Thus at z=1, the temperature for the solid and liquid regions (TL and TS) are equal to the melting temperature, TM [Eq. (12)]. Next, a linear perturbation model of the system is developed by considering a slight change in temperature due to a slight melting of the solid. Such a problem is formulated by introducing perturbations of the basic state de®ned via a conventional barred-primed scheme of notations as TS T S T 0 S ,
TL T L T 0 L
p p p 0 ,
23
Z 1 Z 0:
24
18a,b
The surface boundary conditions become dTL BL TL , dz
21
and
v v 0,
and TL TS 0:
at z 0,
dp Ra T L k dz
20
where z=Z0 is the position of the planar interface and z=L is the thickness of the combined solid±liquid system. For a small A, A ÿ1 may be considered the equivalent Biot number for heat transfer from the liquid to the solid; also, A measures the amount of solid present in the system [9]. The basic state under consideration is that of a stagnant layer of liquid, with a hydrostatic pressure distribution, and a purely conductive temperature
Introduction of the perturbation expressions, Eqs. (23) and (24), into the governing Eqs. (14)±(16), followed by linearization in the disturbance quantities, yields the perturbation model for the solid, 0 1 @T S aS r2 T 0 S 25 sH @t aL and for the liquid @ T 0L w 0 r2 T 0 L , @t
r v0 0
26
and ga
@v0 rp 0 @t
1 v 0 ÿ m~ r2 v 0 Ra T 0 L k 0 Da
27
where w' is the vertical component of the perturbation velocity v '. The linearized boundary conditions for the perturbation model are at z 1,
T 0 S ÿZ 0
@ T S , @z
T 0 L ÿZ 0
0 @Z @ T 0S @ T 0L r St Pr A ÿ @t @z @z
@ T L @z
28
29
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
and u 0 1 ÿ r
@Z @x
0
,
@Z 0 v 0 1 ÿ r , @y
@Z w 0 1 ÿ r : @t
Substitution of the normal mode into the linear perturbation model [Eq. (32)] results in the normal-mode disturbance equations for the solid 30
at z 0, T
0
S
and at z L,
31
0:
@ 2 1 2 r w ÿ r w r2 r2 w Ra r2H TL 0 @t Da
@ TL ÿ w ÿ r 2 TL 0 @t r2 TS 0
aS 2 D ÿ a2 g aL
36
sf w^ D2 ÿ a2 f
The solution to the preceding set of equations results in a sucient condition to ascertain the stability boundary of the liquid layer, as would be expected from the linear analysis. The linear perturbation system, [Eqs. (25)±(31)], is examined to determine the critical and necessary conditions for incipient nonlinear convection. After eliminating the pressure from the liquid momentum equation [Eq. (27)] by taking its curl, the curl is taken again, yielding
ÿga
sg
and, for the liquid
The surface boundary conditions become dT 0 L BL T 0 L , dz
3343
(32)
37
and 2
2
2
2 2
ÿga s D ÿ a w^
1 D2 ÿ a2 w^ Da
38
2
ÿ D ÿ a w^ Ra a f where an ordinary dierential operator, D, has been de®ned as D
d dz
and D2
d2 : dz2
39
The boundary conditions transform as g L 0,
g 1
Z^ A
40
for the solid region, and ^ Df 1 ADg 1 w 1 1 ÿ rsZ
^ Dw 1 1 ÿ rsZ
^ ^ w 0 0 Dw 0 0 Df 0 BL f 0 where the primes have been dropped. Next, to ®nd normal modes solution to the linear stability problem it is noted that the variables become separable under this assumption, yielding solutions of the form ^ Zg est F x,y fTS ,TL ,w,Zg f g z,f z,w,^
33
where the time constant, s=sr+isi, contains sr, the growth rate, and si, the frequency of the disturbance, and F is the plan form function which determines the cellular structure of the ¯uid motion and satis®es the two-dimensional wave or membrane equation Fxx Fyy ÿa2 F
34
or r2H F
2
a F 0,
where
r2H
@ @ : @ x 2 @ y2
35
It may be noted that `a' is the separation constant, mathematically, in addition to being the nondimensional overall wave number, physically.
41 in the liquid region.
4. Solutions 4.1. Solid region solution It is noticed that the general solution for the solidi®ed layer may be obtained independent of the liquid equations, resulting in g z C1 sinh az C2 cosh az:
42
Use of the boundary conditions, [Eq. (41)], yields the temperature pro®le in the solid as g z
ÿZ sinh a 1 A ÿ z : A sinh aA
43
It was necessary ®rst to solve for the solid region, independent of the liquid region, in order to express all of the governing equations in the liquid region in terms of the temperature perturbation variable, f.
3344
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
4.2. Liquid region solution The ®rst two (liquid momentum and energy) equations in Eq. (32) are combined, thus eliminating the vertical velocity component, forming a single-variable perturbation equation in terms of the normal mode perturbation temperature, f. The combined expression becomes 1 1 2 2 2 ÿ ga s ÿ s D ÿ a f ga s m~ s Da Da 44 D2 ÿ a2 2 f ÿ m~ D2 ÿ a2 3 f Ra a2 f 0: Next, it is assumed that the principle of exchange of stability is valid, or that all possible values of the growth rate allowable are real. The implication of s being real is that the neutral or marginal stability curve is characterized by a growth rate, s, of zero. Accordingly, the single-variable temperature perturbation equation for the liquid [Eq. (44)] reduces to a sixth-order, linear, homogeneous equation 1 2 2 3 D ÿ a f ÿ D2 ÿ a2 2 f Ra a2 f 0 45 Da with six linear, homogeneous boundary conditions (determined by applying the solid solution) and evaluating the continuity and thermal energy equations at z=(0, 1), given as D2 ÿ a2 f 0 0, D D2 ÿ a2 f 0 0, Df 0 BL f 0
D2 ÿ a2 f 1 0 D D2 ÿ a2 f 1 0
and Df 1 a Coth aAf 1 0: 46
The linear, homogeneous, sixth-order system [Eqs. (45) and (46)] constitutes an eigenvalue problem. The system obviously has the trivial solution, f=0, for any values of a, BL, Ra, and Da, but, in general, also possesses non-trivial solutions depending on an eigenvalue, say l, and has a general solution obtained as a linear combination of six linearly independent solutions. In this case, a, BL, and Da become ®xed parameters, with Ra serving as the eigenvalue, producing non-trivial continuous solutions for the perturbation temperature eigenfunction, f(z ). Use of the six boundary conditions and the general solution, f(z ), yields six simultaneous linear equations. This de®nes the eigenvalue problem in which Ra serves as the eigenvalue for a particular set of values for a, BL, Da and A. The critical Rayleigh number for the onset of convection has a corresponding critical wave number, both of which being
the minimum values satisfying each value as an eigenvalue. The solution for the stability problem is obtained for given values of BL, Da and A, by determining iteratively the lowest characteristic values of a and Ra. 5. Results and discussion The linear stability analysis for the system in this work has two experimentally veri®able results: the critical temperature dierence for the onset of convection and the corresponding critical wavelength. For a range of boundary conditions varying from a constant temperature to a constant heat ¯ux, the incipient conditions for convection for a pure liquid solidifying in a porous medium are determined. 5.1. Eect of varying solid thickness and boundary conditions Since the solid portion of the porous medium is typically held at temperatures below the melt temperature, initial freezing begins with the contact between the medium and the superheated liquid. Physically, the removal of the latent heat takes place rapidly, creating a thin coating on the porous medium, the degree to which the permeability is aected due to such a phenomenon being quite small. Therefore, any change in permeability will be due to the solidi®cation of the interstitial ¯uid within the slightly coated matrix. In the model, cooled from above, and/or heated from below, solidi®cation occurs only at the interface; interstitial solidi®cation does not occur. Therefore, Figs. 2±8, the critical conditions for incipient convection at ®xed permeabilities are noted. As the solidi®ed portion of the system grows, or as, the aspect ratio of the layer, A, increases for a set of ®xed values of the Darcy number, Da, and the Biot number, BL, the critical Rayleigh number of the system decreases, rendering the system less stable. It may be noted, for a system with no solidi®cation (A=0), constant temperature boundary, BL=1000, and a dense medium, DaE10ÿ6, that the classical problem of convection in a porous medium between two in®nite walls, the so-called `Lapwood problem', and the solution of Rac=4p 2=39.45 are recovered. Departure from the classical problem of rigid boundaries was determined as a function of solid thickness and the results are designated on Figs. 2±7 as the Lapwood problem recovering the results of Karcher et al. Results for other conditions may be compared with those resulting from the modi®ed classical Lapwood case. Figs. 2±4 illustrate the eect of the solid thickness ratio on the critical Rayleigh number for ®xed permeability for dierent boundary conditions, compared to the modi-
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
3345
Fig. 2. Critical Rayleigh±Darcy number vs. solid thickness ratio for Da=10ÿ6.
®ed Lapwood solution, which is denoted by a thin dashed line. The ®gures also show that as the boundary condition at the bottom wall is modi®ed from that of a ®xed temperature (BL 4 1), to that of an insulated boundary (BL 4 0), the system becomes less stable. For all values of the Darcy number considered, simultaneously increasing both the solid thickness ratio and the Biot number has a strictly stabilizing in¯uence. Figs. 2±7 indicate further that the critical Rayleigh±
Darcy and wave numbers approach their lower limits asymptotically as A 4 1. Such a destabilizing eect was previously documented for solidi®cation in a bulk liquid-layer [9] as well as in a dense liquid saturated porous medium [6]. Considering the limiting case of an isothermal lower wall (BL=1000) and a densely packed medium (Da=10ÿ6), there is a 16.44% decrease in the critical Rayleigh±Darcy number as the (normalized) solid thickness, A, increases monotonically from 0±1;
Fig. 3. Critical Rayleigh±Darcy number vs. solid thickness ratio for Da=10ÿ4.
3346
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
Fig. 4. Critical Rayleigh±Darcy number vs. solid thickness ratio for Da=10ÿ3.
Fig. 5. Critical wave number vs. solid thickness ratio for Da=10ÿ6.
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
3347
Fig. 6. Critical wave number vs. solid thickness ratio for Da=10ÿ4.
for the less dense case of Da=10ÿ3 (and an isothermal bottom wall), as the solid thickness ratio increases monotonically, the critical Rayleigh number decreases by 6.4%. Such a comparison at two dierent permeabilities suggests that as the system becomes more dense it is more susceptible to instability, a conclusion
that may appear to be counterintuitive at ®rst. However, it should be noted that Fig. 8 presents a plot of the critical Rayleigh±Darcy number as a function of the Darcy number. It shows that the medium will be less stable as its permeability approaches that of a liquid layer. Physically, the liquid layer presents less
Fig. 7. Critical wave number vs. solid thickness ratio for Da=10ÿ3.
3348
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
Fig. 8. Critical Rayleigh number vs. Darcy number a variety of boundary conditions (A=1).
resistance to the growth and development of an initial disturbance; whereas the porous medium acts to constrict such ¯ow initiation. Thus, as the Darcy number is increased towards the clear liquid limit, non-Darcian eects become more pervasive and the critical conditions deviate more from the predicted limiting values for the porous case. The Brinkman-extended Darcy model employed in this work allows for the comparison between the critical phenomena in a homogeneous liquid vis-aÁ-vis that in a saturated porous layer. Thus, an increase in the amount of solid in the system as well as in the Darcy number of the medium has a destabilizing in¯uence. The results in Figs. 5±7 indicate that the critical wave numbers tend to become decreased fractions of the critical wave number for the corresponding classical cases as the boundary condition approaches that of an insulated boundary, Theoretically, it is known that in the case of upper and lower insulated boundaries, a vanishing wave number may be assumed. The present results indicate that as the lower boundary becomes less conductive, it acts to constrain the critical wave number as expected from theory. However, in contrast, as the Darcy number increases, the critical wave number tends towards a value that represents a higher fraction of the classical result. Therefore, the largest critical wave number corresponds to a case with maximized Biot (constant temperature boundary) and Darcy numbers. In the BeÂnard case (pure liquid layer), which may be viewed as the limiting case of a porous
medium with an in®nite permeability (Daw0), the largest critical wave number corresponds to the limiting case of a constant temperature boundary and a minimum solid thickness (A=0). In other words, the presence of the porous medium and the solid act to constrain the size of the cell at the onset of convection. 5.2. Eect of ®xed boundary conditions In a typical physical situation the boundary during the process has ®xed heat transfer characteristics, usually that of an arbitrarily conducting boundary. By comparing Figs. 2±5, it is noted that decreasing the Biot number or moving towards a constant heat ¯ux boundary has a stabilizing eect on the system (the critical Rayleigh±Darcy number increases). Figs. 2±4 show the critical Rayleigh±Darcy number, Rac, dependence on the solid thickness for a variety of medium permeabilities for a ®xed thermal boundary condition. Fig. 8 also indicates that for all boundary conditions the system is more stable as the medium becomes more dense (Da 4 0), with the constant heat ¯ux boundary among these being the least stable. For Da 4 1, the permeability of the medium may be considered in®nite as the dierence in critical Rayleigh numbers for this case and the classical BeÂnard model is negligible. On the other hand, a signi®cant increase in Rc is observed for Da>10ÿ2. This introduces a transition region for values of Da where the medium is too dense to be treated as a bulk liquid and too porous to
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
model as a porous medium thus introducing nonDarcian behavior [10]. 5.3. Solidi®cation of an aluminum preform The linear stability analysis for the system in this work has two experimentally veri®able physical results: the critical temperature dierence for the onset of convection and the corresponding critical wavelength. This analysis has demonstrated, for a range of boundary conditions varying from a constant temperature to a constant heat ¯ux, how the incipient conditions for convection for a pure liquid solidifying in a homogeneous system may be determined. Speci®c calculations may be performed for a speci®c system of pure aluminum solidifying in an alumina preform in a gravitational ®eld, possessing the typical parameter values for such aluminum melt systems as [12], rL 2:39 103
kg , m3
mL 1:3 10ÿ3 Pa s, b 23:1 10ÿ6 Kÿ1 ,
kL 93
W , mK
DTL 208C,
aL 37 10ÿ6
m2 s
Tf 933:6 K
L 0:01 m: 47
3349
Fig. 9 indicates that it is possible to establish a relation between the critical Rayleigh numbers based on solid thickness ratio and the medium permeability as determined from the physical properties of the system to the modi®ed Lapwood solution from the theoretical analysis. It is shown that when the temperature dierence across the liquid layer is greater than 100 K, the aluminum system is inde®nitely unstable compared to the modi®ed Lapwood solution. If the system has a temperature dierence of 50 K, the system is unstable depending on the solid thickness ratio. To maintain a stable system other parameters would have to be adjusted to compensate for the eects of the temperature gradient. For the given properties, the system is always stable for a temperature dierence of less than 10 K.
6. Closing remarks The stability analysis performed for a solidifying liquid±saturated porous medium allows calculations of the critical Rayleigh and wave numbers for the incipient convective motion. In the asymptotic limits of A 4 0, Da 4 0, and BL 4 1, the critical Rayleigh±
Fig. 9. Critical conditions for the solidi®cation of aluminum for various temperature dierences across the liquid layer compared to critical conditions for the modi®ed Lapwood solutions.
3350
C.Mackieet al. / Int. J. Heat Mass Transfer 42 (1999) 3337±3350
Darcy number (4p 2) and wave number (p ) of the classical Lapwood porous medium analysis are recovered from the model examined in this work. The presence of a solidi®cation boundary acts to destabilize the system. As the solid thickness increases the system becomes less stable as ascertained by a decrease in the Rayleigh number marking the onset of convection. Also, as the solid thickness is increased, the nondimensional wave number decreases, thus rendering smaller convection cells. The number of cells is directly related to the solidi®cation rate and hence the heat transfer rate. As the cell size decreases, a greater heat transfer rate should be expected. The largest wave number corresponds to the case with the highest Biot (constant temperature lower boundary) and Darcy numbers. Furthermore, as the lower boundary condition is modi®ed from an isothermal surface to an insulated boundary, the system becomes less stable. In comparison to the classical BeÂnard problem or that of solidi®cation of a pure liquid layer, the presence of the porous medium renders the porous Lapwood model more stable than the bulk liquid analog. As the system becomes more dense or as the permeability decreases, the system becomes more stable. Also it was found that the presence of the porous medium and the solidi®ed layer thickness act to constrain the size of the cell at the onset of convection.
References [1] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Dover, New York, 1961. [2] P.G. Drazin, W.H. Reid, Hydrodynamic Stability, Cambridge University Press, New York, 1981. [3] E.L. Koschmieder, in: BeÂnard Cells and Taylor Vortices, Cambridge University Press, New York, 1993, pp. 1±35. [4] E.R. Lapwood, Convection of a ¯uid in a porous medium, Proc. Camb. Phil. Soc. 44 (1948) 508±521. [5] A. Bejan, D.A. Nield, in: Convection in Porous Media, Springer-Verlag, New York, 1992, pp. 8±9, 305. [6] C. Karcher, U. MuÈller, Convection in a porous medium with solidi®cation, Fluid Dynamic Research 15 (1995) 25±26. [7] J.A. Weaver, R. Viskanta, Freezing of liquid-saturated porous media, ASME J. Heat Transfer 108 (1986) 654± 659. [8] C. Beckermann, R. Viskanta, Natural convection solid/ liquid phase change in porous media, Int. J. Heat Mass Transfer 31 (1) (1988) 35±46. [9] S.H. Davis, U. MuÈller, C. Dietsche, Pattern selection in single-component systems coupling BeÂnard convection and solidi®cation, J. Fluid Mech. 144 (1984) 133±151. [10] K.L. Walker, G.M. Homsy, A note on convective instabilities in Boussinesq ¯uids and porous media, ASME J. Heat Transfer 99 (1977) 338±339. [11] J.C. Ward, Turbulent ¯ow in porous media 1964, ASCE J. Hydraul. Div. 90 (1964) 6±12. [12] K.M. Fisher, The eects of ¯uid ¯ow on the solidi®cation of industrial castings and ingots, Physicochemical Hydrodynamics 2 (1981) 311±326.