Melting of ice-aluminum balls system

Melting of Ice-Aluminum Balls System []

I

S. Cheil~ah R. Viskanta Heat Transfer Laboratory, School of Mechanical Engineering, Purdue University, West Lafayette, Indiana

Illlll

mlll l

illl

I

IIMelting of ice-porous media (aluminum balls) system contained in a rectangular test cell has been studied both experimentally and numerically. When the superheat across the liquid region was small, the flow in the porous media was weak and the interface was almost planar. For larger superheats, natural convection flow, the interface velocity, and shape were found to depend on the imposed temperature and the permeability of the porous medium. The measured temperature distributions were compared with the predictions of a numerical model that considered both conduction in the solid and natural convection in the liquid region. The model is based on volumetric averaging of the macroscopic transport equations, with phase change assumed to occur volumetrically over a small temperature range. Both Brinkrnan and Forchheimer extensions were added to the Darcy equations. The effect of density inversion of water on the fluid flow and heat transfer has been modeled. Reasonably good agreement has been found between the experimental data and numerical predictions. Heat conduction in the high thermal conductivity porous matrix was found to influence considerably the melting front shape and motion.

Keywords: porous media, convection, heat conduction

INTRODUCTION Solid/liquid phase change in saturated porous media occurs in a wide variety of systems in nature and engineering. Applications include freezing and thawing of soils [1], artificial freezing of ground as a structural support and as a water barrier for construction and mining purposes [2], excavation of frozen soil [3], latent heat-of-fusion energy storage [4], heat transfer in soil around the heat exchanger coils of a ground-based heat pump [5, 6], and metallurgy [7]; other applications have been discussed at a recent International Symposium on Cold Region Heat Transfer [8]. Despite these varied and many applications, relatively little attention has been given to the study of solid/liquid phase change of liquid saturated porous media [9]. The related problem of natural convection in porous media in the absence of phase change has been investigated both experimentally and numerically by many researchers, and reviews are available [10, 11]. There have been a number of theoretical [2, 12-17] and experimental [ 15-17] studies concerned with melting of frozen porous media. Goldstein and Reid [12] studied theoretically thawing of frozen porous media using a mapping technique based on the theory of complex variables by solving the energy equation in the frozen region without knowing the shape of the region. Thawing of soil using electrical heating [13] and microwave energy [3] has been investigated. Both finite element and boundary integral methods were used to model thawing of frozen soil in the absence of convection [14]. Melting of porous media-ice system contained in a cylindrical capsule was studied using a heat-conduction-based model [15]. The

phase-change phenomena,

natural

melting of ice around a horizontal tube embedded in a porous medium (glass beads) has been investigated by Okada and Fukarnoto [16]. Melting of a solidified glass beads-gallium system from a vertical wall has been studied numerically employing a two-dimensional enthalpy-based model [17]. It was conclusively established both theoretically [16, 17] and experimentally [15-17] that natural convection in the liquid region considerably influences the shape and the motion of the melting front. This article reports on an experimental and theoretical melting study of an aluminum balls-ice system. Beckermann and Viskanta [17] have found that when the effective thermal conductivity of the phase change medium is large, heat conduction in the solidified region is significant and affects the melting front shape and motion. Therefore, it is of interest to study the opposite case when the porous matrix has a large thermal conductivity and the phase change medium (ice-water) has a much smaller thermal conductivity. For the system considered, the thermal conductivity ratio is about 200. In addition, water has a maximum density at about 4°C, and the effect of the density inversion of water on the melting process in the presence of natural convection in the liquid region may be significant and needs to be accounted for. The macroscopic transport equations, volumetrically averaged, with phase change assumed to occur volumetrically over a small temperature range [17], are used. In the liquid region, both the Brinkman and Forchheimer extensions were included in the Darcy's equation. The effect of density inversion of water on the fluid flow and heat transfer were also modeled. The experimental temperature data and the deduced solid/liquid interface positions were compared with the predictions of the numerical model.

Address correspondence to Professor R. Viskanta, Heat Transfer Laboratory, School of Mechanical Engineering, Purdue University, West Lafayette, IN 47907.

Experimental Thermal and Fluid Science 1990; 3:222-231 ©1990 by Elsevier Science Publishing Co., Inc., 655 Avenue of the Americas, New York, NY 10010

222

0894-1777/90/$3.50

Melting of Ice-Aluminum Balls System 223

,///////,

E~ERIMENTS

,,///,

Test Cell Melting experiments were performed in a rectangular test cell with inner dimensions of 205 mm in length, 203 mm in height, and 127 mm in width. The top, bottom, front, and back sides were made of Plexiglas (12.5 mm thick). Two Plexiglas plates separated by a 6 mm air gap were used on the front and back of the test cell to minimize the effect of heat gains and condensation of moisture. A 7.5 mm wide and 172 mm long slot was cut in the top plate. The test cell was filled with water and aluminum balls through this opening. A lid fit snugly into this slot. Two 11 mm diameter holes were cut on this lid to bring out the thermocouple wires. The entire test cell was covered with 50 mm thick Styrofoam on all sides. The test cell was placed on an iron plate fitted with leveling screws. Two copper heat exchangers with milled passages for flow of coolant constituted the left and right end-walls of the test cell. The flow passages inside the heat exchanger were milled in such a way that the maximum temperature variation along the surface of the heat exchanger was within :L-0.2°C. Six thermocouples were placed along the surface of the heat exchanger to check the uniformity of the temperature along the faces. Because of the high thermal conductivity of the ice-aluminum system, the initial temperature was uniform to within 0.2°C. Two thermocouple rakes with 21 copper-constantan thermocouples (0.8 mm apart) in each and supported on two half rings epoxyed onto the surface of the heat exchangers were located at ~1 and 2~ the height from the bottom and along the centerplane of the test cell. This was a compromise between the resolution desired and the disturbance of the system by additional rakes. The temperature readings were automatically recorded, at preset intervals, using a data logger connected to a VAX microcomputer.

