# Experimental and numerical study of heat transfer in horizontal concentric annulus containing phase change material.

INTRODUCTION

In many engineering applications such as thermal energy storage using phase change materials (PCMs), the melting and solidification are important phenomena. Melting of the phase-change material gives rise to natural convection, while the flow structure could significantly affect phase change process. The essential feature of melting (or solidification) of PCM is the existence of the moving interface between two phases with time. The convection influences the morphology of this solid-liquid interface. This problem was studied as early as in 1831 by Lame and Clapeyron (Huntler and Kutter, 1989). The sequence of papers written by Stefan has given his name to this type of problem as "Stefan problem." The problem of heat transfer with phase change can be formulated, considering either the temperature or the enthalpy as the dependent variable. When temperature is considered as dependent variable, the energy equations for both the phases are to be written independently and then coupled for the interface. In fact this technique requires the knowledge of the interface position explicitly for the determination of temperature, complicating the formulation and method of solution. The problem has been solved by a method of immobilization of the interface (Duda et al., 1975). When the enthalpy is considered as the dependent variable the knowledge of the interface position is not required and a single enthalpy equation for the whole domain suffices. Huntler and Kutter (1989) demonstrated the enthalpy method for heat transfer problems with moving boundaries associated with phase change. Pannu et al. (1980) calculated temperature and velocity profiles, rates of heat transfer, movements of melting interface in vertical and horizontal cylinders by finite difference method. The integral method has been applied to study the effects of natural convection on the melting of solid around the horizontal cylinder (Yao and Chen, 1980; Yao and Cherney, 1981). The perturbation and numerical solution of melting around the heated horizontal cylinder for an isothermal boundary condition were studied previously (Reiger et al., 1982; Prusa and Yao, 1984a,b) and the melting process in a rectangular enclosure driven by coupling of heat conduction in the solid phase and natural convection in the melt of the PCM has been analyzed (Bernard et al., 1986). Viswanath and Jaluria (1993, 1995) carried out a numerical modelling of the macroscopic melting processes in rectangular enclosures using an enthalpy--porosity formulation and has been used to circumvent the need to track the moving interface at every time instant. In this case a single energy equation is valid over all regions including the interface. The SIMPLE algorithm was employed to implement the solution. A numerical study of melting of PCM around a horizontal circular cylinder of constant wall temperature and in presence of natural convection in melt region was presented by Ismail and da Silva (2003a,b). In a recent literature (Lamberg et al., 2003), the comparison of experimental data with numerical results was obtained where only energy equation was solved by modelling convective source term with an effective heat capacity method. Ettouney et al. (2005) studied heat transfer characteristics of PCM in a vertical annulus and developed empirical correlations involving Nusselt number, Rayleigh number, Steffan number, and Fourier number.

Ng et al. (1998) in their study employed the finite element method to simulate the convection dominated melting of a PCM in a cylindrical horizontal annulus heated isothermally from an inside wall. The effects of Rayleigh number on the melting rate as well as the evolution of the flow pattern were examined. The authors observed that the increasing Rayleigh number promotes the heat transfer rate and the multiple cellular pattern was found to occur at high Rayleigh number (>[10.sup.6]). Khillairkar et al. (2000) studied melting in a concentric horizontal annulus of arbitrary geometric arrangement such as external square tube with circular tube inside (Type A) and circular external tube with a square tube inside (Type B). The authors studied the effects of the Rayleigh numbers as well as the heating of either the inside surface or the outside surface and in some cases both the surfaces at a temperature above the melting point of the PCM. To account for the physics of the time wise evolution of flow at the solid--liquid interface the well known enthalpy--porosity model was employed in the fixed grid method. For both the horizontal annuli of Type A and Type B, it is observed that the effects of heating both the surfaces is the same as the heating of inside surface or the outside surface separately until there is an interaction between the two melt-zone. It is observed that the melting rate is faster due to good mixing between melt zones. This suppresses the thermal stratification attained in the horizontal annuli of both Type A and B. The thermal stratification occurs in the upper part of the cavity due to the fact that the energy charged to the system is mainly carried upward by the free convection. Thus, the energy is used to raise the temperature of the melt instead of melting the PCM. So, the energy storage in the system is more in the form of sensible heat than in the form of latent heat. Balikowski and Mollendorf (2007) recently studied the performance of PCM in a horizontal annulus of a double pipe heat exchanger in a water circulating loop. PCM was placed in the annular region of a double pipe heat exchange with water circulated in the inside pipe. Experiments were performed in which the PCM would absorb or release heat at various temperatures and water flow rates. The effect of different water flow rates on the heat transfer rates were examined. The difference between temperature of hot water and melting point of the PCM was between 2 and 6[degrees]C. To avoid any thermal penalty and the flow rate of hot water were so maintained that there was virtually no difference between the inlet and outlet hot water temperatures.

In this paper, heat transfer is studied experimentally and numerically for PCM (paraffin wax) encapsulated in the horizontal annulus of two coaxial circular cylinders with isothermal and adiabatic condition at the inner and the outer cylinder, respectively. Two-dimensional unsteady Navier-Stokes' equations with energy equations were solved by finite volume method (FVM) using SIMPLE scheme to track the temperature at different points and moving solid--liquid interface numerically. Enthalpy equation was considered as the energy equation. To solve this phase change problem in a complex geometry of an annulus formed by concentric horizontal tubes the enthalpy--porosity method has been used. Comparison of time-temperature history at different points obtained from physical experimentation and numerical experimentation is made. Analysis of the resulting interaction between the flow and the heat transfer in the annular region is considered through the variation of [epsilon] (eccentricity) and the inclination angle of eccentricity ([psi]). The numerical results obtained are displayed through streamlines, isotherm contours and velocity distribution in computational domain. The variation of net circulation and analysis of the heat flux along the outer surface of the inner cylinder is done for different eccentric positions.

PHYSICAL PROCESS

In this experiment the hot water at constant temperature was passed through the inner tube. The heat is conducted through the metal wall to the surrounding wax. The innermost wall layer of wax in the annulus adjacent to the metal wall melts quickly as the temperature of water is greater than the melting point of wax, and a layer of melt surrounds the metal tube. In this situation subsequent heat transfer will take place through the layer of molten wax. The layer of molten wax is initially stagnant, and gets heated. Due to difference in density in its body an internal circulation starts for heat transfer to take place. Thus, the mechanism during melting is primarily natural convection and conduction in the liquid phase and only conduction in the solid phase. In developing the natural convection model for the horizontal annulus, the following assumptions could be made:

(1) The temperature distribution along the axis is negligible.

(2) The thermal resistance of the pipe wall is negligible.

(3) The heat loss from the outer cylinder is negligible.

The first assumption reduces the problem to a two-dimensional one, simplifying the mathematical model. The other two assumptions defined the boundary conditions.

EXPERIMENTAL CONFIGURATION AND PROCEDURE

PCM Used in the Experiment

The PCM used in the experiment is paraffin wax. Its thermophysical properties are given in Table 1.

Experimental Set-Up

The set-up (Figures 1 and 2) consists of a hollow steel cylinder with internal diameter 10.16 cm, length 30.48 cm and the internal pipe of diameter 1.905 cm is placed coaxially and horizontally.

The annular space so formed contained the PCM. The configuration is similar to a concentric pipe heat exchanger. On the body of the outer cylinder eight thermo-couple pockets, four at each side diametrically opposite to each other were constructed. Then the thermo-couples were placed in the annulus maintaining equal distance starting from the wall of the inner tube to the inner wall of the outer cylinder at equal distances along the length of the cylinder as shown in the Figure 1. The stem connected to the outer steal cylinder allows the molten wax to be filled in the annulus. When the mass of the wax in the annulus melts during the heating cycle the volume expansion of the wax is accommodated in the stem and avoids any spillage of the wax. This assembly containing PCM in the annulus was put in a Perspex box to reduce heat loss. The void space in the Perspex box was then filled with glass wool loosely allowing voids to be filled up partially by air increasing the overall insulation property of the space separating the ambient from the outer surface of the outer cylinder. Two T-joints were attached to the two ends of the inner tube for placing two more thermo-couples to measure the inlet and outlet temperatures of the flowing water.

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

Experimental Procedure

