A Finite Volume Simulation of Crystal Growth in Water Falling Films over a Cold Isothermal Surface

A Finite Volume Simulation of Crystal Growth in Water Falling Films over a Cold Isothermal Surface

Abdulmajed Khalifa

Mechanical Engineering Department University of Zawia

Zawia, Libya

Ramadan Gennish

Mechanical Engineering Department University of Zawia

Zawia, Libya

Shaban Jolgam

Mechanical Engineering Department University of Zawia

Zawia, Libya

Ragab Khalil

Mechanical Engineering Department University of Benha

Benha, Egypt

Mussa Radwan

Mechanical Engineering Department University of Tripoli

Tripoli, Libya

AbstractA numerical model of ice crystal growth in a cooled laminar falling film has been developed using a control volume technique. A source term has been included in the energy equation to represent the phase change energy. The effective values of thermal conductivity, viscosity, thermal diffusivity, and specific heat are considered as functions of the volumetric concentration of ice crystals in the falling film. Film temperature distributions and ice crystal diameters in the film are calculated numerically. It is found that the wall Nusselt number is a function of film thickness and is essentially independent of the Stefan number. The ice crystal growth rate was found to be primarily a function of the Stefan number and film thickness. The ice crystal Nusselt number and super-cooling of the solution are dependent on the Stefan number, while the Reynolds number (film thickness) has a lower influence on these parameters.

KeywordsHeat transfer enhancement; ice slurry; cooling; laminar flow; ice crystal growth