Test Materials and Procedure Aluminum balls having diameters of 3.2, 4.8, and 12.2 mm constituted the porous media. The properties used for the numerical study were for with a chemical composition as close as could be obtained to the chemical composition of the beads used in this study [18]. The test cell was filled with aluminum balls with rakes kept in position. Once distilled and degasified water was carefully siphoned into the test cell without introducing air bubbles into the system, the entire test cell was placed inside the freezer compartment of a refrigerator to freeze the water. A mixture of ethyl alcohol and water was circulated through the heat exchangers from two constanttemperature baths, and the porous matrix and ice were allowed to warm to the desired initial temperature. The melting was initiated by switching one of the constant temperature baths to a third one that was already heated separately to the desired temperature. The natural convection flow in the unfrozen porous media and the flow structure were observed in parallel flow visualization experiments. ANALYSIS The physical system modeled consists of a rectangular cavity with two vertical walls maintained at two different temperatures and the top and bottom surfaces insulated. The cavity is filled with a mixture of water and uniform aluminum balls,

(a) T

,'////,/,,,,'/

ky,t/

~/(//

////.,

,

,/

(b) Figure 1. Schematic of test apparatus, top view (a) and fromview (b): (1) heat exchanger, (2) Plexiglas wall, (3) air gap, (4) phase change region (porous medium-water-ice), (5) thermocoupie rake, (6) Styrofoam insulation, and (7) constant temperature baths. frozen, and then at time, t > 0, a uniform temperature greater than the fusion temperature is imposed on the left wall. Melting is initiated at this wall, and the interface moves from left to right (Fig. 1). The following simplifying assumptions are made in the analysis: (1) the porous medium is isotropic, homogeneous, and has uniform porosity; (2) the porous matrix and the phase change material (PCM) are in local thermal equilibrium; (3) the flow is two-dimensional, laminar, and incompressible; (4) the volume change caused by freezing is negligible; (5) the phase change medium has a definite fusion temperature; (6) the porous matrix and the solid are stationary; (7) all thermophysical properties are independent of temperature; and (8) the momentum flux caused by the interface motion is negligible. It is recognized that the porosity is not uniform (discussed in a later section) in the test cell and that ice contracts (approximately 8%) on melting. But the consideration of these factors will complicate the solution procedure of the already difficult problem. The other assumptions are standard and realistic for the melting process. The complicated geometry, heat transfer, and flow pattern in the system prohibit the solution of the governing equations at the microscopic level. Hence, the equations are averaged over a small volume element. In general, in such a volume element, there can be both solid and liquid PCM along with the porous matrix. In the melt region, the PCM is entirely liquid and 3' = 1, ~ -- ~b, and, in the solid region, 3' = 6 = 0. Fuller details of the development of the volume-averaged governing equations can be found in the literature [10, 11]. Even though the phase change is assumed to occur at a discrete temperature, in a volume element containing both the porous matrix as well as the solid and liquid, the average temperature may be slightly higher or lower than Ty. Hence, it is assumed that

22,1 S. Chellaiah and R. Viskanta both solid and liquid may exist simultaneously in a volume element, if its temperature lies with a small temperature difference AT, on either side of the fusion temperature. A more detailed discussion of the model is given elsewhere [ 17]. The two-dimensional conservation of mass, momentum, and energy equations for freezing of liquid saturated porous media are, respectively, ~. O = 0 m 00 + ~p~( 0 . O-~-

(l)

v)O = - V p + ~ v : O

dictions and measured temperatures and interface location for melting of ice-aluminum balls system initially at saturation temperature [15]. The effective thermal conductivity of the liquid region was taken to be the value determined experimentally by Lindemarm [23]. It should be stressed that the effective thermal conductivity of liquid saturated porous media under static conditions was assumed to be the same as that for a transient system with a flowing fluid in which dispersion effects may be present. The value of the permeability, K, was calculated from the Kozeny-Carman equation,

o~C

d2~b 3

K -OT

oc-~- + ow~(O. VT) = V. (ke~Vr) - 4,olZ~f-~t (3) The momentum balance takes into account the unsteady term, Forcl~eimer's extension, and Brinkman's extensions to the Darcy equations. Equation (3) is the volume-averaged energy equation for the general control volume containing the porous medium solid-liquid mixture. With the velocity set to zero ( 0 = 0), Eq. (3) is also appropriate for the frozen region. The boundary conditions for temperature are: T=Th

atx=0forally

(4)

T = T c a t x = L for all y

0T -- =0aty=0andy=Hforallx

Oy

(5)

There is no slip at the walls. At the free surface the boundary conditions for the velocity are: Ou

-- =v =0aty=Hfor0
oy

175(1 - ~b)2

(9)

The value of the inertia coefficient C in Forchheimer's extension was found to be nearly constant [24], and a value of 0.55 is used in the calculations. The model equations were solved with use of the SIMPLER algorithm [25]. Computational details are given in an earlier article [17]. After conduction of numerical sensitivity studies with different grids and time steps, a uniform grid of 26 × 26 nodal points was chosen as a compromise between cost and accuracy. Extensive grid independence studies could not be conducted because of the high cost and large computer time requirements. A dimensionless, time step of z = 3.13 × 10 -5 (t = 1 s) was utilized when the transients are large, and the step was increased to r = 1.57 × 10 -5 (t = 5 s) at later times in the simulation. After some preliminary investigation with different values of A0, a value of A0 = 0.075 was chosen. The simulations were terminated when the significant variables (velocities and temperatures) agreed up to three decimal places and the residual mass was less than 10 - 6 .

(6) RESULTS AND DISCUSSION

Fort <0,

Experimental Results

u = v = O and T = Tc.

The buoyancy force is og', with p denoting the local density, corresponding to the local temperature. For fluids having a linear density-temperature relationship, the usual simplification of the buoyancy term can be done. But the density-temperature relationship for water is nonlinear, and it attains a maximum value at 3.98°C. Several equations of state (second, third, and fourth degree polynomials) [19-21] have been proposed for the density of water as a function of temperature. In this study, we use the approximation for density of water suggested by G-ebhart and Mollendorf [20], p = pro(1 -o~ I(T - T i n ] q)

(7)

where ~0 = 9 . 2972 × 1 0 - 6 ( ° C ) - q , q = 1.8948. The definition of the mean thermal capacitance of the mixture 05c) is given by = 4~[3'01ct + (1 - "y)OsCs] + (1 - 4,)ppCp

(8)

There is considerable uncertainty in the choice of a model for predicting the effective thermal conductivity of an ice-aluminum system. Therefore, the effective conductivity of the frozen porous media (solid region) was calculated using the Veinberg model [22], for it was found to give the best comparison (as compared to different models) between pre-

A number of experiments with the several different size beads and several different superheats were conducted. For space limitations, only few of them are discussed here. Additional results are presented elsewhere [ 18]. The experimental conditions are summarized in Table 1. All properties of ice, water, and aluminum balls were taken at the fusion temperature 0°C. Since water undergoes a density inversion, the Rayleigh number was calculated as suggested in the literature [21]. From Table 1, it is seen that the porosity ~ varies with ball size. For a system of infinite volume, randomly packed with uniform size spheres, the porosity is constant and is independent of the ball size [26]. Since the test cells used in experiments are finite in size, there is a considerable variation of porosity near the walls, especially with the larger-size balls. With 3.2 mm balls, the porosity is close to the theoretically expected value for randomly packed beds. The porosity of the system with the two smaller-size balls are nearly the same, because there is not much difference in their size to affect the packing density. But, with the larger-size balls, the porosity is higher. There are a large number of parameters that affect the transport processes during melting of frozen porous media. The parameters include the Stefan number (Ste), the Darcy number (Da), the Rayleigh number of a porous medium (Ra*), the subcooling parameter (S), and the aspect ratio (A). Some of

Melting of Ice-Aluminum Balls System 225 Table 1. Summary of Experimental Conditions for Melting of Ice-Aluminum Balls System (,4 = 0.75)

Exp.

d (mm)

Th (*C)

Tc (*C)

~

K x 109 (m2)

Da x 107

S

Ste

Ra*

1 2 3 4 5 6 7

12.2 12.2 12.2 12.2 12.2 4.8 3.2

4.4 8.4 19.3 22.3 20.3 21.8 26.2

-0.2 -0.2 -0.2 -0.2 -6.4 - 0.2 -0.2

0.431 0.431 0.431 0.431 0.431 0.398 0.395

210.3 210.3 210.3 210.3 210.3 22.9 9.9

50.04 50.04 50.04 50.04 50.04 54.44 2.34

0.0 0.0 0.0 0.0 0.04 0.0 0.0

0.061 0.107 0.247 0.282 0.257 0.276 0.320

328.6 945.7 4629.8 5962.0 5010.9 621.9 383.4

these parameters are not independent of each other. For example, the temperature difference (Th - T / ) appears in both the definitions of the Stefan and the Rayleigh numbers. As a consequence, it is difficult to be comprehensive. Some typical results are presented in the article, and additional results are available elsewhere [18]. Effect o f Ball S i z e

The effect of ball size can be investigated by comparing the time-temperature traces and time-interface location traces for experiments 4, 6, and 7. It should be noted that the subcooling parameter for the three experiments is the same (S = 0), while the Stefan number for experiment 7 is slightly higher than that for experiments 4 and 6 (Table 1). The temperature distribution at different times along the two rakes for experiment 4 are shown in Fig. 2. The initial subcooling is - 0 . 2 ° C and is negligible. The temperaturetime traces reveal that the melting along the top rake is faster than that along the bottom. This is same as that observed in the melting of pure PCMs in the absence of a porous matrix [27]. The effect of thermal convection can be seen as early as t = 0.5 h, when the melt thicknesses at the top and bottom are 7.1 and 5.5 cm, respectively. At this time, the temperature distribution along the top and bottom rakes is characterized by two regions of sharp drop. Near the top of the hot wall, there is a sudden drop of 5°C in a distance of 1.45 cm, and along the bottom rake the corresponding value is 10°C. Near the solid-liquid interface, there is a sharp decrease of 15°C at the bottom in a distance of 1.2 cm of PCM. As the experiment

25[

7

T

on[ ~[-

II

0

V tCh> l O

~ 3

0.$ ~

4

5

i0

(a)

7

1

15

T

V


0.~5 k

1 2

O

t.O

4

t.0

3

I\\'~

0

x (era)

5

I0

15

n. 2'~ 0.5

20

(b)

Figure 2. Temperature distribution at different times for experiment 4 (d = 12.2 mrn, ~b = 0.431, S = 0.0, Ste = 0.282, Ra* = 5962.0): (a) bottom rake, 71 = 0.333 and (b) top rake, ~? = 0.667.

progresses and more PCM melts (and the region of convection increases), the temperatures show a smoother decline resembling that observed in a liquid saturated porous media in the presence of natural convection without phase change [28]. The melting proceeds faster along the top of the heated wall, and all the PCM has melted along the top rake by the end of 1.5 h, whereas along the bottom rake, it takes 2 h to melt all the PCM. At t = 2 h, the temperature distribution along the top rake is characterized by small steplike nonuniformities. Along the bottom rake, such variations have larger amplitudes. It is speculated that these nonuniformities are caused by instabilities in the transient, evolving flow structure and developing secondary flow. The results for experiment 6 are for an intermediate Rayleigh number (Ra* = 621.9) between those for experiments 4 and 7, and therefore are not included here for the sake of brevity. The results show that up to t = 0.5 h conduction is the predominant mode of heat transfer and, after convection sets in the rate of melting along the top, is higher than that along the bottom. The approximate melt thicknesses at t = 1 h along the top and bottom rakes are 9.2 and 8.2 cm, respectively. At the same time for experiment 4 (Fig. 2), the respective melt thicknesses are about 13.1 cm along the top and 10 cm along the bottom rakes. With decrease in ball size, the permeability and, hence, the Darcy number decreases, and the onset of convection is delayed and the flow is also weakened. The high thermal conductivity of the porous metal matrix, also mitigates the effect of convection. At t = 2.5 h, all the ice had melted along the top rake, and about 4 cm thick ice was left along the bottom rake, The temperature distributions at different times for experiment 7 are shown in Fig. 3. The average initial subcooling was - 0 . 2 ° C . At different times, the temperature in the melt region decreases monotonically with distance from the hot wall, suggesting very weak convection currents. The small size of the balls reduces the permeability and delays the onset of natural convection. The effect of convection is seen at t = 1.5 h, when the respective melt thicknesses along the top and bottom are 10 and 11.2 cm. A comparison of the melt thicknesses for experiments 4, 6 and 7 is given in Fig. 4. These melt thicknesses were determined using the fact that the solid/liquid interface is said to have crossed a particular thermocouple when the temperature recorded by the thermocouple is equal to 0°C. The temperatures in the solid region, about 1.5 cm away from the front, were in the range of 0.05 and 0. I°C below the fusion temperature. Therefore, the thermocouples in this region recorded temperatures very close to each other. This presented inaccuracies in ascertaining the exact location of the front from the temperature measurements. Also, the temperatures in the frozen region were measured at distances 8 mm apart. There-

226 S. Chellaiah and R. Viskanta 30

~

I

--

,

25

..... *.. Oh)

~

0 0-~~3

I

T

~

.

2o

o

~.~

201 t (h)

.

~3 Is

I

t

I

"-~

I

~ 2

O ; 0,25~

4

1.5

~

~.o

4

7

I0

z

5

II

4

2

~

/

!

3

Oi 5

15

I0

5

0

(a)

I0

15

20

t//,'

iment 7 (d = 3.2 mm, ~b = 0.395, S = 0.0, Ste = 0.32, Ra* = 383.4): (a) bottom rake, 7/ -- 0.333 and (b) top rake, ,7 = 0.667. fore, the location of the melting front at any instant of time was found by interpolation of the temperatures over the 8 mm region between the two adjacent thermocouples that the front crossed at times on either side of the time instant under consideration. It is estimated that the uncertainty in determining the front position is /:4 mm. It is evident from Fig. 4 that the strength of convection (as revealed by the difference in melting rates along the two rakes) decreases with decrease in the aluminum ball size, that is, experiments 4, 6, and 7, for reasons mentioned earlier. In the first 0.4 h, after the start of the melting, for experiment 7, the melting rates are the highest because of the higher Stefan number. After natural convection currents are established, the melting rates for experiment 7 are the lowest even though the Stefan number was higher than that for experiments 4 and 6. This clearly indicates that the additional temperature potential (provided by this increase in Stefan number) is not sufficiently large enough to induce flow and augment the melting rate

20

I

1

r

Exp. 4 ~

6 l

'

/

i5

\l z~

/ ¢/,.'." lo

---

?=0.333 --0.667

0

I 2

0 t

/

---

,7--o.333

4

6

(b)

x (cr.)

Figure 3. Temperature distribution at different times for exper-

3

/

I :3

4

(h)

Figure 4, Effect of bead size on the rate of melting for experiments 4 (d = 12.2 mm), 6 (d = 4.8 ram), and 7 (d = 3.2 mm).

0

2

8

t (h)

Figure 5. Effect of Rayleigh number on the rate of melting for experiments 1 (Ra* = 328.6), 2 (Ra* = 945.7), 3 (Ra* = 4629.8) and 4 (Ra* = 5962.0). to be higher than that for experiments 4 or 6. It should be noted that the effective thermal conductivity is almost the same (except for small variations caused by porosity changes) for all the experiments, whereas the role played by convection is system dependent. The lower permeability delays the onset of convection and mitigates the contribution of convection in accelerating the melting process.

Effect of Rayleigh Number With the other parameters (ie, aspect ratio, Stefan number and Darcy number) kept constant, the effect of Rayleigh number on the melting process can be examined by varying the hot wall temperature (Th). It should be noted that changes in Th also produce changes in the Stefan number. The effect of Rayleigh number can be studied by comparing the timeinterface position traces (Fig. 5) of experiments 1, 2, 3, and 4 with respective Rayleigh numbers of 328.6, 945.7, 4629.8, and 5962.0. For experiment 1 (smallest Ra*), the melting is controlled by conduction up to about 4 h (as revealed by the equal melting rates at the two heights). Then convective flow sets in, and the melting is accelerated along the bottom. This is in contrast to that observed in other experiments and is due to the density inversion of water (as explained in the following section). A comparison of the melting rates of experiments 3 and 4 reveals that a small increase of about 3°C in the hot wall temperature considerably accelerates the melting process. The times taken to melt a 17 cm layer (the location of farthest thermocouple in both the rakes from the hot wall) of ice along the top rakes for experiments 3 and 4 are 120 and 87 min, respectively. The corresponding times for bottom rake are 155 and 120 min. Thus, a 3°C increase in the hot wall temperature significantly accelerates the rate of melting. This is due to the increase in both the conductive heat transfer as well as the convective heat transfer. The temperature distributions for experiment 3 (not shown) and experiment 4 (Fig. 2) are qualitatively the same. But, near the completion of the experiment (t = 2 h for experiment 4 and t -- 2.5 h for experiment 3), the temperature distributions

Melting of Ice-Aluminum Balls System 22"/ 25F

7

T

20,

r

T

t (h)

!

2oF|

~3 ~.~0.5 4 5 ~

11

15~

1.0

t. 5 2.s

T

1

t~

1

(I

3

t~'t

e

~..s

2

k,.

I

I

Exp. 5

n~lq

//

|

4 x.o / 6

I

1 =.

,,,'/;;;;" 5

/

/

/// //f

---

,F~' 0

5

I0

(a)

15

5 x (cm)

I0

15

,7

•

o.333

~ 0.667 ;"

_

"

20

(b)

Figure 6. Temperature distribution at different times for experiment 5 (d = 12.2 ram, (~ = 0.431, S = 0.04, Ste = 0.257, Ra* = 5010.9): (a) bottom rake, ~/ = 0.333 and Co) top rake, ~/-- 0.667. show steplike variations. It is speculated that these nonuniformities are caused by secondary flow ceils present in the evolving flow field. From Fig. 5, it is observed that an increase in the Rayleigh number by more than an order of magnitude (for experiment 4, Ra* is about 20 times that for experiment 1), significantly increases the melting rate. At t = 87 min, the respective melt thicknesses along the top rake for experiments 1 and 4 are 3.2 and 17 cm, while the corresponding values for bottom rake are 3.2 and 13 cm. Effect o f Subcooling P a r a m e t e r When the PCM is subcooled, the sensible heat of the system must be added to bring it to the fusion temperature before it is melted. The effect of subcooling on the melting process can be studied by comparing the results of experiments 3 and 5. The temperature distributions at different times along the two thermocouple rakes for experiment 5 are shown in Fig. 6. Again, they show that melting along the top of the test cell is faster than that along the bottom. The heat transfer is purely by conduction up to t = 0.25 h, as shown by two linear temperature profiles, one in the liquid and other in the solid region. At t = 0.5 h, the effect of convection is seen to increase melting at the top. The temperature distributions for experiment 3 with an average initial subcooling of - 0 . 2 ° C (not shown) have revealed that the effect of natural convection can be noticeable even at t = 0.25 h [18]. At t = 0.5 h, the natural convection has accelerated the melting along the top rake, so that the respective melt thicknesses along the top and bottom are 6.1 and 4.6 cm. The melt thickness along the top rake at t = 2 h for experiments 3 and 5 are 17 cm and 15.5 cm, respectively. The corresponding values along the bottom rake at the same time are 14 and 11.6 cm. Thus, subcooling is seen to retard the melting process. A comparison of the melting rates for the two experiments along the two rakes is shown in Fig. 7. The interface locations at different times can be easily inferred from this figure. Effect o f Density Inversion o f Water Water attains its maximum density at 4°C, and natural convection currents will be accordingly influenced depending on

I 1

I 2 t (h)

I 3

4

Figure 7, Effect of subcooling parameter on the rate of melting for experiments 3 (S = 0) and 5 (S -- 0.04). the hot wall temperature in comparison to 4°C isotherm. The effect of this density inversion on the melting process is explained by comparison of experiments 1 and 4. The temperature distributions at different times for experiment 1 are shown in Fig. 8. The subcooling of PCM in this experiment is negligible, and the hot wall temperature is 4.42°C, which is slightly above the inversion point. Figure 8 shows that melting occurs at higher rates along the bottom rake than along the top rake. Water near the hot wall is heavier and warmer than that near the interface. Therefore, it descends along the hot wall and rises along the interface, forming a counterclockwise rotating cell. Thus, the hottest liquid is at the bottom, and the coldest is at the top of the test cell. This causes the melting rates along the bottom rake to be higher than that along the top rake. In the melting of ordinary PCMs, the flow is in the clockwise direction. The hot liquid rises along the hot wall, and the cooler melt descends along the interface, causing higher melting rates at the top than at the bottom [27]. In experiment 4 (Fig. 2), the hot wall is at 22.25°C. The clockwise flow driven by a temperature difference of 4°C in Io

{

I

i

I

tG'o 1 o 8 4 ~ ~ ?

6

I

I t d-o

J

'°

3 4 5 6

2 3 6 5

1 2 l 4

7

6

7

4 ~6

~5 5

4

2

0 1

0

I i0

5

(a)

I 15

0

l

I

I

5

i0

15

20

fb) Figure 8. Temperature distribution at different times for experiment 1 (d = 12.2 mm, ~ -- 0.431, S = 0.0, Ste = 0.061, Ra* = 328.6): (a) bottom rake, ~/ -- 0.333 and (b) top rake, = 0.667. x

(era)

228

S. Chellaiah and R. Viskanta I

I

i

]

m

I

i t O~)

t (h)

~g

] [

~ ~.~

o.~ 8

g / 4

/

j4

7

3 5

L

L

I

5

i0

15

(a)

0 ×

(cm)

I

5

I

1

ia

i5

]

20

(b)

Figure 9. Temperature distribution at different times for experiment 2 (d = 12.2 mm, ~b = 0.431, S = 0.0, Ste = 0.107, Ra* = 945.7): (a) bottom rake, ~/ = 0.333 and (b) top rake, = 0.667.

(a) i

the region bounded by the interface and the 4°C isotherm is practically overpowered by the counterclockwise flow driven by 18°C temperature difference in the region bounded by the hot wall and the 4°C isotherm. Therefore, melting occurs at a higher rate along the top. In experiment 2 (Fig. 9), the hot wall is at 8.4°C. The melting rates are nearly the same along the two rakes up to about 2.5 h. After this time, the convective flow sets in, and melting is faster along the bottom. The reason for this can be deduced from the temperature distributions. After the convective flow has been established, say at t = 6 h, the temperature distribution along the bottom rake reveals that a drop from 8.4°C to 4°C occurs in a small layer (5.6 cm wide), whereas the drop from 4°C to 0°C, occurs in a 9 cm wide layer. Thus, the counterclockwise cell in the region bounded by 0 and 4°C isotherms is larger and overpowers the clockwise flow between the hot wall and 4°C isotherms. Therefore, melting progresses at a higher rate along the bottom. In summary, when the hot wall is at 8°C or lower, for 12.2 mm aluminum balls, a unicellular counterclockwise flow is established, and more melting occurs at the bottom than at the top of the test cell. An increase in Th to higher values causes the melting to be faster along the top as in the case of pure PCMs.

Numerical Simulation of Melting Melting of the ice-aluminum balls system in experiment 5 was simulated using the enthalpy based model. As the initial temperature of the ice-aluminum system was - 6 . 4 ° C , the conduction heat transfer in the frozen region caused by this subcooling was also accounted for in addition to the temperature driven natural convection currents in the liquid. The predicted flow and temperature fields at r = 0.0169 (t = 1.5 h) are shown in Fig. 10. Results at other times are available elsewhere [18]. The influence of convection is clearly evident, as revealed by the higher melting rates along the top. In this simulation, the hot wall is at a temperature of 20.3°C. The flow in the region bounded by the interface and the 4°C isotherm is weak since only a temperature difference of 4°C drives the flow. Therefore, it is overpowered by the clockwise flow driven by a temperature difference of 16°C;

i,

,

.# /

,

, ,

r

,

/

,.//

t£

(b)

Figure 10. Flow and temperature field for experiment 5 at t : 1.5 h (d = 12.2 mm, ~ = 0.431, S = 0.04, Ste = 0.257, Ra* = 5010.9): (a) streamlines (¢/ kg/m s) and (b) isotherms

(K). between the hot wall and the 4°C, isotherm and only single clockwise rotating cell is predicted. The liquid rises along the hot wall, turns right at the top of the test cell, and the hot liquid impinges on the interface. The liquid is cooled as it flows downward. Therefore, the melting rate is faster at the top of the test cell than at the bottom. The streamline and temperature distributions are similar to those for pure natural convection in a porous media and qualitatively compare well with those of Chan et al [29]. Quantitative comparison is not possible since the experimental conditions are different. The flow patterns are also similar, although the flow is weaker because of the presence of the porous matrix, to those observed in the case of melting of pure substances [27]. At the start of the melting, there is a very thin layer of liquid near the hot wall [18], and the heat transfer is primarily by conduction. The liquid region is small, and its aspect ratio is very large. With the progression of melting, the temperature difference between the hot wall and the interface induces buoyancy-driven convection, causing higher melting rates at the top of the test cell than at the bottom (Fig. 10). The interface becomes curved, and the mean aspect ratio continu-

Melting of Ice-Aluminum Balls System 229 ,o , , , ~ , 0.8

AIIIII, ( I l,',~,~ L Tlll4IIll - _i_ .- -.-.- . . n ~0 ,:~

'

\',,',,'~.~ ~ ,a'.--~""-~

o.6 '2,,"~

0

., ,~

° ko \

?. " . t \

-~.

o

oz

'

&

'.. 0.2

'

~,-

04

~,

06

os

o

(a)

oz

o4

¢

o6

o8

,o

(b)

Figure 11. Comparison of experimental and predicted tempera-

ture distribution at different times for experiment 5 (d = 12.2 ram, ~ = 0.431, S = 0.04, Ste = 0.257, Ra* = 5010.9): (a) top rake, ~1 = 0.667 and (b) bottom rake, ~ = 0.333. ally decreases and influences the flow structure. In the solid region, the heat transfer is by conduction, driven by a temperature difference of 6.4°C between the cold wall and the melting front. The average Nusselt number at the hot wall as a function of time has been predicted [18]. When conduction is the only mode of heat transfer, the heat flux at the hot wall is inversely proportional to the square root of time [27]. Therefore, at early times in the melting process, the heat transfer at the wall is large, and the Nusselt number is high. But as convection sets in, the heat transfer rate decreases, and the Nusselt number reaches a quasi-steady value. This pattern agrees with the published f'mdings of melting from a vertical wall in the absence of porous matrix [27]. C o m p a r i s o n o f Predicted Results W i t h Experimental Data

A comparison between the predicted and the measured temperature distributions for experiment 5 at different times and along the two different vertical locations are presented in Fig. 11. In general, the agreement between the two results is good, but the temperatures in the liquid region are always overpredieted. The maximum difference between the measured and predicted temperatures is about 12% of the total temperature difference across the test cell and is within the experimental uncertainty. Any heat gains from the ambient to the test cell would tend to increase and not lower the temperatures;

0.8

,

,

o6r

/

J

/

o 2 [/'-. 0

,

°

"o . ..-"

, .... ,~,

r

i

i

i

0

0...5

LO

~5

r

20

therefore, the gains can be dismissed as the reason for the discrepancy. Certainly, the porosity and permeability are not uniform within distances of a few ball diameters from the melting front and the hot wall. The effect of wall channeling caused by the porosity variation near the wall is believed to be an important reason for increasing the temperature gradients near the hot boundary. The predicted temperature gradients based on uniform porosity and permeability are smaller than the measured ones, particularly at later times in the melting process. Figure 12 shows a comparison of the predicted and measured dimensionless interface positions. Reasonably good agreement is noted between the two results. The model prodiets higher rates of melting along the top rake than the bottom rake as observed in the experiment, but the model overpredicts the front locations at later times. This is because the temperature in the liquid region is predicted to be higher than that determined experimentally. There are a number of factors (some of which have been discussed earlier) that may contribute to the discrepancy between data and predictions; they include the following: (1) contraction of the liquid, (2) inappropriate permeability and effective thermal conductivity models for the porous media, (3) imprecise location of the front in the region where both the liquid and solid phases exist simultaneously, and (4) numerical inaccuracies caused by insufficiently fine grid. The major difficulty in the case of aluminum-water porous media (where the thermal conductivities of PCM and porous matrix differ by more than two orders of magnitude) is in predicting the effective thermal conductivity of the system. In this simulation, the effective thermal conductivity measured by Lindemann [23] was used. The value was based on a single experiment, and the accuracy of the experimental data has not been determined by Lindemann. Also, the porosity varies exponentially from a value of unity at the hot wall to about 0.431 in the region about 3--4 ball diameters away [30]. As the permeability is dependent on the porosity, this translates to higher permeability and lower resistance to fluid flow in this region. Additionally, the viscosity of water decreases by a factor of 1.8 from 1.9 x 10 -3 N • s/m 2 at 0°C to 0.996 x 10 -3 N . s/m 2 at 20.3°C near the hot wall [31]. These factors tend to produce steeper temperature gradients near the hot wall. The model has ignored these effects and, therefore, overpredicts temperatures in the liquid region. The calculated temperature gradients based on uniform porosity and permeability as well as constant viscosity and temperature are smaller than the measured ones, particularly at later times in the melting process. Finally, the model has not accounted for the contraction in volume of ice on melting. The calculation of the permeability in the liquid region (particularly where both the liquid and solid are present) may not be accurate [17]. For 12 mm balls, the permeability is large, and the flow is considerably influenced by its value. For example, the control volumes are relatively large (8.5 ram) to be represented by a node and assigned a single temperature. All of these factors can affect the flow structure and temperature distributions in the solid and liquid regions, and alter the solidification front shape and motion. PRACTICAL SIGNIFICANCE/USEFULNESS

" I00

Figure 12. Comparison of experimental and predicted interface

locations for experiment 5 (d = 12.2 nun, ~ = 0.431, S = 0.04, Ste = 0.257, Ra* = 5010.9).

This fundamental research was conducted to clarify the role of natural convection on the local rate of melting of a frozen porous medium consisting of ice as the phase change medium

230 S. Chellaiah and R. Viskanta and high thermal conductivity solid (aluminum balls) as the porous matrix. Physical situations of this type arise in geophysical, excavation, and latent heat-of-fusion energy storage applications. The aluminum balls-ice system was chosen for the study because the effective thermal conductivity of the phase change medium is large, and heat conduction in the frozen region may be significant and affect the melting rate and the front shape. The experiments reported here were made in an attempt to gather fundamental data that would provide basic understanding of the phenomena and could be used as basis for validating theoretical models. The experimental results are the first of the kind in which the effects of beat conduction in both the frozen and melted regions, natural convection of the melt, and density inversion in the liquid (water) are considered simultaneously. The results show that in spite of the high effective thermal conductivity of the porous medium, natural convection impacts significantly the melting rate and the shape of the melted region. The Rayleigh number for the porous medium is the most important parameter that controls the intensity of natural convection flow in the melted region. The results demonstrate that for sufficiently large Rayleigh numbers, natural convection is as important as conduction in the melted region and must be accounted for in an analysis of the melting rate and phase change front shape. For relatively small temperature differences between the heated wall and the fusion temperature of water (ie, Th - T f < 10°C), the density inversion of water can influence in a major way natural convection circulation in the melt and the shape of the melted region. CONCLUSIONS An experimental and numerical study of the ice-aluminum balls system has been performed. A number of different experiments using three different size balls and subcoolings have been conducted. The experiments have provided conclusive evidence that natural convection in the liquid region causes the melting front to become nonplanar and increases the rate of melting. The intensity of natural convection in the melt region depends on the Rayleigh number of the porous medium

Ra*. An enthalpy-based numerical model that considers both diffusion in the solid and liquid regions and buoyancy-driven convection in the liquid has been used to simulate the melting of liquid saturated porous media. The numerical predictions were compared with measured temperatures and interface positions, and good correspondence has been found. The possible reasons for the discrepancy between model predictions and the data have been discussed. The computational resources needed to obtain solutions were excessive. There is a need for flow visualization and nonintrusive diagnostics for temperature and melting front position measurement in porous media, both in the absence and in the presence of phase change. There is also a need for developing efficient, accurate, cost-effective numerical algorithms for solving the model equations for two- and three-dimensional solid/liquid phase change of porous media in the presence of natural convection. The work reported in this article was supported in part by the Heat Transfer Program of the National Science Science Foundation under Gram No. CBT-8313573.

NOMENCLATURE A C c Da d g" Ahf

H K k L Ra* S s Ste T t 0 u v v x y

aspect ratio, H/L, dimensionless Forchheimer's constant, C -- 0.55, dimensionless specific heat J/(kg -1 K - l ) Darcy number, K/L 2, dimensionless ball diameter (mm) gravitational acceleration, (m/s 2) latent heat of fusion, J kg - t height of liquid level (m) permeability (m 2) thermal conductivity (W m - i K - l) length of cavity (m) Rayleigh number for a porous medium, g~oKL(Th T f )q /i, l OLl, dimensionless supercooling parameter, cs(T f - T c ) / z ~ x h f , dimensionless interface position from cold wall (m) Stefan number, c j (Th - T f ) / A h f , dimensionless temperature (K) time(s) velocity vector, ~u + f o , m/s velocity component in the x-direction (m/s) volume (m 3) velocity component in the y-direction (m/s) horizontal distance from cold wall see Fig. 1, (m) vertical coordinate see Fig. 1, (m)

Greek Symbols c~ thermal diffusivity, k / p c (m 2 s - l ) /3 coefficient of thermal expansion ( K - I ) 7 volume fraction of liquid phase change material (PCM) in V f , l i t / V f , dimensionless volume fraction of liquid PCM in the volume element, V t / V = 07, dimensionless ~" vertical coordinate, s/L, dimensionless ~/ vertical coordinate, y / L , dimensionless 0 temperature, (T - T c ) / ( T h - T e L dimensionless dynamic viscosity ( N . s/m 2) v kinematic viscosity (m 2 s - l ) horizontal coordinate, x / L , dimensionless r time, tort/L 2, dimensionless ~b porosity or volume fraction of PCM in the volume element, V y / V , dimensionless ~b stream function, kg/(m s)

Subscripts eft f c h 1 m p s

effective fluid cold hot liquid maximum porous matrix solid REFERENCES

1. Lunardini, V. J., Heat Transfer in Cold Climates, Van Nostrand Reinhold Co., New York, 1981.

Melting of Ice-Aluminum Balls System 2. Sanger, F. J., Ground Freezing and Construction, ASCE Mech. Foundations Div., 94, 131-158, 1968. 3. Misnyck, Yu. M., Shonin, O. B., Ryabets, N. I., and Petrov, V. S., Application of Microwave Energy for Accelerating Excavation in Frozen Soil, in Final Proceedings of Fourth International Conference on Permafrost, pp. 268-272, National Academy Press, Washington, D.C., 1983. 4. ME Staff, Seasonal Thermal Energy Storage, Mech. Eng., 105, 28-34, 1983. 5. Metz, P. D., A Simple Computer Program to Model ThreeDimensional Underground Heat Flow with Realistic Boundary Condition, J. Solar Energy Eng., 105, 42-49, 1983. 6. Svec, O., Goodrich, L. E., and Planar, J. H. L., Heat Transfer Characteristics of Inground Heat Exchangers, J. Energy Res., 7, 263-278, 1983. 7. Fisher, K. M., The Effects of Fluid Flow on the Solidification of Castings and Ingots, PhysicoChem. Hydrodyn., 2, 311-326, 1981. 8. Cheng, K. C., Lunardini, V. J., and Seki, N., Proceedings of the 1987 International Symposium on Cold Regions Heat Transfer, ASME, New York, 1987. 9. Aung, W., and Yener, Y., Research Directions in Natural Convection, in Natural Convection Fundamentals and Applications, S. Kakac, W. Aung and R. Viskanta, Eds., pp. 1155-1171, Hemisphere, Washington, D.C., 1985. 10. Combarnous, M. A., and Bories, S. A., Hydrothermal Convection in Saturated Porous Media, in Advances in Hydroscience, Vol. 10, pp. 231-307, Academic Press, New York, 1976. 11. Cheng, P., Heat Transfer in Geothermal Systems, in Advances in Heat Transfer, T. F. Irvine, Jr. and J. P. Hartnett, Eds., Vol. 14, pp. 1-106, Academic Press, New York, 1978. 12. Goldstein, M. E., and Reed, R. L., Effect of Fluid Flow on Freezing and Thawing of Saturated Porous Media, Proc. Royal Soc. London, Series A, 365, 45-73, 1978. 13. Jumikis, A. R., Electrical Thawing of Frozen Soils, in Final Proceedings of Fourth International Conference on Permafrost, pp. 333-338, National Academy of Sciences Press, Washington, D.C., 1983. 14. Hromadka II, T. V., Guymon, G. L., and Berg, R. L., Comparison of Two- Dimensional Domain and Boundary Integral Geothermal Models with Embankment Freeze-Thaw Field Data, in Final Proceedings of Fourth International Conference on Permafrost, pp. 509-513, National Academy of Sciences Press, Washington, D.C., 1983. 15. Weaver, J. A., and Viskanta, R., Melting of Liquid-Saturated Porous Media, Int. J. Heat Mass Transfer, 29, 1943-1951, 1986. 16. Okada, M., and Fukamoto, T., Melting Around a Horizontal Pipe Embedded in a Porous Medium, Trans. Jap. Soc. Mech. Engrs., Series B, 48, 2041-2049, 1982. 17. Beckermann, C., and Viskanta, R., Natural Convection Solid/Liquid

231

Phase Change in Porous Media, Int. J. Heat Mass Transfer, 31, 35-46, 1988. 18. Chellaiah, S., An Experimental and Numerical Investigation of Solid/Liquid Phase Change of Water Saturated Porous Media in the Presence of Natural Convection, Ph.D. Thesis, School of Mech. Eng., Purdue University, West Lafayette, Ind., 1988. 19. Gebhart, B., and Mollendorf, J. C., Buoyancy-Induced Flows in Water Under Conditions in which Density-Extrema May Arise, J. Fluid Mech., 89, 673-707, 1978. 20. Nguyen, T. H., Vasseur, P., and Robillard, L., Natural Convection between Horizontal Concentric Cylinders with Density Inversion of Water for Low Rayleigh Numbers, Int. J. Heat Mass Transfer, 25, 1559-1568, 1982. 21. G-ebhart, B., Jaluria, Y., Mohafan, R. L., and Sammakia, B., Buoyancy-Induced Flows and Transport, Hemisphere, Washington, D.C., 1978. 22. Veinberg, A. K., Permeability Electrical Conductivity, Dielectric Constant and Thermal Conductivity of a Medium with Spherical and Ellipsoidal Inclusions, Sovient Physics Doklady, 11,593-595, 1967. 23. Lindemann, W., Experimental and Theoretical Studies on Design and Calculations for Latent Heat Storage, M.S.M.E. Thesis, School of Mech. Eng., Purdue University, West Lafayette, Ind., 1986. 24. Ward, J. C., Turbulent Flow in Porous Media, J. Hydraul. Div. ASCE, 90, 1-12, 1964. 25. Patankar, S. V., Numerical Heat and Fluid Flow, Hemisphere, Washington, D.C., 1980. 26. Benenati, R. F., and Brosilow, C. B., Void Fraction Distribution in Beds of Spheres, AIChE Journal, 8, 359-361, 1962. 27. Viskanta, R., Natural Convection in Melting and Solidification, in Natural Convection Fundamentals and Applications, S, Kakac, W. Aung, and R. Viskanta, Eds., pp. 845-877, Washington, D.C., 1985. 28. Inaba, H., and Seki, N., An Experimental Study of Transient Heat Transfer Characteristics in a Porous Layer Enclosed Between Two Opposing Surfaces with Different Temperatures, Int. J. Heat Mass Transfer, 24, 1854-1857, 1981. 29. Chart, B. K. C., Ivey, C. M., and Barry, J. M., Natural Convection in Enclosed Porous Media with Rectangular Boundaries, J. Heat Transfer, 92, 21-27, 1970. 30. Chandrasekara, B. C., and Vortmeyer, D., Flow Model for Velocity Distribution in Fixed Porous Beds Under Isothermal Conditions, Warme-und Stoffubertragung, 12, 105-111, 1979. 31. Weast, R. C., Editor, Handbook of Chemistry and Physics, 60th ed., pp. El0, F5-6, F449-451, D210-211, CRC Press Inc., Cleveland, Ohio, 1984.

Received October 31, 1988; revised May 1, 1989