For heating, at first the water bath was set to a desired fixed temperature. When the water reached a steady state temperature hot water was passed through the inner tube at a constant rate. The temperatures of different thermo-couples placed radially in the annulus were recorded in real time in the computer by HP-bench-link software with the help of the data logger (data-acquisition unit-HP 34970 A) connected via line driver to get the heating curves for different radial points. The recording of temperature data was done at 5 s interval and thermo-couple used was T-type with accuracy 0.5-1[degrees]C. The estimated accuracy of the calibration of the flow meter and the thermo-couples is about [+ or -]5%. The resulting accuracy of the inference of the heat transfer is about [+ or -]10%. The estimated accuracy of time and length measurement is about [+ or -]0.1 s and 2-3 mm, respectively. After both the thermo-couple placed at the outermost position of the annulus reached the melting point of the paraffin wax, the experiment was terminated.

MATHEMATICAL MODEL

For the numerical simulation of melting process two methods are generally adopted. These are transformed grid (T-grid) method and enthalpy--porosity methods. In this paper, the latter method with fine meshes is employed for the numerical solution of melting of paraffin wax in the horizontal concentric annulus. The accuracy of interface tracking has been achieved by considering finer meshes for this complex geometry.

Governing Equations and Enthalpy Formulation

In the enthalpy--porosity formulation, enthalpy has been taken as the dependent variable in place of temperature, along with other primitive variables u, v, and P. In this scheme the liquid mass fraction (f), defined as the ratio of the liquid mass to the total mass in a given computational cell plays an important role. Then the enthalpy content of that computational cell will be:

H = fL + [C.sub.p]T (1)

where the enthalpy content of the solid has been taken as reference enthalpy. Here the heat capacity ([C.sub.p]) can vary with the change of phases.

Now, in other way round the liquid mass fraction can be obtained from the enthalpy in the following way:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

The continuity and momentum equations for the dependent variable u, v, P, considering Boussinesq' approximation, are the following:

[partial derivative]u/[partial derivative]x + [partial derivative/[partial derivative = 0 (3)

[rho] ([partial derivative]u/[partial derivative]t + u [partial derivative]u/[partial derivative]x + v [partial derivative]u/[partial derivative]y) = [partial derivative]p/[partial derivative]x + [mu] ([[partial derivative].sup.2]u/[partial derivative][x.sup.2] + [[partial derivative].sup.2]u/[partial derivative][y.sup.2]) + [S.sub.u] (4)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

and the energy equation is:

[rho] ([partial derivative](cT)/[partial derivative]t + u [partial derivative](cT)/[partial derivative]x + v [partial derivative](cT)/[partial derivative]y) = k ([[partial derivative].sup.2]T/[partial derivative][x.sup.2] + [[partial derivative].sup.2]T/[partial derivative][y.sup.2] - [rho]L [partial derivative]f/[partial derivative]t (6)

Source Terms

In the Equations (4) and (5) the source terms [S.sub.u] and [S.sub.v] have been chosen in a manner such that they can suppress the velocity fields in the solid phase. When the PCM is in solid phase it takes large values, suppressing the velocity field. Also, with the increase of 'f' it should go smaller and finally be 0. One of the most popular and successful methods has been given by Viswanath and Jaluria (1993, 1995). They formulated the source term from the Karman-Kozeny equation for flow in porous media, where the equations resemble a Darcy-type law (where the velocities are proportional to the pressure gradient). Following are the source terms:

[S.sub.u] = -C (1-[f.sup.2]/[f.sup.3]+b) and [S.sub.v] = -C (1-[f.sup.2]/[f.sup.3]+b) v

The constants C and b can be suitably chosen to simulate the flow characteristics in porous media. From the equations it is clearly evident that when the liquid fraction is 0 then the source term dominate the equation and force the velocities to become zero in the region. The source term becomes negligible when liquid fraction is 1.

Boundary Conditions and Solution Method

In conformity with the experimental conditions the following boundary conditions have been imposed.

(1) No-slip condition at the surfaces of the tubes (inner and outer), forming the annulus.

(2) Isothermal condition at the inner surface of the annulus and adiabatic condition at the outer surface of the annulus.

Solution of the above mentioned Equations ((3)-(6)), after discretising by FVM, have been done by the SIMPLE algorithm (Patankar, 1980), introducing power law scheme and the deferred correction method (Ferziger and Peric, 2002) in the interpolation scheme for the convection term. After proper solution of these equations at each time step, the enthalpy and liquid mass fraction have been updated by the Equations ((1)-(2)).

Code Validation

Our code has been validated with the benchmark problem of natural convection in a square enclosure with a hot and cold wall. Comparisons were made with the previous works on phase-change materials, specifically the vertical annulus case. In all those cases the results were in good agreement.

Time-Efficiency of the Code

The code is also efficient in the sense of time-efficiency. Calculation for our problem at a grid-size of 60 x 120, in a P4 workstation with 512 MB RAM have a ratio of experiment time and numerical simulation time of 20:1.

NUMERICAL EXPERIMENTATION

For a better insight in the natural convection in phase change heat transfer in an annulus, the effects of the eccentricity as well as the angle of inclination of eccentricity are considered to be of great importance. It is already observed that the effects of Rayleigh number and Prandtl number on the natural convective heat transfer is well established. In this study the eccentricity is varied from 0 to +0.5. Similar such positions of eccentricity is varied among the different angle of inclination ([psi]) between [pi]/2 to [pi]/2 (e.g. -90[degrees], -60[degrees], -30[degrees], 0[degrees], 30[degrees], 60[degrees], 90[degrees]) as shown in Figure 3.

[FIGURE 3 OMITTED]

Due to the eccentric position the annulus becomes enlarged as well as narrowed in the direction of angle of inclination giving rise to the different fluid dynamic phenomena. Eminent thermal effects are also observed in both the parts of the annulus across the vertical midplain along the centre of the inner cylinder. In the numerical experiments the dimensions of the annulus assumed to be of 1 m.(i.d.) and 5 m (o.d.), which are larger than that of experimental set-up. The assumed physical properties of the material filled in the annulus are given in Table 2.

The boundary conditions of the inner and outer surfaces of the annulus are assumed to be isothermal and adiabatic, respectively. The above numerical experimental conditions are not varied as the dependence of heat transfer and fluid dynamics on Rayleigh and Prandtl numbers are well established.

The angular positions on the surface of the inner cylinder are shown in Figure 4. Due to the formation of thermal plume by the natural convective flow in the melt phase, the dimensionless heat flux value changes from the bottom point to the top point on the surface of the inner cylinder. The dimensionless heat flux is found by, -d[THETA]/dR at the surface where [THETA] = [T.sub.H]-T/[T.sub.H]-[T.sub.M], R = r/[r.sub.o] -[r.sub.i], [T.sub.H] = temperature of the inner surface, T = temperature of any point in the annulus, [T.sub.M] = melting temperature of the PCM, r = radius of any point on the annulus, [r.sub.i] = inner radius, [r.sub.o] = outer radius of the annulus.

[FIGURE 4 OMITTED]

The numerical simulation is done for a number of grid sizes in both the radial and circular directions. The mesh size is chosen after the performance of several tests as 200 x 200. For the accuracy of the numerical results and the convergence the error level is taken to be [10.sup.-5]. All the simulations are carried up to a time 30 000 s ([F.sub.0] = 3.08).

RESULTS AND DISCUSSION

The real time temperatures at different points away from the outer surface of the inner tube, that is, at four cross-sectional points in upper and lower position of the annulus, were measured. The position of the censors are clearly marked in Figure 2, which is a cross-sectional view of the horizontal cylindrical annulus. The positions of the temperature censors attached to the outer surface of the inner tube at an equal distance from the two ends along the length of the cylinder are shown in Figure 1. As temperature difference of outlet and inlet water was very small (<1[degrees]), it is assumed that the heating takes place in isothermal condition. The temperature difference between these two point at similar radial position but at different longitudinal points are very small and conforms the assumption 1 of the Physical Process section. The experiment was performed for different temperature of water (75, 80, 85, and 90[degrees]C). From the experiment the time--temperature history for eight points was obtained. The data for 15 000 s or more were recorded till the temperature of the outermost node reached the melting point. One representative graph of time-temperature history for hot water at 80[degrees]C at the points U3 and L3 is given below in Figure 5.

The time--temperature history for the wax in the upper and lower half of the annulus shows a substantial difference. The "time--temperature history" of a point in the upper half of the annulus is quite different from the same point in the lower half having same radial position. This is because as soon as the hot fluid is allowed to flow through the inner tube, the tube wall reaches a higher temperature than the melting point of the wax in the annulus (in all our experimental cases) forming a molten layer of wax around the inner tube. The molten layer grows, receiving heat from the hot fluid flowing inside the inner tube. Due to natural convection, buoyancy-driven upward flow of the molten wax enhances the supply of heat to the wax in the upper part of the annulus and solid--liquid interface moves faster in the upper half of the annulus. The experimental result supports this fact by showing a higher growth rate of melting front in the upper half of the annulus. It is evident from the previous discussion and the experimental result that the heat transfer phenomena here consist of both the convective and conductive mode of energy transfer. To explain the physical phenomena, time--temperature history inside the annulus for 15 000 s was simulated, by our mathematical model. Other variables like velocities, pressure and interface positionwere also being tracked for this time-span. Velocity vector profiles at three different time points 5000, 10 000, and 15 000 s by simulated data are given in Figure 6, in which only one-half of the annulus is shown, as both the halves are symmetrical showing similar features. White parts show the solid parts of the paraffin wax and the rest is molten portion having different velocities in different part with the progress of the melting.

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

The comparisons of experimental and simulated time-temperature history are given in the Figures 7 and 8. The figures indicate a good agreement between experimental and simulated data in these cases. The discrepancies may be due to sudden dislodging of a piece or pieces of semi-molten wax in the molten wax. In the positions nearer to the outer wall (U2) the simulated temperature is higher than the experimentally predicted one. This may be due to heat loss from the outer cylindrical surface even with the insulation provided around. Similar accuracy in prediction by simulated data of this order was prevalent in the cases of other hot fluid temperatures (80, 85, and 90[degrees]C).

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

Effects of Eccentricity

The effect of eccentricity [epsilon] on the kinematic and thermal field is studied for Ra = [5.7710.sup.7], Pr = 46.5 for the inclination angles [phi] =, -90[degrees], -60[degrees], -30[degrees], 0[degrees], 30[degrees], 60[degrees], 90[degrees] and [epsilon] = 0.0, 0.25, 0.50.

For the concentric annulus the streamlines and the isotherms are kidney shaped or crescent shaped and symmetric along the vertical midline along the centre of the inner cylinder. On both sides of the inner cylinder, the fluid flows up along the heated inner cylinder wall and down along the melt front. At a time so far elapsed with the fluids in the central core vortices are almost stagnant, and both sides have equal maximum stream function values, that is, [[psi].sub.max].

For an eccentric annulus, the flow pattern is also generally kidney shaped but is asymmetric along the midline. The position and formation of the crescent on two sides change with respect to the eccentricity [epsilon]. As the inner cylinder moves from the centre to the right side the extent of the left kidney gradually increases as the gap increases. The centre of the closed loop moved upward in the annulus. But the right crescent gradually reduces in size until the two crescent shaped flow region appears in both the left and right sections of the annulus when the eccentricity reaches a critical value. It is also observed that the eccentricity generally increases the strength of the flow circulation, that is, the value of the maximum stream function on the "intrusion" side of the flow region and decreases the flow strength of the "expanded" side of the flow. The imbalances of the flow strength between the two sides lead to a net circulation around the eccentric annulus excepting in the vertical eccentric cases.

[FIGURE 9 OMITTED]

The primary attention will be focused upon a problem of heat transfer by convection in a system of coaxial cylinder at fixed temperatures leading to approximately double increase of a flow compared with that determined only by heat conductivity. In the concentric case the strength of buoyancy-driven fluid motion as well as the enhancement in the convective heat transfer with the increase in Rayleigh number are well established. As the Rayleigh number becomes sufficiently high the inner cylinder is surrounded by the thermal boundary layer which starts at the bottom of the cylinder and leaves at the top. The inner cylinder local heat transfer coefficient varies with the angular position as the boundary layer develops from its bottom and is shed as a plume at the top. This plume heats the outer surface of the melt, causing the heat transfer coefficient at that position to be much higher than that for the inner one at the top. As the flow moves down along the outer surface of the melt, there decreases very rapidly. It is also found that the thermal stratification in the lower zone of the annulus weakens the fluid motion in that zone. One single recirculating vortex is found in a half domain and the strength of the flow obviously increases with the Rayleigh number. The effects of the eccentricity and its inclination angle on the heat transfer is considered with a fixed Rayleigh number. The value of the dimensionless eccentricity is varied from 0 to +0.5 and the inclination angle of inclination is varied from -90[degrees] to 90[degrees] by an interval of 30[degrees], by fixing [epsilon]. For a fixed inclination angle at the vertical axis of symmetry that is [phi] = -90[degrees] or 90[degrees], when the inner cylinder moves downward it causes a stronger fluid motion in the upper zone and also increases the temperature gradient in the lower zone, resulting in higher heat transfer between two cylinders, Figure 9 in comparison to the concentric case. On the other hand when the inner cylinder moving upward provides least favourable circumstances for the development of natural convection (Figure 10). Here the concentric cell centre has moved towards [phi] = 90[degrees]. The circulation cell is less powerful than in the concentric case. Clearly the top narrow mid gap between the two cylinders inhibits the convective motion in the fluid. Thus, the isotherms shows maller temperature increase in comparison to the concentric case. The bottom of the annulus near [phi] = -90[degrees] exhibits an enlarged stagnant region of fluid. In each case the motion of the fluid is clockwise on the right-hand side of the annulus. The region of almost stagnant fluid exists near the lower zone and this zone increases when the inner cylinder moves upward along the vertical axis increasing resistance to flow. Since the eccentricity due to the position of the inner cylinder in the upper side of the horizontal axis provides the least favoured circumstances for the development of natural convection, both the size and strength of flow of momentum are markedly reduced.

[FIGURE 10 OMITTED]

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

The heat flux for the natural convective heat transfer from the surface of the inner cylinder to the molten PCM changes along the surface from the bottom to the top point ([theta] = 0[degrees] -180[degrees]). The heat flux has the maximum value at the bottom and decreases to a minimum value at the top point, that is almost at 180[degrees] (e.g. Figures 11-14). In the case of concentric annulus, as there is virtually no net circulation, the variation of heat flux with position on the surface of the inner cylinder is almost symmetric about the vertical midline. It is observed that with the change in the vertical eccentricity the heat flux value reaches the maximum for [epsilon] = 0.5 in the direction of [phi] = -90[degrees]. On the other hand at the same value of the eccentricity in the direction of [phi] = 90[degrees] the heat flux at the bottom point on the surface of the inner cylinder reaches a very low value. The changes in the angle of inclination of the eccentricity the heat flux variation along the surface of the inner cylinder are not symmetric about the vertical midline due to the presence of net circulation. For a particular angle of inclination when the inner cylinder is in the first and the fourth quadrant the maximum thermal flux value are different. The maximum thermal flux is higher in the case of first quadrant than that in fourth quadrant. It is also observed that there is a shift of angular position for the occurrence of the minimum heat flux for the cases of same inclination angle but at different angle (Figures 12-14).

[FIGURE 13 OMITTED]

[FIGURE 14 OMITTED]

Effects of Angle of Inclination of Eccentricity

The net circulation [[psi].sub.0], being the difference between the [[psi].sub.max] values of the "intrusion" and "expanded" zones of the annulus, changes with the inclination angle of the eccentricity. When the inclination angle is -[pi]/2 or [pi]/2, that is, the case of vertical eccentricity, the net circulation is absent or very small. The net circulation is also symmetric along the axis of [phi] = -30[degrees] and has a maximum value. This also indicates the imbalance of fluid convection between the two sides of the inner cylinder is the reason for the net circulation around the annulus. A new phenomenon identified by this study is that when the inclination angle of eccentricity is not vertical, a net circulation is identified around the inner cylinder. For a given [epsilon] the magnitude of net circulation is expected to approach maximum when the angle of inclination is horizontal. The strength of the net circulation reaches a maximum for a given cylindrical annuli with the eccentricity in the range of 0 < [epsilon] < + 0.5, for the angle of eccentricity [phi] = -30[degrees] (Figure 15). This is due to the fact that the natural circulation of the molten phase creates an unmelted or solid phase front as an eccentric zone around the molten phase in addition to its eccentric position with respect to outer cylinder. The maximum value of the dimensionless heat flux changes with the eccentricity as well as the angle of inclination of eccentricity (Figure 16). It is observed that the maximum value of dimensionless heat flux reaches the highest value at [epsilon] = 0.5 and [phi] = 90[degrees].

[FIGURE 15 OMITTED]

[FIGURE 16 OMITTED]

CONCLUSIONS

A computational model for the prediction of the thermal behaviour of paraffin wax during its melting within an annulus formed between two horizontal cylinders has been presented. The model rests on solving Navier-Stokes equation with the consideration of the isothermal condition in the inner core and adiabatic condition at the outer surface of the annulus. The velocity vector field and the isotherms obtained from the computation confirm the conductive and conductive-convective heat transfer phenomena existing in the solid and the molten phases, respectively. The performance of the present model is verified with the experimental findings. The predicted result shows good agreement with the experimental results.

The results of the mathematical experiment obtained by varying both the eccentricity as well as their angle of inclination identify some new phenomena of flow. When the inclination angle of eccentricity is not in the vertical direction a net circulation is identified around the surface of the inner cylinder. For a given eccentricity the magnitude of the net circulation approaches maximum value when its angle of inclination approaches -[degrees]30? (Figure 15).

The analysis of heat flux along the surface of the inner cylinder indicates that for a fixed angle of inclination heat flux reaches a maximum value for the highest eccentricity ([epsilon] = 0.5) in the fourth quadrant and reaches a minimum value when the eccentricity is the maximum in the first quadrant. Similar trends are followed for all other angle of inclinations (Figures 11-14).

The maximum heat flux value increases monotonically both with angle of inclination and eccentricity (Figure 16).

Manuscript received March 25, 2007; revised manuscript received October 17, 2007; accepted for publication January 13, 2008.

REFERENCES

Balikowski, J. R. and J. C. Mollendorf, "Performance of Phase Change Materials in a Horizontal Annulus of a Double--Pipe Heat Exchanger in a Water Circulating Loop," ASME J. Heat Transf. 129, 265-272 (2007).

Bernard, C., D. Gobin and A. Zanoli, "Moving Boundary Problem: Heat Conduction in the Solid Phase of a Phase-Change Material during Melting Driven by Natural Convection in the Liquid," Int. J. Heat & Mass Transf. 29, 1669-1681 (1986).

Duda, J. L., M. F. Malone, R. H. Notter and J. S. Ventras, "Analysis of Two-Dimensional Diffusion-Controlled Moving Boundary Problems," Int. J. Heat Mass Transf. 18(7), 901-910 (1975).

Ettouney, H., H. El-Dessouky and A. Al-Ali, "Heat Transfer During Phase Change of Paraffin Wax Stored in Spherical Shells," J. Solar Energy Eng. 127, 357-365 (2005).

Ferziger, J. H. and M. Peric, "Computational Methods for Fluid Dynamics," Springer-Veralg Inc, (2002).

Huntler, L. W. and J. R. Kutter, "The Enthalpy Method for Heat Conduction Problems with Moving Boundaries," ASME J. Heat Transf. 111, 239-242 (1989).

Ismail, K. A. R. and M. D. E. da Silva, "Melting of PCM Around a Horizontal Cylinder With Constant Surface Temperature," Int. J. Therm. Sci. 42, 1145-1152 (2003a).

Ismail, K. A. R. and M. D. E. da Silva, "Numerical Solution of the Phase Change Problem Around a Horizontal Cylinder in the Presence of Natural Convection in the Melt Region," Int. J. Heat Mass Transf. 46, 1791-1799 (2003b).

Khillairkar, D. B., Z. X. Gong and A. S. Majumder, "Melting of Phase Change Material in a Concentric Horizontal Annuli of Arbitrary Cross Section," Appl. Therm. Eng. 20, 893-912 (2000).

Lamberg, P., R. Lehtiniemi and A. M. Henell, "Numerical and Experimental Investigation of Melting and Freezing Processes in Phase Change Material Storage," Int. J. Therm. Sci. (2003).

Ng, K. W., Z. X. Gong and A. S. Majumder, "Heat Transfer Convection-Dominated Melting of a Phase Change Material in a Horizontal Annulus," Int. Commun. Heat Mass Transf. 25(5), 634-640 (1998).

Pannu, J., G. Joglekar and P. A. Rice, "Natural Convection Heat Transfer to Cylinders of Phase Change Material Used for Thermal Torage," AIChE Symp. Ser. 76, 47-55 (1980).

Patankar, S. V., "Numerical Heat Transfer and Fluid Flow," Hemisphere, Washington, D.C. (1980).

Prusa, J. and L. S. Yao, "Melting Around a Horizontal Heated Cylinder: Part-I Perturbation and Numerical Solutions for Constant Heat Flux Boundary Condition," ASME J. Heat Transf. 106(2), 376-384 (1984a).

Prusa, J. and L. S. Yao, "Melting Around a Horizontal Heated Cylinder: Part-II Numerical Solution for Isothermal Boundary Condition," ASME J. Heat Transf. 106(2), 469-472 (1984b).

Reiger, H., U. Projan and H. Beer, "Analysis of the Heat Transport Mechanisms during Melting Around a Horizontal Circular Cylinder," Int. J. Heat Transf. 25(1), 137-147 (1982).

Viswanath, R. and Y. Jaluria, "A Comparison of Different Solution Methodologies for Melting and Solidification Problems in Enclosures," Numer. Heat Transf. Part B 24, 77-105 (1993).

Viswanath, R. and Y. Jaluria, "Numerical Study of Conjugate Transient Solidification in an Enclosed Region," Numer. Heat Transf. Part A 27, 519-536 (1995).

Yao, L. S. and F. F. Chen, "Effects of Natural Convection in the Melted Region Around a Heated Horizontal Cylinder," ASME J. Heat Transf. 102, 667-672 (1980).

Yao, L. S. and W. Cherney, "Transient Phase-Change Around a Horizontal Cylinder," Int. J. Heat Mass Transf. 24(12), 1971-1981 (1981).

Ritabrata Dutta, (1) Arnab Atta (2) and Tapas Kumar Dutta (2) *

(1.) Indian Statistical Institute, Kolkata 700108, India

(2.) Department of Chemical Engineering, Jadavpur University, Kolkata 700032, India

Arnab Atta's present address is Department of Chemical Engineering, IIT Delhi, New Delhi, India.

* Author to whom correspondence may be addressed. E-mail address: proftapas_dutta@yahoo.com

In many engineering applications such as thermal energy storage using phase change materials (PCMs), the melting and solidification are important phenomena. Melting of the phase-change material gives rise to natural convection, while the flow structure could significantly affect phase change process. The essential feature of melting (or solidification) of PCM is the existence of the moving interface between two phases with time. The convection influences the morphology of this solid-liquid interface. This problem was studied as early as in 1831 by Lame and Clapeyron (Huntler and Kutter, 1989). The sequence of papers written by Stefan has given his name to this type of problem as "Stefan problem." The problem of heat transfer with phase change can be formulated, considering either the temperature or the enthalpy as the dependent variable. When temperature is considered as dependent variable, the energy equations for both the phases are to be written independently and then coupled for the interface. In fact this technique requires the knowledge of the interface position explicitly for the determination of temperature, complicating the formulation and method of solution. The problem has been solved by a method of immobilization of the interface (Duda et al., 1975). When the enthalpy is considered as the dependent variable the knowledge of the interface position is not required and a single enthalpy equation for the whole domain suffices. Huntler and Kutter (1989) demonstrated the enthalpy method for heat transfer problems with moving boundaries associated with phase change. Pannu et al. (1980) calculated temperature and velocity profiles, rates of heat transfer, movements of melting interface in vertical and horizontal cylinders by finite difference method. The integral method has been applied to study the effects of natural convection on the melting of solid around the horizontal cylinder (Yao and Chen, 1980; Yao and Cherney, 1981). The perturbation and numerical solution of melting around the heated horizontal cylinder for an isothermal boundary condition were studied previously (Reiger et al., 1982; Prusa and Yao, 1984a,b) and the melting process in a rectangular enclosure driven by coupling of heat conduction in the solid phase and natural convection in the melt of the PCM has been analyzed (Bernard et al., 1986). Viswanath and Jaluria (1993, 1995) carried out a numerical modelling of the macroscopic melting processes in rectangular enclosures using an enthalpy--porosity formulation and has been used to circumvent the need to track the moving interface at every time instant. In this case a single energy equation is valid over all regions including the interface. The SIMPLE algorithm was employed to implement the solution. A numerical study of melting of PCM around a horizontal circular cylinder of constant wall temperature and in presence of natural convection in melt region was presented by Ismail and da Silva (2003a,b). In a recent literature (Lamberg et al., 2003), the comparison of experimental data with numerical results was obtained where only energy equation was solved by modelling convective source term with an effective heat capacity method. Ettouney et al. (2005) studied heat transfer characteristics of PCM in a vertical annulus and developed empirical correlations involving Nusselt number, Rayleigh number, Steffan number, and Fourier number.

Ng et al. (1998) in their study employed the finite element method to simulate the convection dominated melting of a PCM in a cylindrical horizontal annulus heated isothermally from an inside wall. The effects of Rayleigh number on the melting rate as well as the evolution of the flow pattern were examined. The authors observed that the increasing Rayleigh number promotes the heat transfer rate and the multiple cellular pattern was found to occur at high Rayleigh number (>[10.sup.6]). Khillairkar et al. (2000) studied melting in a concentric horizontal annulus of arbitrary geometric arrangement such as external square tube with circular tube inside (Type A) and circular external tube with a square tube inside (Type B). The authors studied the effects of the Rayleigh numbers as well as the heating of either the inside surface or the outside surface and in some cases both the surfaces at a temperature above the melting point of the PCM. To account for the physics of the time wise evolution of flow at the solid--liquid interface the well known enthalpy--porosity model was employed in the fixed grid method. For both the horizontal annuli of Type A and Type B, it is observed that the effects of heating both the surfaces is the same as the heating of inside surface or the outside surface separately until there is an interaction between the two melt-zone. It is observed that the melting rate is faster due to good mixing between melt zones. This suppresses the thermal stratification attained in the horizontal annuli of both Type A and B. The thermal stratification occurs in the upper part of the cavity due to the fact that the energy charged to the system is mainly carried upward by the free convection. Thus, the energy is used to raise the temperature of the melt instead of melting the PCM. So, the energy storage in the system is more in the form of sensible heat than in the form of latent heat. Balikowski and Mollendorf (2007) recently studied the performance of PCM in a horizontal annulus of a double pipe heat exchanger in a water circulating loop. PCM was placed in the annular region of a double pipe heat exchange with water circulated in the inside pipe. Experiments were performed in which the PCM would absorb or release heat at various temperatures and water flow rates. The effect of different water flow rates on the heat transfer rates were examined. The difference between temperature of hot water and melting point of the PCM was between 2 and 6[degrees]C. To avoid any thermal penalty and the flow rate of hot water were so maintained that there was virtually no difference between the inlet and outlet hot water temperatures.

In this paper, heat transfer is studied experimentally and numerically for PCM (paraffin wax) encapsulated in the horizontal annulus of two coaxial circular cylinders with isothermal and adiabatic condition at the inner and the outer cylinder, respectively. Two-dimensional unsteady Navier-Stokes' equations with energy equations were solved by finite volume method (FVM) using SIMPLE scheme to track the temperature at different points and moving solid--liquid interface numerically. Enthalpy equation was considered as the energy equation. To solve this phase change problem in a complex geometry of an annulus formed by concentric horizontal tubes the enthalpy--porosity method has been used. Comparison of time-temperature history at different points obtained from physical experimentation and numerical experimentation is made. Analysis of the resulting interaction between the flow and the heat transfer in the annular region is considered through the variation of [epsilon] (eccentricity) and the inclination angle of eccentricity ([psi]). The numerical results obtained are displayed through streamlines, isotherm contours and velocity distribution in computational domain. The variation of net circulation and analysis of the heat flux along the outer surface of the inner cylinder is done for different eccentric positions.

PHYSICAL PROCESS

In this experiment the hot water at constant temperature was passed through the inner tube. The heat is conducted through the metal wall to the surrounding wax. The innermost wall layer of wax in the annulus adjacent to the metal wall melts quickly as the temperature of water is greater than the melting point of wax, and a layer of melt surrounds the metal tube. In this situation subsequent heat transfer will take place through the layer of molten wax. The layer of molten wax is initially stagnant, and gets heated. Due to difference in density in its body an internal circulation starts for heat transfer to take place. Thus, the mechanism during melting is primarily natural convection and conduction in the liquid phase and only conduction in the solid phase. In developing the natural convection model for the horizontal annulus, the following assumptions could be made:

(1) The temperature distribution along the axis is negligible.

(2) The thermal resistance of the pipe wall is negligible.

(3) The heat loss from the outer cylinder is negligible.

The first assumption reduces the problem to a two-dimensional one, simplifying the mathematical model. The other two assumptions defined the boundary conditions.

EXPERIMENTAL CONFIGURATION AND PROCEDURE

PCM Used in the Experiment

The PCM used in the experiment is paraffin wax. Its thermophysical properties are given in Table 1.

Experimental Set-Up

The set-up (Figures 1 and 2) consists of a hollow steel cylinder with internal diameter 10.16 cm, length 30.48 cm and the internal pipe of diameter 1.905 cm is placed coaxially and horizontally.

The annular space so formed contained the PCM. The configuration is similar to a concentric pipe heat exchanger. On the body of the outer cylinder eight thermo-couple pockets, four at each side diametrically opposite to each other were constructed. Then the thermo-couples were placed in the annulus maintaining equal distance starting from the wall of the inner tube to the inner wall of the outer cylinder at equal distances along the length of the cylinder as shown in the Figure 1. The stem connected to the outer steal cylinder allows the molten wax to be filled in the annulus. When the mass of the wax in the annulus melts during the heating cycle the volume expansion of the wax is accommodated in the stem and avoids any spillage of the wax. This assembly containing PCM in the annulus was put in a Perspex box to reduce heat loss. The void space in the Perspex box was then filled with glass wool loosely allowing voids to be filled up partially by air increasing the overall insulation property of the space separating the ambient from the outer surface of the outer cylinder. Two T-joints were attached to the two ends of the inner tube for placing two more thermo-couples to measure the inlet and outlet temperatures of the flowing water.

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

Experimental Procedure

For heating, at first the water bath was set to a desired fixed temperature. When the water reached a steady state temperature hot water was passed through the inner tube at a constant rate. The temperatures of different thermo-couples placed radially in the annulus were recorded in real time in the computer by HP-bench-link software with the help of the data logger (data-acquisition unit-HP 34970 A) connected via line driver to get the heating curves for different radial points. The recording of temperature data was done at 5 s interval and thermo-couple used was T-type with accuracy 0.5-1[degrees]C. The estimated accuracy of the calibration of the flow meter and the thermo-couples is about [+ or -]5%. The resulting accuracy of the inference of the heat transfer is about [+ or -]10%. The estimated accuracy of time and length measurement is about [+ or -]0.1 s and 2-3 mm, respectively. After both the thermo-couple placed at the outermost position of the annulus reached the melting point of the paraffin wax, the experiment was terminated.

MATHEMATICAL MODEL

For the numerical simulation of melting process two methods are generally adopted. These are transformed grid (T-grid) method and enthalpy--porosity methods. In this paper, the latter method with fine meshes is employed for the numerical solution of melting of paraffin wax in the horizontal concentric annulus. The accuracy of interface tracking has been achieved by considering finer meshes for this complex geometry.

Governing Equations and Enthalpy Formulation

In the enthalpy--porosity formulation, enthalpy has been taken as the dependent variable in place of temperature, along with other primitive variables u, v, and P. In this scheme the liquid mass fraction (f), defined as the ratio of the liquid mass to the total mass in a given computational cell plays an important role. Then the enthalpy content of that computational cell will be:

H = fL + [C.sub.p]T (1)

where the enthalpy content of the solid has been taken as reference enthalpy. Here the heat capacity ([C.sub.p]) can vary with the change of phases.

Now, in other way round the liquid mass fraction can be obtained from the enthalpy in the following way:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

The continuity and momentum equations for the dependent variable u, v, P, considering Boussinesq' approximation, are the following:

[partial derivative]u/[partial derivative]x + [partial derivative/[partial derivative = 0 (3)

[rho] ([partial derivative]u/[partial derivative]t + u [partial derivative]u/[partial derivative]x + v [partial derivative]u/[partial derivative]y) = [partial derivative]p/[partial derivative]x + [mu] ([[partial derivative].sup.2]u/[partial derivative][x.sup.2] + [[partial derivative].sup.2]u/[partial derivative][y.sup.2]) + [S.sub.u] (4)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

and the energy equation is:

[rho] ([partial derivative](cT)/[partial derivative]t + u [partial derivative](cT)/[partial derivative]x + v [partial derivative](cT)/[partial derivative]y) = k ([[partial derivative].sup.2]T/[partial derivative][x.sup.2] + [[partial derivative].sup.2]T/[partial derivative][y.sup.2] - [rho]L [partial derivative]f/[partial derivative]t (6)

Source Terms

In the Equations (4) and (5) the source terms [S.sub.u] and [S.sub.v] have been chosen in a manner such that they can suppress the velocity fields in the solid phase. When the PCM is in solid phase it takes large values, suppressing the velocity field. Also, with the increase of 'f' it should go smaller and finally be 0. One of the most popular and successful methods has been given by Viswanath and Jaluria (1993, 1995). They formulated the source term from the Karman-Kozeny equation for flow in porous media, where the equations resemble a Darcy-type law (where the velocities are proportional to the pressure gradient). Following are the source terms:

[S.sub.u] = -C (1-[f.sup.2]/[f.sup.3]+b) and [S.sub.v] = -C (1-[f.sup.2]/[f.sup.3]+b) v

The constants C and b can be suitably chosen to simulate the flow characteristics in porous media. From the equations it is clearly evident that when the liquid fraction is 0 then the source term dominate the equation and force the velocities to become zero in the region. The source term becomes negligible when liquid fraction is 1.

Boundary Conditions and Solution Method

In conformity with the experimental conditions the following boundary conditions have been imposed.

(1) No-slip condition at the surfaces of the tubes (inner and outer), forming the annulus.

(2) Isothermal condition at the inner surface of the annulus and adiabatic condition at the outer surface of the annulus.

Solution of the above mentioned Equations ((3)-(6)), after discretising by FVM, have been done by the SIMPLE algorithm (Patankar, 1980), introducing power law scheme and the deferred correction method (Ferziger and Peric, 2002) in the interpolation scheme for the convection term. After proper solution of these equations at each time step, the enthalpy and liquid mass fraction have been updated by the Equations ((1)-(2)).

Code Validation

Our code has been validated with the benchmark problem of natural convection in a square enclosure with a hot and cold wall. Comparisons were made with the previous works on phase-change materials, specifically the vertical annulus case. In all those cases the results were in good agreement.

Time-Efficiency of the Code

The code is also efficient in the sense of time-efficiency. Calculation for our problem at a grid-size of 60 x 120, in a P4 workstation with 512 MB RAM have a ratio of experiment time and numerical simulation time of 20:1.

NUMERICAL EXPERIMENTATION

For a better insight in the natural convection in phase change heat transfer in an annulus, the effects of the eccentricity as well as the angle of inclination of eccentricity are considered to be of great importance. It is already observed that the effects of Rayleigh number and Prandtl number on the natural convective heat transfer is well established. In this study the eccentricity is varied from 0 to +0.5. Similar such positions of eccentricity is varied among the different angle of inclination ([psi]) between [pi]/2 to [pi]/2 (e.g. -90[degrees], -60[degrees], -30[degrees], 0[degrees], 30[degrees], 60[degrees], 90[degrees]) as shown in Figure 3.

[FIGURE 3 OMITTED]

Due to the eccentric position the annulus becomes enlarged as well as narrowed in the direction of angle of inclination giving rise to the different fluid dynamic phenomena. Eminent thermal effects are also observed in both the parts of the annulus across the vertical midplain along the centre of the inner cylinder. In the numerical experiments the dimensions of the annulus assumed to be of 1 m.(i.d.) and 5 m (o.d.), which are larger than that of experimental set-up. The assumed physical properties of the material filled in the annulus are given in Table 2.

The boundary conditions of the inner and outer surfaces of the annulus are assumed to be isothermal and adiabatic, respectively. The above numerical experimental conditions are not varied as the dependence of heat transfer and fluid dynamics on Rayleigh and Prandtl numbers are well established.

The angular positions on the surface of the inner cylinder are shown in Figure 4. Due to the formation of thermal plume by the natural convective flow in the melt phase, the dimensionless heat flux value changes from the bottom point to the top point on the surface of the inner cylinder. The dimensionless heat flux is found by, -d[THETA]/dR at the surface where [THETA] = [T.sub.H]-T/[T.sub.H]-[T.sub.M], R = r/[r.sub.o] -[r.sub.i], [T.sub.H] = temperature of the inner surface, T = temperature of any point in the annulus, [T.sub.M] = melting temperature of the PCM, r = radius of any point on the annulus, [r.sub.i] = inner radius, [r.sub.o] = outer radius of the annulus.

[FIGURE 4 OMITTED]

The numerical simulation is done for a number of grid sizes in both the radial and circular directions. The mesh size is chosen after the performance of several tests as 200 x 200. For the accuracy of the numerical results and the convergence the error level is taken to be [10.sup.-5]. All the simulations are carried up to a time 30 000 s ([F.sub.0] = 3.08).

RESULTS AND DISCUSSION

The real time temperatures at different points away from the outer surface of the inner tube, that is, at four cross-sectional points in upper and lower position of the annulus, were measured. The position of the censors are clearly marked in Figure 2, which is a cross-sectional view of the horizontal cylindrical annulus. The positions of the temperature censors attached to the outer surface of the inner tube at an equal distance from the two ends along the length of the cylinder are shown in Figure 1. As temperature difference of outlet and inlet water was very small (<1[degrees]), it is assumed that the heating takes place in isothermal condition. The temperature difference between these two point at similar radial position but at different longitudinal points are very small and conforms the assumption 1 of the Physical Process section. The experiment was performed for different temperature of water (75, 80, 85, and 90[degrees]C). From the experiment the time--temperature history for eight points was obtained. The data for 15 000 s or more were recorded till the temperature of the outermost node reached the melting point. One representative graph of time-temperature history for hot water at 80[degrees]C at the points U3 and L3 is given below in Figure 5.

The time--temperature history for the wax in the upper and lower half of the annulus shows a substantial difference. The "time--temperature history" of a point in the upper half of the annulus is quite different from the same point in the lower half having same radial position. This is because as soon as the hot fluid is allowed to flow through the inner tube, the tube wall reaches a higher temperature than the melting point of the wax in the annulus (in all our experimental cases) forming a molten layer of wax around the inner tube. The molten layer grows, receiving heat from the hot fluid flowing inside the inner tube. Due to natural convection, buoyancy-driven upward flow of the molten wax enhances the supply of heat to the wax in the upper part of the annulus and solid--liquid interface moves faster in the upper half of the annulus. The experimental result supports this fact by showing a higher growth rate of melting front in the upper half of the annulus. It is evident from the previous discussion and the experimental result that the heat transfer phenomena here consist of both the convective and conductive mode of energy transfer. To explain the physical phenomena, time--temperature history inside the annulus for 15 000 s was simulated, by our mathematical model. Other variables like velocities, pressure and interface positionwere also being tracked for this time-span. Velocity vector profiles at three different time points 5000, 10 000, and 15 000 s by simulated data are given in Figure 6, in which only one-half of the annulus is shown, as both the halves are symmetrical showing similar features. White parts show the solid parts of the paraffin wax and the rest is molten portion having different velocities in different part with the progress of the melting.

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

The comparisons of experimental and simulated time-temperature history are given in the Figures 7 and 8. The figures indicate a good agreement between experimental and simulated data in these cases. The discrepancies may be due to sudden dislodging of a piece or pieces of semi-molten wax in the molten wax. In the positions nearer to the outer wall (U2) the simulated temperature is higher than the experimentally predicted one. This may be due to heat loss from the outer cylindrical surface even with the insulation provided around. Similar accuracy in prediction by simulated data of this order was prevalent in the cases of other hot fluid temperatures (80, 85, and 90[degrees]C).

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

Effects of Eccentricity

The effect of eccentricity [epsilon] on the kinematic and thermal field is studied for Ra = [5.7710.sup.7], Pr = 46.5 for the inclination angles [phi] =, -90[degrees], -60[degrees], -30[degrees], 0[degrees], 30[degrees], 60[degrees], 90[degrees] and [epsilon] = 0.0, 0.25, 0.50.

For the concentric annulus the streamlines and the isotherms are kidney shaped or crescent shaped and symmetric along the vertical midline along the centre of the inner cylinder. On both sides of the inner cylinder, the fluid flows up along the heated inner cylinder wall and down along the melt front. At a time so far elapsed with the fluids in the central core vortices are almost stagnant, and both sides have equal maximum stream function values, that is, [[psi].sub.max].

For an eccentric annulus, the flow pattern is also generally kidney shaped but is asymmetric along the midline. The position and formation of the crescent on two sides change with respect to the eccentricity [epsilon]. As the inner cylinder moves from the centre to the right side the extent of the left kidney gradually increases as the gap increases. The centre of the closed loop moved upward in the annulus. But the right crescent gradually reduces in size until the two crescent shaped flow region appears in both the left and right sections of the annulus when the eccentricity reaches a critical value. It is also observed that the eccentricity generally increases the strength of the flow circulation, that is, the value of the maximum stream function on the "intrusion" side of the flow region and decreases the flow strength of the "expanded" side of the flow. The imbalances of the flow strength between the two sides lead to a net circulation around the eccentric annulus excepting in the vertical eccentric cases.

[FIGURE 9 OMITTED]

The primary attention will be focused upon a problem of heat transfer by convection in a system of coaxial cylinder at fixed temperatures leading to approximately double increase of a flow compared with that determined only by heat conductivity. In the concentric case the strength of buoyancy-driven fluid motion as well as the enhancement in the convective heat transfer with the increase in Rayleigh number are well established. As the Rayleigh number becomes sufficiently high the inner cylinder is surrounded by the thermal boundary layer which starts at the bottom of the cylinder and leaves at the top. The inner cylinder local heat transfer coefficient varies with the angular position as the boundary layer develops from its bottom and is shed as a plume at the top. This plume heats the outer surface of the melt, causing the heat transfer coefficient at that position to be much higher than that for the inner one at the top. As the flow moves down along the outer surface of the melt, there decreases very rapidly. It is also found that the thermal stratification in the lower zone of the annulus weakens the fluid motion in that zone. One single recirculating vortex is found in a half domain and the strength of the flow obviously increases with the Rayleigh number. The effects of the eccentricity and its inclination angle on the heat transfer is considered with a fixed Rayleigh number. The value of the dimensionless eccentricity is varied from 0 to +0.5 and the inclination angle of inclination is varied from -90[degrees] to 90[degrees] by an interval of 30[degrees], by fixing [epsilon]. For a fixed inclination angle at the vertical axis of symmetry that is [phi] = -90[degrees] or 90[degrees], when the inner cylinder moves downward it causes a stronger fluid motion in the upper zone and also increases the temperature gradient in the lower zone, resulting in higher heat transfer between two cylinders, Figure 9 in comparison to the concentric case. On the other hand when the inner cylinder moving upward provides least favourable circumstances for the development of natural convection (Figure 10). Here the concentric cell centre has moved towards [phi] = 90[degrees]. The circulation cell is less powerful than in the concentric case. Clearly the top narrow mid gap between the two cylinders inhibits the convective motion in the fluid. Thus, the isotherms shows maller temperature increase in comparison to the concentric case. The bottom of the annulus near [phi] = -90[degrees] exhibits an enlarged stagnant region of fluid. In each case the motion of the fluid is clockwise on the right-hand side of the annulus. The region of almost stagnant fluid exists near the lower zone and this zone increases when the inner cylinder moves upward along the vertical axis increasing resistance to flow. Since the eccentricity due to the position of the inner cylinder in the upper side of the horizontal axis provides the least favoured circumstances for the development of natural convection, both the size and strength of flow of momentum are markedly reduced.

[FIGURE 10 OMITTED]

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

The heat flux for the natural convective heat transfer from the surface of the inner cylinder to the molten PCM changes along the surface from the bottom to the top point ([theta] = 0[degrees] -180[degrees]). The heat flux has the maximum value at the bottom and decreases to a minimum value at the top point, that is almost at 180[degrees] (e.g. Figures 11-14). In the case of concentric annulus, as there is virtually no net circulation, the variation of heat flux with position on the surface of the inner cylinder is almost symmetric about the vertical midline. It is observed that with the change in the vertical eccentricity the heat flux value reaches the maximum for [epsilon] = 0.5 in the direction of [phi] = -90[degrees]. On the other hand at the same value of the eccentricity in the direction of [phi] = 90[degrees] the heat flux at the bottom point on the surface of the inner cylinder reaches a very low value. The changes in the angle of inclination of the eccentricity the heat flux variation along the surface of the inner cylinder are not symmetric about the vertical midline due to the presence of net circulation. For a particular angle of inclination when the inner cylinder is in the first and the fourth quadrant the maximum thermal flux value are different. The maximum thermal flux is higher in the case of first quadrant than that in fourth quadrant. It is also observed that there is a shift of angular position for the occurrence of the minimum heat flux for the cases of same inclination angle but at different angle (Figures 12-14).

[FIGURE 13 OMITTED]

[FIGURE 14 OMITTED]

Effects of Angle of Inclination of Eccentricity

The net circulation [[psi].sub.0], being the difference between the [[psi].sub.max] values of the "intrusion" and "expanded" zones of the annulus, changes with the inclination angle of the eccentricity. When the inclination angle is -[pi]/2 or [pi]/2, that is, the case of vertical eccentricity, the net circulation is absent or very small. The net circulation is also symmetric along the axis of [phi] = -30[degrees] and has a maximum value. This also indicates the imbalance of fluid convection between the two sides of the inner cylinder is the reason for the net circulation around the annulus. A new phenomenon identified by this study is that when the inclination angle of eccentricity is not vertical, a net circulation is identified around the inner cylinder. For a given [epsilon] the magnitude of net circulation is expected to approach maximum when the angle of inclination is horizontal. The strength of the net circulation reaches a maximum for a given cylindrical annuli with the eccentricity in the range of 0 < [epsilon] < + 0.5, for the angle of eccentricity [phi] = -30[degrees] (Figure 15). This is due to the fact that the natural circulation of the molten phase creates an unmelted or solid phase front as an eccentric zone around the molten phase in addition to its eccentric position with respect to outer cylinder. The maximum value of the dimensionless heat flux changes with the eccentricity as well as the angle of inclination of eccentricity (Figure 16). It is observed that the maximum value of dimensionless heat flux reaches the highest value at [epsilon] = 0.5 and [phi] = 90[degrees].

[FIGURE 15 OMITTED]

[FIGURE 16 OMITTED]

CONCLUSIONS

A computational model for the prediction of the thermal behaviour of paraffin wax during its melting within an annulus formed between two horizontal cylinders has been presented. The model rests on solving Navier-Stokes equation with the consideration of the isothermal condition in the inner core and adiabatic condition at the outer surface of the annulus. The velocity vector field and the isotherms obtained from the computation confirm the conductive and conductive-convective heat transfer phenomena existing in the solid and the molten phases, respectively. The performance of the present model is verified with the experimental findings. The predicted result shows good agreement with the experimental results.

The results of the mathematical experiment obtained by varying both the eccentricity as well as their angle of inclination identify some new phenomena of flow. When the inclination angle of eccentricity is not in the vertical direction a net circulation is identified around the surface of the inner cylinder. For a given eccentricity the magnitude of the net circulation approaches maximum value when its angle of inclination approaches -[degrees]30? (Figure 15).

The analysis of heat flux along the surface of the inner cylinder indicates that for a fixed angle of inclination heat flux reaches a maximum value for the highest eccentricity ([epsilon] = 0.5) in the fourth quadrant and reaches a minimum value when the eccentricity is the maximum in the first quadrant. Similar trends are followed for all other angle of inclinations (Figures 11-14).

The maximum heat flux value increases monotonically both with angle of inclination and eccentricity (Figure 16).

NOMENCLATURE [C.sub.p] specific heat at constant pressure E distance between centres of two cylinders [F.sub.0] = t[alpha]/ Fourier number 4[r.sup.2.sub.i] f liquid mass fraction g gravitational acceleration H enthalpy k thermal conductivity of PCM L latent heat of fusion of PCM r' = [r.sub.o] - [r.sub.i] characteristics length Nu Nusselt number p pressure Pr = [C.sub.p][mu]/k Prandtl number R = r/[r.sub.o] -[r.sub.i] radial coordinate [r.sub.o] outer radius [r.sub.i] radius of Inner cylinder Ra = [C[rho]g[beta][r'.sup.3] Rayleigh number [DELTA]T/k[upsilon u, v velocities of x and y coordinate t time T temperature [T.sub.M] melting temperature of PCM [T.sub.H] temperature of the inner cylinder [alpha] = k/[rho][C.sub.p] thermal diffusibility [beta] thermal expansion coefficient T - [T.sub.M] temperature difference [THETA] = T - [T.sub.M]/ non-dimensional temperature [T.sub.H] - [T.sub.M] [psi] stream function [epsilon] = E/r' eccentricity [phi] inclination angle of eccentricity [upsilon] = [mu]/[rho] kinematic viscosity [rho] density [mu] viscosity [S.sub.u], [S.sub.v] source term b, c constant -d[THETA]/dR dimensionless heat flux

Manuscript received March 25, 2007; revised manuscript received October 17, 2007; accepted for publication January 13, 2008.

REFERENCES

Balikowski, J. R. and J. C. Mollendorf, "Performance of Phase Change Materials in a Horizontal Annulus of a Double--Pipe Heat Exchanger in a Water Circulating Loop," ASME J. Heat Transf. 129, 265-272 (2007).

Bernard, C., D. Gobin and A. Zanoli, "Moving Boundary Problem: Heat Conduction in the Solid Phase of a Phase-Change Material during Melting Driven by Natural Convection in the Liquid," Int. J. Heat & Mass Transf. 29, 1669-1681 (1986).

Duda, J. L., M. F. Malone, R. H. Notter and J. S. Ventras, "Analysis of Two-Dimensional Diffusion-Controlled Moving Boundary Problems," Int. J. Heat Mass Transf. 18(7), 901-910 (1975).

Ettouney, H., H. El-Dessouky and A. Al-Ali, "Heat Transfer During Phase Change of Paraffin Wax Stored in Spherical Shells," J. Solar Energy Eng. 127, 357-365 (2005).

Ferziger, J. H. and M. Peric, "Computational Methods for Fluid Dynamics," Springer-Veralg Inc, (2002).

Huntler, L. W. and J. R. Kutter, "The Enthalpy Method for Heat Conduction Problems with Moving Boundaries," ASME J. Heat Transf. 111, 239-242 (1989).

Ismail, K. A. R. and M. D. E. da Silva, "Melting of PCM Around a Horizontal Cylinder With Constant Surface Temperature," Int. J. Therm. Sci. 42, 1145-1152 (2003a).

Ismail, K. A. R. and M. D. E. da Silva, "Numerical Solution of the Phase Change Problem Around a Horizontal Cylinder in the Presence of Natural Convection in the Melt Region," Int. J. Heat Mass Transf. 46, 1791-1799 (2003b).

Khillairkar, D. B., Z. X. Gong and A. S. Majumder, "Melting of Phase Change Material in a Concentric Horizontal Annuli of Arbitrary Cross Section," Appl. Therm. Eng. 20, 893-912 (2000).

Lamberg, P., R. Lehtiniemi and A. M. Henell, "Numerical and Experimental Investigation of Melting and Freezing Processes in Phase Change Material Storage," Int. J. Therm. Sci. (2003).

Ng, K. W., Z. X. Gong and A. S. Majumder, "Heat Transfer Convection-Dominated Melting of a Phase Change Material in a Horizontal Annulus," Int. Commun. Heat Mass Transf. 25(5), 634-640 (1998).

Pannu, J., G. Joglekar and P. A. Rice, "Natural Convection Heat Transfer to Cylinders of Phase Change Material Used for Thermal Torage," AIChE Symp. Ser. 76, 47-55 (1980).

Patankar, S. V., "Numerical Heat Transfer and Fluid Flow," Hemisphere, Washington, D.C. (1980).

Prusa, J. and L. S. Yao, "Melting Around a Horizontal Heated Cylinder: Part-I Perturbation and Numerical Solutions for Constant Heat Flux Boundary Condition," ASME J. Heat Transf. 106(2), 376-384 (1984a).

Prusa, J. and L. S. Yao, "Melting Around a Horizontal Heated Cylinder: Part-II Numerical Solution for Isothermal Boundary Condition," ASME J. Heat Transf. 106(2), 469-472 (1984b).

Reiger, H., U. Projan and H. Beer, "Analysis of the Heat Transport Mechanisms during Melting Around a Horizontal Circular Cylinder," Int. J. Heat Transf. 25(1), 137-147 (1982).

Viswanath, R. and Y. Jaluria, "A Comparison of Different Solution Methodologies for Melting and Solidification Problems in Enclosures," Numer. Heat Transf. Part B 24, 77-105 (1993).

Viswanath, R. and Y. Jaluria, "Numerical Study of Conjugate Transient Solidification in an Enclosed Region," Numer. Heat Transf. Part A 27, 519-536 (1995).

Yao, L. S. and F. F. Chen, "Effects of Natural Convection in the Melted Region Around a Heated Horizontal Cylinder," ASME J. Heat Transf. 102, 667-672 (1980).

Yao, L. S. and W. Cherney, "Transient Phase-Change Around a Horizontal Cylinder," Int. J. Heat Mass Transf. 24(12), 1971-1981 (1981).

Ritabrata Dutta, (1) Arnab Atta (2) and Tapas Kumar Dutta (2) *

(1.) Indian Statistical Institute, Kolkata 700108, India

(2.) Department of Chemical Engineering, Jadavpur University, Kolkata 700032, India

Arnab Atta's present address is Department of Chemical Engineering, IIT Delhi, New Delhi, India.

* Author to whom correspondence may be addressed. E-mail address: proftapas_dutta@yahoo.com

Table 1. Thermo-physical properties of paraffin wax Melting point, [T.sub.m] 42[degrees]C Specific heat of solid (40[degrees]C), [C.sub.ps] 2.0 kJ/kg K Specific heat of liquid (60[degrees]C), [C.sub.pl] 2.15 kJ/kg K Latent heat of melting, L 30 kJ/kg Thermal conductivity of solid, [K.sub.s] 0.173 W/m K Thermal conductivity of liquid, [K.sub.l] 0.08 W/m K Density of solid (at 40[degrees]C), [[rho].sub.s] 857 kg/[m.sup.3] Density of liquid (at 60[degrees]C), [[rho].sub.l] 790 kg/[m.sup.3] Viscosity, [mu] 0.78 N s/[m.sup.2] Thermal expansion coefficient, [beta] 0.001 Table 2. Thermo-physical properties of paraffin wax Melting point, [T.sub.m] 27.4[degrees]C Specific heat of solid (29[degrees]C), [C.sub.ps] 1.97 kJ/kg K Specific heat of liquid (25[degrees]C), [C.sub.pl] 2.15 kJ/kg K Latent heat of melting, L 228 kJ/kg Thermal conductivity of solid, [K.sub.s] 0.081 W/m K Thermal conductivity of liquid, [K.sub.l] 0.08 W/m K Density of solid (at 40[degrees]C), [[rho].sub.s] 814 kg/[m.sup.3] Density of liquid (at 60[degrees]C), [[rho].sub.l] 775.5 kg/[m.sup.3] Viscosity, [mu] 3.7 Ns/[m.sup.2] Thermal expansion coefficient, [beta] 0.003

Printer friendly Cite/link Email Feedback | |

Author: | Dutta, Ritabrata; Atta, Arnab; Dutta, Kumar Tapas |
---|---|

Publication: | Canadian Journal of Chemical Engineering |

Date: | Aug 1, 2008 |

Words: | 6409 |

Previous Article: | Laminar, transitional and turbulent flow of Herschel--Bulkley fluids in concentric annulus. |

Next Article: | On the drift-flux analysis of flotation and foam fractionation processes. |