1. INTRODUCTION

   Thermal energy storage (TES) plays a significant role in shifting cooling loads in buildings and consequently reduces energy costs due to the reduction in electrical power demand in peak hours. Amongst the various thermal energy storage techniques, latent heat thermal storage (LHTS) is particularly attractive due to its ability to provide a high energy storage density at a constant temperature corresponding to the phase transition temperature of the heat storage material. Ice slurries are interesting secondary refrigerants compared to single-phase fluids since they produce high heat capacities due to latent heat storage. An important advantage of this high heat capacity is the possibility of cold storage that can take place either as a static process, in which heat transfer takes place via a solid surface, or a dynamic process. In the dynamic process (using ice slurry), a large amount of cold energy can be transported

   with less pumping work due to fluidity. Also, the dynamic system responds quickly to the change in heat load because the ice particles have a large surface area. The formation of ice slurries and their removal phenomena on a horizontal surface has been studied by many researchers [1-3]. The researchers in

   1. predicted that the change in size and shape is supposed to be the explanation for the difference in pressure drops in tubes of freshly produced ice slurry and those of the same ice slurry after storage. In [5], a numerical model was developed to approximate the growth rate of ice crystals in a falling film flowing down a cooled vertical plate. A source term was included in the energy equation to represent the phase-change energy. The governing equations were solved using the finite difference technique. The growth rate of ice crystals, Nusselt numbers, and overall heat transfer coefficients between the film and cooled plate were calculated for various film thicknesses, Stefan numbers (Ste), and initial ice crystal diameters. They found that the growth rate of ice crystals and the Nusselt number depend on film thickness, initial ice number concentration, and Stefan number. Their results of heat transfer coefficients were in reasonable agreement with experimental data. Pronk et al. in [6] studied the time-dependent behavior of different ice slurries during storage. They investigated mechanisms such as attrition, agglomeration, and Ostwald ripening that cause changes in the crystal size distribution, which have great importance for ice slurry systems since they influence other parameters such as pressure drop and the characteristics of heat transfer. Kozawa et al. in [7] investigated the effects of ice content (ice packing factor, IPF) and mass flow rate of supplying ice-slurry on the ice-storing characteristics in a tank. It was found that by raising IPF and reducing the mass flow rate of supplying ice-slurry, the ice-rich layer was hard to spread horizontally in the tank. Cliche and

      Lacroix in [8] presented a mathematical model for the optimization of ice formation in a laminar falling film. Numerical experiments were conducted to study the effects of the mass flow rate, the inlet water temperature, and the length and temperature of the cold plate on ice accumulation.

      An ice slurry with enhancing additives was studied as a secondary refrigerant in cold storage systems using a water/oil (W/O) type emulsion [9-14]. It was concluded that a high- performance ice slurry with low adhesion to the cooled wall could be formed by the water-oil type emulsion. Matsumoto et al. in [15] reported that a high ice packing factor (IPF) ice slurry could be formed without adhesion of ice to the cooling wall by cooling and stirring a functional fluid of 10 vol% silicone oil and 90 vol% water with a small amount of additive (silane- coupler) in a resin beaker. Matsumoto et al. in [9] carried out experiments to clarify the characteristics of the suspension formation process in an ice storage system using a water/oil mixture. From the thermal analysis of the formed substance, it was found that the substance was not ice but a compound of ice and an additive with a reduction in the freezing point about 7 oC below 0 oC. Swanson and Nelson in [16] developed new equipment that can precisely measure the growth rates of ice crystals and droplets at temperatures as low as -60 oC. Using this equipment, preliminary studies have indicated the advantage of following individual faces of multiple crystals. They improved significantly thermodynamic controls than it was with previous instruments, where they were able to grow multiple crystals under identical thermodynamic conditions.

      The present work is focused on the growth of phase change

      • The density of the fluid and the ice crystals are assumed to be constant and equal.

        Radiant and evaporative heat losses were neglected while convective heat transfer was allowed at the film-free surface.

        Based on the assumptions mentioned above, the plate is maintained at a constant temperature that is lower than the freezing temperature of the solution. The heat is transferred as a latent heat from the suspended ice crystals to the surrounded super-cooled fluid and then to the plate surface. As a result of this mechanism of pumping latent heat to the surrounding fluid, more ice formed (increasing ice crystal size).

        Using a source term to simulate the phase change, the governing steady-state energy equation can be expressed as follows:

        = ( ) +

        (1)

        The flow is assumed to be fully developed, and the velocity profile suggested by Levich in [18] is applied in the present model as follows:

        = (1 ) 2

        (2)

        The source term s in (1) is due to the latent heat released during the Solidification process. The physical model is exposed to the following boundary conditions:

        (, 0) = ; (0, ) = = ;

        (,)

        materials, in particular ice crystals, to produce an enhanced heat transfer medium in a dynamic cold storage system. The adopted

        |

        =

        = ((, ))

        finite volume of this study combines some of the numerical models that have been devoted to laminar fully developed falling films. The ice crystal temperature and concentration are predicated. The influences of the Stefan number and film thickness on the wall and ice crystal Nusselt numbers, ice crystal growth rate, and super-cooling of the solution are investigated.

