Strojniški vestnik - Journal of Mechanical Engineering 51(2005)7-8, 456-461 UDK-UDC 536.2 Izvirni znanstveni članek - Original scientific paper (1.01) A Numerical Study of Thermosolutal Melting Thomas J Scanlon1 and Matthew T Stickland1 1Department of Mechanical Engineering, University of Strathclyde, Glasgow, G1 1XJ, Scotland * e-mail tscanlon@mecheng.strath.ac.uk Abstract This paper describes the numerical investigation into the melting of a pure ice block into an aqueous solution of sodium carbonate (Na2CO3). The numerical study is concerned with capturing the evolving solid-fluid interface during phase change while solving sequentially the double-diffusive conditions resulting from the combined effects of thermal and solutal buoyancy in the flow field. The results show that, with a relatively simple phase change model incorporated into a fixed grid, finite volume numerical formulation, reasonable concurrence may be obtained in comparison with published experimental data. Introduction The melting of ice in mixtures began as a problem for the melting of glaciers in sea water motivated by oceanographers. The melting of glaciers in sea water has assumed further importance recently because of the apparent increasing temperature in the earth’s environment. Therefore, it is important to analyse correctly the characteristics of the melting of ice in mixtures. Such natural convection phenomena, driven by the twin competing buoyancy forces due to concentration and temperature gradients, have been of interest to many authors during the latest decades. The thermosolutal natural convection is seen to play an important role in controlling the characteristics of solid-liquid phase change in mixtures [1-2]. It has been found that the time evolution of the solid-liquid interface and the local equilibrium conditions existing at the interface are strongly coupled to the interface heat and mass transfer due to thermosolutal convection. In solidification processes, changes in the time evolution of this coupling is one of the principal criteria which governs the front velocity and any micro or macro segregation effects. Such effects are thus considered to impinge directly on the homogeneity of the forming solid [3-8]. A significant number of experimental and numerical studies have been devoted to double-diffusive convection associated with melting or solidification. In the case of solidification, the existence of a mushy zone adjacent to the solid-liquid interface is a specific problem that will not be considered in this paper. In the context of melting, a series of fundamental experimental and theoretical studies have been executed concerning buoyancy-induced flows driven by combined thermal and solute transport near a vertical melting ice surface in saline water [9-12]. Considering melting in a cavity [13-21], an analytical and experimental study has been carried out [13] which focused upon the different modes of double-diffusive convection arising from a horizontal melting interface. Using the same configuration, a quantitative model and scaling analyses of melting or dissolving driven by compositional convection has also been developed [14-16]. When the phase change interface is vertical, multicellular flow structures are seen to arise in the cavity that strongly influence the melting characteristics [17-21]. Numerical Methodology The aim of the computational is to develop, as far as possible, a reliable model of the thermo-fluid behaviour of a block of pure ice melting into an aqueous solution of sodium carbonate under conditions of natural convection as detailed in figure 1. The computational model would be based on a coupled thermo-fluid system with laminar natural convection heat and mass transfer with phase change occurring. The aim would be to establish the critical modelling parameters encountered in such flows. The computational model is undertaken within the framework existing commercial computational fluid dynamics analysis software FLUENT which 456 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 456-461 Nomenclature gßT ?TH 3 RaT = av = thermal Rayleigh number ? = thermal diffusivity of the fluid (m2/s) ? = kinematic viscosity of the fluid (m2/s) SfajrtintSMO tet Thertnasollitit Cells ^ Spill' Al | ."Sirij,- .','.JIU Thermal C*M I J ." . i. j:v-v-T'.: I ¦ ¦¦: rfa*iflfYuM OHIHI Figure 1 Tyical flow structures involved in thermosolutal melting. offers current capabilities for coupled thermal and fluid analysis but cannot address the details of the thermo-physical modelling being required here. Such features need to be incorporated using user subroutines. Such work has recently been described by Scanlon and Stickland [22-27]. This numerical methodology proposes that the governing equations are solved in a manner such that if the temperature falls below the freezing isotherm then the convection terms in the equations of motion are effectively disengaged. Variations in the specific heat of the material are incorporated in order to account for the phase change according to the graph shown in figure 2, where ?hfus is the latent heat of fusion of ice, ?T is the finite temperature bandwidth over which phase change is assumed to occur and PCZ is an acronym for the phase change zone. Cp* k P C Solid Z Liquid ?T «- T The model of Scanlon and Stickland [22-27] has been shown to be a stable, computationally efficient, uncomplicated numerical technique that reasonable solutions to complex problems natural convection with phase change, as figures 3 and 4. provides involving shown in I C E i* *-¦¦¦ .- ^ .5C* - — |B Figure 3 Experimental results (particle tracers) for the freezing of pure water on a vertical wall [24]. Figure 2 Variation in specific heat with temperature Figure 4 Numerical prediction (streamlines) for the freezing of pure water o a vertical wall [24]. A numerical study of thermosolutal melting 457 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 456-461 The experimental work of Bénard et al [20], to which the numerical solution will be compared, is detailed in figure 5 Table 1 Parameters considered in the numerical study Hm Umitim,1« ¦i P 211 "TI ICbUS) L = 0.07 m Figure 5 Experimental set up of Bénard et al [20]. The experimental cell consists of a rectangular cavity (height H = 20 cm, width L = 7 cm and depth P = 20 cm) designed to generate a two-dimensional flow. The two opposing walls are heat exchangers which are differentially heated at constant and uniform temperatures controlled by two independent thermal loops where a water-glycol mixture is circulated. For the experimental conditions considered in this paper, pure water (concentration C = 0) is contained within distance e. the rest of the cell contains a still solution of Na2CO3 of concentration Ci = 2 % by mass. This aqueous solution is separated from the pure water by a thin metal plate. Initially, the whole experimental system is maintained a few tenths of a degree under 0 oC after which a rectangular block of pure ice will have formed and is attached to the cold wall. The initial dimensions of the ice block are 3 cm thickness (distance e) and height Hi = 14 cm. This implies an initial length for the aqueous solute of 4 cm. For an initial solute concentration of 2 % (by mass), the corresponding equilibrium liquidus temperature is –0.66 oC, so the solute mixture remains in its liquid phase. At time t = 0 the metal barrier is removed very slowly such that the fluid is not disturbed and the hot exchanger is raised quickly to a temperature of Th = 10 oC which is maintained throughout the experiment. The cold exchanger temperature is kept at Tc = 0 oC. For the numerical solution the bottom wall is considered to be adiabatic while at the top surface symmetric solutions are assumed. Thermosolutal convection is characterized by five dimensionless parameters; the Prandtl number Pr, the Lewis number Le, the thermal Rayleigh number RaT, the buoyancy ratio N and the aspect ratio A (A = Hi/e in Figure 5). This leads to the conditions outlined in Table 1 which are considered in this study. Such a range of parameters may be varied through the initial concentration of the liquid and the temperature difference Pr 11 Le 190 RaT 3.8 x 108 N -8 A 4.667 between the heat exchangers. The buoyancy ratio N represents the effects of solutal buoyancy to thermal buoyancy where N is defined as: N = ßT?T (1) For an initial concentration Ci = 2 % by mass, a temperature difference ?T = 10 K and a thermal expansion coefficient ßT = 3.275 x 10-4 K-1 this produces a solutal expansion coefficient ßS = -1.31002 wt percent-1. The experimental conditions are such that the initial solute mass fraction ranges between 2 % and 5 %. Under such conditions, the density is observed to be a linear function of temperature and concentration such that the standard linear Boussinesq approximation can be applied and the density inversion effects found in pure water are negated. In this manner the density may be considered as a linear function of temperature and concentration according to [20]: ? = ?0(1-ßT?T-ß SC) (2) where ?0 = 999.8 kg/m3 is the density of pure water at T = 0.01 oC. In order to incorporate mass transfer an additional transport equation is required for the solute transport with appropriate initial and boundary conditions. For the initial solute conditions a concentration is applied to the ice (C = 0) and solute (Ci = 0.02) regions respectively. As the numerical solution progresses, it is only when melting of an initially solid cell in the ice zone occurs that the nascent fluid cell will allow the solute to be advected into the melt flow. For the numerical solution a computational mesh of 180000 rectangular cells was employed and a time step of 4 s was adopted. For the pressure-velocity coupling the conventional SIMPLE scheme was used. Temporal discretisation was based on an implicit method and for convective discretisation the QUICK scheme was employed. This led to total computational times of approximately 30 hours on a Pentium4 3GHz PC with 512 Mb RAM. Results and Discussion The following results show the comparison between the numerical work (contours of solute concentration) and the experimental results of Bénard et al [20] (observations by fluorescence – lighter regions = lower concentration): 458 Scanlon T.J. - Stickland M.T. Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 456-461 Fig. 6 Experiment (left), Numerical (right), t = 7 min Fig. 9 Experiment (left), Numerical (right), t = 25 min Fig. 7 Experiment (left), Numerical (right), t = 15 min Fig. 10 Experiment (left), Numerical (right), t = 30 min Fig. 8 Experiment (left), Numerical (right), t = 20 min Fig. 11 Experiment (left), Numerical (right), t = 35 min A numerical study of thermosolutal melting 459 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 456-461 In figures 6-11 we can observe the temporal evolution of the thermosolutal melting process over a period of 35 minutes. It is evident that the competing effects of thermal and solutal buoyancy in the melt flow result in the non-linear solid-fluid interface profiles observed. The general form of this interface is of a progressive ‘waisting’ effect around the central core of ice with thicker upper and lower solid regions. From the numerical study it was observed that horizontal fingers of warmer fluid existed which traversed the fluid region and effectively ‘bit’ into the ice via the melting process. Such multicellular regions of varying concentration appear to have been captured qualitatively by the numerical techniques employed in this study. Conclusions A numerical study of thermosolutal melting has been carried out. The numerical techniques employed have attempted to capture the time-evolving interface of pure ice melting into an aqueous solution of Sodium Carbonate and the resulting double-diffusive flow with its competing influences due to thermal and solutal buoyancy. The numerical results show a reasonable qualitative concurrence with published experimental work with the general trends of solid-fluid interface evolution and thermosolutal melt flow characteristics being captured. References [1] Huppert, H. E. and Turner, J. S., Ice Blocks Melting into a Salinity Gradient, Journal of Fluid Mechanics, Vol. 100, pp. 367-384, 1980. [2] Huppert, H. E., The Fluid Mechanics of Solidification, Journal of Fluid Mechanics, Vol. 280, pp. 287-302, 1990. [3] Prescott, J. P. and Incropera, F. P., Convective Transport Phenomena and Macrosegregation During Solidification of a Binary Metal Alloy: I-Numerical Predictions, ASME J. Heat Transfer, 116, pp. 735-741, 1994. [4] Prescott, J. P. and Incropera, F. P., Convective Transport Phenomena and Macrosegregation During Solidification of a Binary Metal Alloy: II- Experiments and Comparisons with Numerical Predictions, ASME J. Heat Transfer, 116, pp. 742-749, 1994. [5] Campbell, T. A. and Koster, J. N., Visualizations of Liquid-Solid Interface Morphologies in Gallium Subject to Natural Convection, J. Cryst. Growth, 140, pp. 414-425, 1994. [6] Yin, H. and Koster, J. N., Double Diffusive Convective Flow and Interface Morphology During Transient Ga-5 percent in Alloy Melting, J. Cryst. Growth, 217, pp. 170-182, 2000. [7] Tanny, J., Experimental Study on the Crystallization of a Binary Melt at the Vertical Boundary of an Enclosure, Int. J. Heat Mass Transf., Vol. 38, pp. 1141-1150, 1995. [8] Wettlaufer, J. S. and Worster, M. G., Natural Convection During Solidification of an Alloy from Above with Application to the Evolution of Sea Ice, Journal of Fluid Mechanics, Vol. 111, pp. 291 – 316, 1997. [9] Joseberger, E. G. and Martin, S., A Laboratory and Theoretical Study of the Boundary Layer Adjacent to a Vertical Melting Ice Wall in Water, Journal of Fluid Mechanics, Vol. 111, pp. 439-473, 1981. [10] Carey, V. P. and Gebhart, B., Transport Near a Vertical Ice Surface Melting in Saline Water: Some Numerical Calculations, Journal of Fluid Mechanics, Vol. 117, pp. 379 – 402, 1982. [11] Carey, V. P. and Gebhart, B., Transport Near a Vertical Ice Surface Melting in Saline Water: Experiments and Low Salinities, Journal of Fluid Mechanics, Vol. 117, pp. 403 – 423, 1982. [12] Sammakia, B. and Gebhart, B., Transport Near a Vertical Ice Surface Melting in Water of Various Salinity Levels, Int. J. Heat Mass Transf., Vol. 26, pp. 1439-1452, 1983. [13] Woods, A. W., Fluid Mixing During Melting, Phys. Fluids A, 3, pp. 1393-1404, 1991. [14] Woods, A. W., Melting and Dissolving, Journal of Fluid Mechanics, 239, pp. 429-448, 1992. [15] Kerr, R. C., Melting Driven by Vigorous Compositional Convection, Journal of Fluid Mechanics, 280, pp. 255-285, 1994. [16] Kerr, R. C., Dissolving Driven by Vigorous Compositional Convection, Journal of Fluid Mechanics, 280, pp. 287-302, 1994 [17] Sugawara, M. and Irvine, T. F., The Effect of Concentration Gradient on The Melting of a Horizontal Ice Plate From Above, Int. J. Heat Mass Transf., Vol. 43, pp. 1591-1601, 2000. [18] Beckerman, C. and Viskanta, R., Double-Diffusive Convection Due to Melting, Int. J. Heat Mass Transf., Vol. 31, pp. 2077-2089, 1988. [19] Beckerman, C. and Viskanta, R., An Experimental Study of Melting of Binary Mixtures With Double-Diffusive Convection in the Liquid, Exp. Therm. Fluid Sci., Vol. 2, pp. 17-26, 1989. [20] Benard, C., Benard, R., Bennacer, R. and Gobin, D., Melting Driven Thermohaline Convection, Physics of Fluids, Vol. 8, pp. 112-130, 1996. [21] Mergui, S., Geoffroy, S. and Benard, C., Ice Block Melting Into a Binary Solution: Coupling of the Interfacial Equilibrium and the Flow Structures, ASME Journal of Heat Transfer, Vol. 124, pp. 1147 – 1157, 2002. [22] Scanlon, T.J., Stickland, M.T. and Robertson, N.A., A Numerical and Experimental Analysis of Natural Convection Ice Melting. Presented at 460 Scanlon T.J. - Stickland M.T. Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 456-461 the 6th UK National Conference on Heat Transfer, Heriot-Watt University, 15-16 September 1999. Published in IMechE Conference Transactions 1999-7-C565/064, pp. 443-452, 1999. [23] Scanlon, T.J., and Stickland, M.T., An Experimental and Numerical Investigation of Natural Convection Melting. Int. Comm. Heat Mass Transfer, Vol. 28, No. 2 (01), pp. 181-190, 2001. [24] Scanlon, T.J. and Stickland, M.T., A Numerical Analysis of Buoyancy-Driven Melting and Freezing , Int. J. Heat Mass Transf., 47, pp.429436, 2004. [25] Stickland, M.T., Scanlon, T.J., Oldroyd, A. and Waddell, P. An Experimental and Computational Analysis of Buoyancy Driven Flows by LaserSheet Tomography, Particle Image Velocimetry and Computational Fluid Dynamics. Proceedings of the 1st Pacific Symposium on Flow Visualisation and Image Processing, Honolulu, USA, February 23-26, Vol. 1, pp.203-208, 1997. [26] Scanlon, T.J. A Numerical Model of the Ice Melting Process. Presented at the 3rd International Conference on Advances in Heat Transfer, Montreal, Canada, May 2000. Published in Advances in Fluid Mechanics III, Ed. M. Rahman and C. A. Brebbia, pp. 313-322, 2000. [27] Scanlon, T.J., Stickland, M.T. and Lacombe, M., A PIV Analysis of Natural Convection Ice Melting, European Association of Laser Anemometry Conference, September 2001, Limerick, Eire. A numerical study of thermosolutal melting 461