2. MATHEMATICAL MODELING

   Figure 1 shows a schematic diagram of the falling film physical model. The film is considered as an ice slurry flowing along with a cooled isothermal vertical plate. The following assumptions are used to approximate the falling film model:

   • The film is smooth, laminar, incompressible, and hydro-dynamically fully developed.

   • The flow is Newtonian, which is approximately the case for an ice crystal concentration of less than 45% [17].

   • The suspended ice crystals in the falling film are assumed to be spheres with no-slip conditions between the ice particles and the surrounding fluid.

   • The physical properties are assumed constant except for effective values of thermal conductivity, viscosity, and specific heat due to ice and fluid interactions.

   • There is no physical interaction between the plate and the ice crystals.

   • The ice particles are assumed to be homogeneously distributed in the flow field.

   As suggested by Vand in [19], the bulk dynamic viscosity of the ice slurries is assumed to be a function of ice crystal volumetric concentration C and expressed as follows:

   = (1 1.182)2.5 ; 0.45

   (3)

   and the effective kinematic viscosity is defined as follows:

   =

   (4)

   The effective specific heat of the slurries is the averaged mass specific heats of the fluid and ice particles, which is defined by Stewart and Stickler in [20] as follows:

   = () + (1 )

   (5)

   where and are the specific heats of the ice particles and fluid, respectively. Also, the effective density of the slurries is the averaged mass densities of the fluid and ice particles

   , which is defined as follows:

   = () + (1 )

   (6)

   where is the volumetric thermal expansion coefficient, assumed to be equal to 1/ T.

   The ice crystal growth throughout the flow field can be obtained from (10) and (11) as:

   2( ) = if

   (14)

   Fig. 1: Schematic diagram of falling film model.

   The bulk thermal conductivity of static dilute suspension

   The volumetric concentration C is the number of ice crystals per unit volume n times the volume of each ice crystal

   3 = 6

   (15)

   3 = 6

   (15)

   crystal as follows:

   can be evaluated from Maxwell's relation, presented by Maxwell in [21] as follows:

   2 + / + 2( 1) = 2 + ( 1)

   (7)

   Note that, the subscripts p and f represent ice particle and fluid, respectively. Due to the interaction between ice particles and the fluid, the effective conductivity of the slurry flow is higher than that predicted by (7). Reference [22] proposed a correlation for the effective thermal conductivity in flowing conditions as:

   For ice slurry, the freezing temperature is a polynomial function of the ice crystal concentration discussed by Burns et al. in [25] as;

   = 0 + 1 + 22

   (16)

   0 = 271.0, 1 = 1.492 and 2 = 2.38

   (17)

   0 = 271.0, 1 = 1.492 and 2 = 2.38

   (17)

   The coefficients 0, 1and 2 are calculated from the phase diagram of NaC1-2O system presented by Fang et al. in [26] as follows:

   where p

   2

   2

   = 1 + ; where d is the ice crystal dia

   (8) meter, =

   = 1 + ; where d is the ice crystal dia

   (8) meter, =

   =

   = , = , = , = , = () ()

   (18)

   = , = , = , = , = () ()

   (18)

   e

   Finally, by defining the following dimensionless variables:

   and e

   = . The experimental results presented by Sohn and

   Chen in [23] were used to evaluate the constants B and m at moderate Peclet numbers and assigned them values of 1.8 and 0.18, respectively. The heat generation rate per ice crystal is given as:

   where L is the length of the plate, t is the film thickness and is the maximum fluid velocity which occurs at the falling film free surface (y=t) and defined as:

   crystal = if

   (9)

   crystal = if

   (9)

   with the help of (18) the dim

   2 = 2

   (19)

   2 = 2

   (19)

   ionless mathematical model

   where hif is the latent heat of the liquid solvent.

   = = () 2 if 2

   (10)

   = = () 2 if 2

   (10)

   Then the source term s in the governing energy equation can be computed from:

   becomes as follows:

   () = ( ) +

   (20)

   () = ( ) +

   (20)

   where:

   ens

   (, 0) = 0; (0, ) = 1 and ( ) = (=1 ) =1	(21a)
   2 = ( ) and = ( )	(21b)
   = 2 2 ; = and 2 = ( )	(21c)
   3 2 = [1 + ( (1 ) ) ] e	(21d)

   (, 0) = 0; (0, ) = 1 and ( ) = (=1 ) =1	(21a)
   2 = ( ) and = ( )	(21b)
   = 2 2 ; = and 2 = ( )	(21c)
   3 2 = [1 + ( (1 ) ) ]/p> e	(21d)

   where n is the number of ice crystals per unit volume.

   Due to the convective heat transfer between the ice crystals and the fluid, the phase change source term can also be written as:

   = ( )

   (11)

   where is the surface area of each ice crystal (2). The heat

   transfer coefficient

   between the ice crystal and the fluid is

   determined from natural convection correlation for spheres as the particle Nusselt number, which is presented by Stewart et al. in [24] as follows:

   = [2 + 0.50.25]

   (12)

   = [2 + 0.50.25]

   (12)

   where Rayleigh number Ra is defined as follows:

   ( )3 = e

   (13)

   = and 3 = ( ) 2 2	(21e)
   ( ) = and if = [2 + 0.50.25]	(21f)
   = (1 1.182)2.5 and = + (1 )	(21g)
   0 1 = + 2 + 2 and =	(21h)

   = and 3 = ( ) 2 2	(21e)
   ( ) = and if = [2 + 0.50.25]	(21f)
   = (1 1.182)2.5 and = + (1 )	(21g)
   0 1 = + 2 + 2 and =	(21h)

   Note that the operator , returns the greater of A and

   B. The neighbors of point P are assigned e, w, n, and s

   corresponding to the east, west, north, and south, respectively.

   The discretized form of the boundary condition at the free surface is:

   where

   =

   = + (e), = and

   (25) = +

   = + (e), = and

   (25) = +

   ()

3. NUMERICAL METHOD

   The governing partial differential equations have been discretized on a uniform grid using a control volume formulation [27]. Equations (20) and (21b) are solved to obtain the dimensionless temperature () and ice crystal diameter (D) distribution throughout the flow domain at steady-state conditions, respectively.

   = + + + +

   (22)

   = + + + +

   (22)

   = (||) and = (| |)

   (23a)

   = (||) and = (| |)

   (23a)

   The discretized form of the energy equation is: where the coefficients are defined as follows:

   = , 0 and = , 0

   (23b)

   = and = + + +

   (23c)

   = () and = ()

   (23d)

   () () + 1 = [ + ] and () () + 1 = [ + ]

   (23e)

   (| |) = 0, (1 0.1| |)5, and = , , and

   (23f)

   The Peclet numbers and the diffusion coefficients, ,

   , are defined respectively as follows:

   = , = and =

   (24)

   The discretized form of the boundary condition at the outflow boundary can be written as follows:

   = + + +

   (26)

   The above algebraic equations, allow the iterative solution for the two unknowns and D. The numerical computations have been performed on a uniform 50×50 grid where a of

   0.02 mm and =2 mm. The solution is converged when the changes in s and Ds at all grid points fall below convergence criteria of 10-6.

   The inlet flow at the top of the plate is assumed to contain a specified uniform number of ice crystals per unit volume of flow. The film thickness t in the y-direction is modeled as 0.5 mm, 0.75 mm and 1.0 mm and the plate length as 1.0 m. The concentration of the ice crystals at the inlet condition, has been set at 4 × 10-7, which corresponds to the inlet ice crystal

   diameter (din = 0.025 mm) to initialize the calculations. The density of the ice crystals is assumed to be 921 kg/m3 [28]. The average physical properties of NaCl-H2O solution (salt mass concentration = 2.5%) are assigned constant values as

   = 4 kJ/kgK, = 1.4 × 10-3 kg/ms, = 1024 kg/m3 and kf = 0.6506 W/mK.

4. RESULTS AND DISCUSSION

   The results of the numerical solution are presented in terms of the dimensionless particle diameter D, temperature profiles, local Nusselt number at the wall, particle Nusselt number, and local super-cooling. The flow for different isothermal plate temperatures (Stefan number) and film thicknesses are taken into account. The isotherms and the ice crystal size distribution in the film for the case of = 267 K, = 271 K, Stefan number = 0.05, film thickness of 0.5 mm and = 0.025 mm are shown in Fig. 2.

   As expected, the film temperature decreases in the flow direction and approaches the plate temperature at the end of the computational domain. The reduction of the film temperature is mainly attributed to the fact that during the solidification process, the thermal energy is transferred to the plate surface in the form of latent and sensible heats. During this process, the ice crystals increase in size due to mass transfer from the solution to the surface of the ice particles leaving behind a more concentrated solution that depresses the freezing temperature of the solution. Also, it is found that the crystal diameter increases along the x-direction and near the plate surface (Fig. 2b). This mass transfer process is more significant near the plate surface and increases in the stream-wise direction. However, the mass transfer decreases in the y-direction due to the convection heat gains near the film free surface. Figure 3 shows the variation of ice concentration in the x-direction within the film thickness for Ste = 0.05, t = 0.75 mm and din = 0.025 mm at different y- locations. The ice concentration increases with the increase in the distance along the x-direction. As can be seen, the greatest

   concentrations are near the plate surface. This may be attributed to the fact that the mass transfer process from the solution to the surface of ice particles is more significant near the plate surface and increases in the stream-wise direction which is responsible for increasing the ice particle diameters, and hence, increasing the ice slurry concentration. This behavior of the super-cooling of the mixture is shown in Fig. 4 for different locations within the flow domain.

   The degree of super-cooling of te mixture (Tm T) is maximum in the vicinity of the wall (y = 0) and increases with the flow direction. At the end of the plate, the gradient of the degree of super-cooling almost diminishes.

   The local Nusselt number is defined as:

   = /

   (27)

   Fig. 2: Variation of dimensionless temperature and ice crystal diameter within the flow film for Ste = 0.05, t = 0.5 mm and din = 0.025 mm.

   Fig. 3: Variation of ice concentration within the film thickness at Ste = 0.05, t

   = 0.75 mm and din = 0.025 mm.

   Fig. 4: Supercooling of the fluid within the film thickness y at different x- locations, for Ste = 0.05, t = 0.5 mm and din = 0.025 mm.

   where is the local heat transfer coefficient and it is defined as follows:

   = /( )

   (28)

   The local heat flux at the surface of the plate, qw is obtained by applying Fourier's law to the control volume adjacent to the wall as follows:

   = (,2 )/()

   (29)

   where ,2 is the finite volume temperature adjacent to the plate surface.

   Figure 5 shows the variation of the local Nusselt number in the flow direction for the three cases of the film thicknesses (0.5, 0.75 and 1.0 mm) at Ste = 0.06 mm and din = 0.035 mm. However, the solution temperature decreases along the x- direction and the temperature gradient at the wall and the heat flux increase due to the enhancement in the effective thermal conductivity.

   The increase in the local Nusselt number in the flow direction for t = 0.75 mm and 1.0 mm can be attributed to the increase in the ice particle concentration and the enhancement in the solution's effective thermal conductivity. It has been observed that for the film thickness of 0.5 mm corresponding to Re = 200, the local Nusselt number increases along the plate surface, approaching its maximum value at the middle of the plate, then decreases slightly to an asymptotic value of 590 at the end of the plate. This is due to the enhancement in the effective thermal conductivity of the ice slurry not being equivalent to the increase in the temperature gradient at the wall, consequently decreasing the heat transfer coefficient along the plate surface. Also, the results showed that the average heat transfer coefficient is almost independent of the Stefan number for different film thicknesses, as depicted in Fig. 6.

   The variation of particle Nusselt number within the film thickness for different x-locations is shown in Fig. 7. It can be seen that the particle Nusselt number increases in the flow direction and reaches an asymptotic value at the middle of the film thickness.

   This can be attributed to the continuous increase in the particle diameter, especially in the streamlines beside the wall. The effect of the Stefan number on the distribution of particle

   Nusselt numbers within the film thickness at t = 0.5 mm and din

   = 0.025 mm is depicted in Fig. 8. The increase in Stefan number increases the rate of growth of particle and consequently enhances particle Nusselt number.

   Fig. 5: Nusselt number variation at the wall at Ste = 0.05 and din = 0.025 mm for different film thicknesses.

   Fig. 6: The average heat transfer coefficient versus Ste at din = 0.025 mm for different film thicknesses.

   Fig. 7: Variation of ice crystal Nusselt number within the film thickness at different x-distances for Ste = 0.05, t = 0.5 mm and din = 0.025 mm.

   The effect of film thickness on the particle Nusselt number is shown in Fig. 9. It can be seen that the increase in the film thickness decreases the particle Nusselt number. This can be attributed to the decrease in heat transfer and particle growth rates.

   The growth rate of ice crystals in the falling film at any y- location is the difference between the final and initial crystal diameters divided by the residence time of an ice crystal ( =

   /). The mass of ice growth rate for a plate of width W

   and length L can be obtained from:

   3 = (() ()3) 6 =2

   (30)

   where and are the number of nodes in the z and y- directions, respectively. Which are defined as follows:

   = and =

   (31)

   Fig. 8 Effect of Stefan number on the distribution of particle Nusselt number

   t = 0.5 mm, x = 0.5 and din = 0.025 mm.

   where and are the length of the control volume in the z and y-directions, respectively.

   Figure 10 shows the variation of the ice crystal growth rate with the Stefan number for different film thicknesses. The increase in Stefan number due to the decrease in the wall temperature enhances the heat transfer rate and consequently increases the ice growth rate. Also, it is found that the ice crystal growth rate decreases with film thickness.

5. CONCLUSIONS

A numerical simulation of ice crystal growth in the falling films for the laminar case has been performed. As expected, the results showed that the Nusselt number at the wall appears to be primarily a function of film thickness and is essentially independent of Stefan number. The ice crystal growth rate was seen to be primarily a function of the Stefan number and film thickness. The ice crystal Nusselt number and super-cooling of the solution are dependent on Stefan number, while Reynolds number (film thickness) has a lower influence on these parameters.

Fig. 9: Effect of film thickness on the particle Nusselt number at Ste = 0.05,

x = 0.5 mm and din = 0.025 mm.

Fig. 10 Variation of ice crystal growth rate with Stefan number for different film thicknesses and din = 0.025 mm.

REFERENCES

1. Kawabe, H. Fukusako, S. Yamada, M. Yanagita, K., 1998. Continuous production characteristics of slush ice by use of a horizontal oscillating cooled wall. Trans. JSRAE. 15 3, 193-201 (in Japanese)

2. Hirata, T., Nagasaka, K., Ishikawa, M., 2000. Crystal ice formation of solution and its removal phenomena at cooled horizontal solid surface (part I: ice removal phenomena). Int. J. Heat and Mass Transfer. 43 (3), 333-339.

3. Hirata, T. Kato, M. Nagasaka, K. Ishikawa, M., 2000. Crystal ice formation of solution and its removal phenomena at cooled horizontal solid surface (part II: ice removal condition). Int. J. Heat and Mass Transfer. 43 (5), 757-765

4. Frei, B., Egolf, P.W. 2000. Viscometry applied to the Bingham substance ice slurry. Proceedings of the Second Workshop on Ice Slurries of the International Institute of Refrigeration. May 2526 2000, Paris, France, 4860.

5. Stewart, W., Kaupang, R., Tharp, C., Wendland, R., Stickler, L., 1993. An Approximate Numerical Model of Falling-Film Ice Crystal Growth for Cool Thermal Energy Storage. ASHRAE Trans. 1993, vol. 99, 347- 355.

6. Pronk, P., Hansen, T.M., Ferreira, C.A., Witkamp, G.J., 2005. Time- dependent behavior of different ice slurries during storage. Int. J. of Refrigeration. 28, 2736.

7. Kozawa, Y., Aizawa, N., Tanino, M., 2005, Study on ice storing characteristic in dynamic-type ice storage system by using supercooled water. Effects of the supplying conditions of ice-slurry at deployment to district heating and cooling system, Int. J. Refrig., vol. 28: p. 73-82.

8. Cliche, A., Lacroix, M., 2006. Optimization of ice making in laminar falling films. Energy Conversion and Management. 47, 22602270.

9. Matsumoto, K., J. L. Sarmiento, and M. A. Brzezinski 2002. Silicic acid leakage from the Southern Ocean: A possible explanation for glacial atmospheric pCO2, Global Biogeochem. Cycles, 16(3), 1031, doi:10.1029/2001GB001442.

10. Matsumoto, K. Sakurai, H., 2006. Study on revention of ice adhesion to cooling wall due to voltage impression in ice storage- discussion on possibility of attraction of oil to wall. Int. J. of Refrigeration. 29, 142 149.

11. Matsumoto, K., Oikawab, K., Okada, M., Teraoka, Y., Kawagoe, T., 2006. Study on high performance ice slurry formed by cooling emulsion in ice storage (discussion on adaptability of emulsion to thermal storage material). Int. J. of Refrigeration. 29, 1010-1019.

12. Matsumoto, K., Suzuk, Y., Okada, M., Teraoka, Y., Kawagoe, T., 2006 Study on continuous ice slurry formation using functional fluid for ice storage- Discussion of optimal operating conditions. Int. J. of Refrigeration. 29, 1208-1217

13. Matsumoto, K., Sakae, K., Oikawac, K., Okada, M., Teraoka, Y., Kawagoe, T., 2007. Study on formation of high performance ice slurry by W/O emulsion in ice storage (discussion on characteristics of propagation of supercooling dissolution). Int. J. of Refrigeration. 30, 1300-1308

14. Matsumoto, K., Sakae, K., Yamauchi, H., Teraoka, Y., 2008. Formation of high performance ice slurry by W/O emulsion in ice storage (effective method to propagate super-cooling dissolution). Int. J. of Refrigeration. 31, 832-840.

15. Matsumoto, K., Okada, M., Kawagoe, T., Kan, C., 2000. Ice storage system with water-oil mixture (formation of suspension with high IPF). Int. J. of Refrigeration. 23 (5), 336-344.

16. Swanson B. and Nelson, J., 2019. Low-temperature triple-capillary cryostat for ice crystal growth studies. Atmos. Meas. Tech., 12, 6143 6152.

17. Rutgers., I.R.,1962. Relative Viscosity of Suspensions of Rigid Spheres in Newtonian Liquids. Rheologica Acta. 24, 305348.

18. Levich, V., Physiochemical hydrodynamics. New Jersey, USA, 1962.

19. Vand, V., 1948. Viscosity of Solutions and Suspension. J. Phys. Coll. Chem. 52, 300321.

20. Stewart, W., Stickler, L., 1991. A Pumpable Ice Slurry for Cool Thermal Storage. Proceedings of the 5th International Conf. on Thermal Energy Storage. Scheveningen, Netherlands, May 13-16, (6-1)- (6-5).

21. Maxwell, J.C., A treatise on Electricity and Magnetism. New York: Dover Publications, USA, 1954.

22. Charunyakorn, P., Roy, S.K., 1991. Forced Convection Heat Transfer in Microencapsulated Phase-Change Material Slurries Flow in Circular Ducts. Int. J. Heat and Mass Transfer. 34, 819832.

23. Sohn, C. and Chen, M., 1984. Microconvective Thermal Conductivity in Disperse Two-Phase Mixtures as Observed in a Low Velocity Couette Flow. Experiment. J. Heat Transfer. 106, 539542.

24. Stewart, W., Kaupang, R., Tharp, C., Wendland, R., Stickler, L., 1993. An Approximate Numerical Model of Falling-Film Ice Crystal Growth for Cool Thermal Energy Storage. ASHRAE Trans. 1993, vol. 99, 347- 355.

25. Burns, A.S., Stickler, L.A, Stewart. W.E., 1992. Solidification of an Aqueous Salt Solution in a Circular Cylinder. ASME Trans. J. Heat Transfer. 114, 3033.

26. Fang, L.J., Cheung, F.B., Linehan J.H., Pedersen D.R., 1984. Selective Freezing of a Dilute Salt Solution on a Cold Ice Surface. J. of Heat Transfer. 106, 385393.

27. Patankar, S.V., Numerical heat transfer and fluid flow. Hemisphere Publishing Co, New York, USA, 1980.

28. Cengel, A., Boles, M.A., Thermodynamics: An Engineering Approach.

McGraw Hill, 2006.