3D numerical and experimental study of gallium melting in a rectangular container

International Journal of Heat and Mass Transfer 67 (2013) 260–271

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

3D numerical and experimental study of gallium melting in a rectangular container O. Ben-David, A. Levy ⇑, B. Mikhailovich, A. Azulay Pearlstone Center for Aeronautical Engineering Studies, Department of Mechanical Engineering, Ben-Gurion University of the Negev, Israel

a r t i c l e

i n f o

Article history: Received 18 December 2012 Received in revised form 6 July 2013 Accepted 19 July 2013 Available online 4 September 2013 Keywords: Gallium melting Convective flow Liquid–solid interface Numerical model 3D effects Ultrasonic velocimetry

a b s t r a c t Gallium melting in a rectangular container with a heated side face has been investigated. The focus of the research is the advancement of the numerical model to 3D status taking into consideration thermal variations of medium properties and the presence of mushy zone, and the model experimental verification using known data and new data obtained on a specially developed experimental setup. This setup is oriented to trace the liquid–solid interface position and profile and melt velocity using ultrasonic Doppler velocity measurements in the liquid phase. The numerical model was based on COMSOL Multiphysics software. 2D and 3D versions were built to calculate the melting regime and the flow process characteristics. Outputs of both versions were compared with known experimental data and new data obtained in the present study. The results revealed a better agreement between 3D computational and experimental data indicating to a profound effect of the container boundaries on the problem under study. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Metal melting in a rectangular container has been studied for many years [1–18], but the respective analytical solutions are known only for the simplest cases. Meanwhile, heat convection and metal properties dependence on temperature play an important role in real melting processes. They complicate the analysis and make it necessary to perform computer simulations and validation by low-temperature physical experiments. Both kinds of investigations complement each other, since the assessment of computer simulation efficiency is realized by comparison with available experimental data, and in a number of cases a certain discrepancy between computed and experimental results is observed [3,18,20]. The latter can be due to either experiment imperfections or computation methods deficiency. The known numerical methods make it possible to solve melting and solidification problems taking into account both the heat conduction mechanism, which is essential at the initial melting stage, and the convective heat transfer playing an important role in the dynamics of the melt and liquid–solid interface [17,18]. In this research, computer simulations and physical low-temperature experiments have been carried out for gallium melt in a rectangular container with a heated side face (Fig. 1). To improve the comprehension of such a process, it is necessary to estimate ⇑ Corresponding author. Address: P.O. Box 653, Beer-Sheva 84105, Israel. Tel.: +972 8 6477092; fax: +972 8 6472814. E-mail address: [email protected] (A. Levy). 0017-9310/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.07.058

the influence of the numerical model dimension (2D and 3D), temperature dependence of the metal physical properties and boundary conditions on the problem solution. Therefore, it is advisable to specify the necessary aspects and to modify computational approach and adapt experimental methods as required. Industrial-scale hot melting experiments are very expensive and complicated. Moreover, as known, until now there are no workable technical means for velocity measurements in high-temperature liquid metals. Therefore, melting experiments at low temperature using gallium or PCM (phase change material) are performed in various forms in order to understand the melting process and heat transfer and flow mechanisms. Gau and Viskanta [1,2], used pure gallium and performed vertical and horizontal melting, measuring the solid–liquid interface temperature and position. For this purpose, thermocouples and special L-shaped probes for draining the container filled with melt were used to expose the solid–liquid interface. These measuring methods are invasive, disturb the melting process and do not provide information about the melt flow. A noninvasive method was applied by Campbel and Koster [3] who used gallium and measured the solid–liquid interface shape and position by real-time radioscopy method. X-ray observation showed a sharp contrast between molten and solid phases with a 3% density difference. Another method can be found in Menon et al. [4], Pal and Joshi [5] and Katsman et al. [6]. All of them used different kinds of PCM and performed melting experiments in glass tubes, which allowed the registration of the solid–liquid interface shape and position by a high-contrast camera. These methods are

O. Ben-David et al. / International Journal of Heat and Mass Transfer 67 (2013) 260–271

261

Nomenclature B cp C Fl ~ F H DH L f(T) k T P t g d V ~ u x,y,z Q I h Ra

computational constant (kg/ms) specific heat (J/kg °C) Carman–Kozeny constant (kg/m3 s) liquid–solid volume fraction force vector (N) enthalpy (J/kg) modified latent heat (J/kg) latent heat (J/kg) enthalpy temperature function (J/kg) thermal conductivity (W/m °C) temperature (°C or K) pressure (Pa) time (s) gravity acceleration (m/s2) diameter (m) volume (m3) velocity vector (m/s) coordinate (m) heat quantity (J) unit matrix container height Rayleigh number, (q  qm)gh3/qmalml

indeed noninvasive, but they cannot provide any quantitative information about the melt flow in the liquid phase. Moreover, they are not suitable for opaque media such as metals. All experimental results show a great influence of convective flow on the melting process and the solid–liquid interface shape. One of noninvasive measurement methods in phase-change experiments is the use of ultrasound instruments. Beginning with medical applications, this method of measuring the velocity of fluid flows found its use in scientific and engineering research, especially when dealing with opaque and aggressive liquids, including liquid metals [7] and [8]. As an example, one can cite the study [9], where mean velocities of the melt flow and liquid–solid interface have been found experimentally during the study of the magnetic field effect on the alloy solidification. Unlike known experiments using UDV (ultrasonic Doppler velocimeter) technique, in our study we define the liquid–solid interface position and profile, its displacement and longitudinal mean velocity during the metal melting process. There are two main approaches to the calculation of moving solid–liquid interface. The first can be found in Ramachandran et al. [10], Gadgil and Gobin [11], Albert and O’Neill [12] and Vynnycky and Kimura [13] who have developed deformable grids which

Fig. 1. Melting process scheme.

Pr Prandtl number ml/al Ste Stefan number cpl(TH  TC)/L TH TC Tint Temperatures of hot and cold walls and initial temperature, accordingly Greek letters l dynamic viscosity (kg/ms) D difference q density (kg/m3) e computational constant m kinematic viscosity (m2/s) K aspect ratio a thermal diffusivity (m2/s) Subscripts m melting s solid l liquid p particle

follow the moving solid–liquid interface. The moving boundary approach may be appropriate when tracking a sharp boundary interface is required. However, specific models describing the interface motion are needed. Since these models use several simplifying assumptions, the accuracy of interface movement and shape is questionable. Following [13], the analytical solution describing the moving solid–liquid boundary ‘‘is not uniformly valid for all time, and tends to overestimate the final thickness of the solid layer’’. Besides, the moving boundary approach does not consider the mushy zone and is more complicated numerically, and therefore the calculation time for the 3D problem can be much longer. An alternative approach is the use of a fixed grid, which simplifies the numerical modeling requirements. Examples of such approach can be found in Morgan [14], Gartling [15] and Voller et al. [16–18]. In this approach, conservation equations are solved for the entire volume, and therefore the velocity of the solid phase should be zero. Gartling [15] employs a simple way of making the viscosity of the metal a function of enthalpy: by decreasing extra enthalpy DH from the latent heat value to zero, he drives the viscosity to a high value. 2. Present study In the present study, the enthalpy-porosity method, which allows a fixed-grid solution of the momentum and energy equations, was adopted. This method was developed by Voller and Prakash [17]. It is based on solid–liquid volume fraction, while the phase change is assumed to occur within a narrow temperature range. In addition, the Darcy source approach is used to simulate motion in the mushy zone, which consists of solid particles dispersed in a liquid. Because of lack of information about the solid particles diameter, the Carman–Kozeny equation [19] for modeling drag forces in porous media is used. The momentum equation contains a source term, which depends on the local liquid fraction Fl and the Carman–Kozeny constant C. When the local liquid fraction equals zero, there is only solid at this location. The Carman–Kozeny constant is mostly influenced by the morphology and viscosity of the mushy zone. Voller and Prakash [17] and Shmueli et al. [20] investigated the influence of this constant on the solution using a

262

O. Ben-David et al. / International Journal of Heat and Mass Transfer 67 (2013) 260–271

number of constant values. It was found that there exist more suitable values for a specific experiment. It was noted that the constant C should be investigated for different geometries, materials and boundary conditions. In our research, the constant C was also examined and used as a gauge constant in order to calibrate the numerical model and the experiments results. In our research, the energy equation is expressed by enthalpy formulation. Voller and Prakash [17] presented the enthalpy of the material as a sum of sensible heat and latent heat, H = cpT + DH, where the latter term is a function of temperature DH = f(T). This function allows estimating a mushy zone region from the phasechange equilibrium diagram of the material. In the present research, the enthalpy-porosity method was adopted, as noted above, but unlike the studies mentioned, the enthalpy-porosity method was used both in 2D and 3D calculations for taking into account boundaries and their effect on the melting process. Moreover, in the present research, temperature dependence of metal physical properties was used unlike other numerical studies, which assume a small change and neglect temperature dependence. At the same time, some details of the melting process remain unclear, such as the number of rolls in the liquid metal flow [21], particularly at the early stage of melting. Various authors obtained different melt flow structures in their computations [18,22,23]. Moreover, scanty particulars are known about 3D simulations of this problem [24], and experimental data about liquid phase flow have not practically been described (principles of liquid metal flow measurements are well-known and employed in alloys solidification process [9,25], but they are not used in melting problems). 3. Mathematical model and numerical analysis 3.1. Theoretical model and setup A model is built for describing gallium melting in a rectangular container (see Fig. 1). The initial temperature TInt < Tmelt, and all the metal is in the solid state. The process starts when the left vertical side wall temperature abruptly becomes higher than the melting temperature, i.e. THot > Tmelt. When considering a phase change in an alloy (or non-pure metal), three regions are present: a solid region, a liquid region and a porous media region which consists of solid particles dispersed in the liquid (see Fig. 1). A local liquid–solid volume fraction can be written for a pure material melting at a specific temperature:

F l ðTÞ ¼

Vl ¼ V



0 T < Tm; 1 T > Tm:

ð1Þ

For alloys melting, which occurs over a certain temperature range, the liquid volume fraction can be calculated by [17]:

8 > 0 V l < TT s F l ðTÞ ¼ ¼ T l T s > V : 1

T < T s; Ts 6 T 6 T l;

ð2Þ

T > T l:

When the temperature is lower than the solidus temperature, the entire cell is solid. At the temperature higher than the liquidus temperature, the cell is liquid, and at the temperature between the solidus and liquidus it is porous. Tl  Ts was defined as a small interval around the melting temperature. Beckermann and Viskanta [26] noted that an unrealistically large Tl  Ts (i.e.>2 °C) would introduce too much artificial ‘smearing’ of the predicted interface, while a too small interval will complicate numerical calculation and require a very fine mesh. Beckermann and Viskanta [26] used the value Tl  Ts/TH  TC = 0.02,

which means in this case that Tl  Ts is within the range of 0.3– 0.5 °C. In the present research, Tl  Ts = 0.5 °C was chosen. The mathematical model of metal melting by a heated vertical wall was examined, for example, by Webb and Viskanta [27]. Here we describe the problem considering mass, momentum and energy balance equations under the assumption of a compressible laminar flow of a Newtonian liquid. In this study we adopted a simplifying assumption of Prakash and Voller [17] with some modifications: continuity equation:

@q þ r  ðq~ uÞ ¼ 0; @t

ð3Þ

momentum equation:

q

    2 @~ u u  rÞ~ u ¼ rP þ r  l r~ u þ r~ uT  lðr  ~ uÞI þ qð~ @t 3 ~ ~ þ F buoyancy þ F drag ;

ð4Þ

energy equation:

q

@H þ qr  ðH~ uÞ ¼ r  ðkrT Þ: @t

ð5Þ

F buoyancy term in the momentum equation is the force stimulatThe ~ ing natural convection in the container, which is described by the temperature dependence of density:

~ F buoyancy ¼ qðTÞ~ g:

ð6Þ

The flow is calculated as compressible in order to take into account small changes of density due to temperature differences. Here, owing to the use of temperature dependence of gallium physical properties, the expansion coefficient is not necessary. The term ~ F drag in the momentum equation describes drag forces between the liquid and solid phases in the mushy zone. Therefore, it becomes zero outside the mushy domain. The pressure drop, which is related to the drag force between the liquid and solid phases in a fluidized bed, can be calculated by the well-known empirical Ergun expression [28]:

rP ¼

150ð1  F l Þ2 ll 2

F 3l dp

! 1:75ql ð1  F l Þ þ juj ~ u: dp

ð7Þ

The first term in the right-hand part of the last equation represents the Darcy term, and the second is the Forchheimer term. Here ~ u is the superficial velocity of the liquid metal, and dp is solid particles diameter in the mushy zone. Due to the low fluid velocity in the mushy zone, the Forchheimer term can be neglected [26], and the Darcy term can be replaced with Carman–Kozeny equation [20], which deals with drag forces in porous media: 2

Cð1  F l Þ ~ ~ F drag ¼  u: e þ F 3l

ð8Þ

Here e is a computational constant used to avoid the division by zero when the fluid volume fraction becomes vanishingly small and tends to zero [20]. When the cell is completely liquid Fl = 1, the source acquires a zero value, and the momentum equation becomes suitable for an actual one-phase fluid flow. When Fl decreases in the mushy zone, the value of the source increases, and the momentum equation corresponds to the Darcy law. When the local liquid fraction equals zero, the source term becomes very large, and the velocity of the liquid diminishes to zero. In such an approach, the velocity of the solid should be equal to zero. C is a constant which depends on the morphology of the mushy zone. Voller and Parkash [17] noted that the constant C can reduce or increase the flow intensity at the mushy zone. This aspect will be discussed later in the conclusions.

O. Ben-David et al. / International Journal of Heat and Mass Transfer 67 (2013) 260–271

263

2D model was 90,000. For the 3D model, two grids and discretization were examined:

Table 1 The investigated cases. Case

DT(°C)

Thot(°C)

Tint(°C)

1 2 3

15 18 25

40 45 50

25 27 25

(a) 32,144 elements with discretization of the first order for the pressure and heat transfer equations and of the second order for the flow equation with 176,610 degrees of freedom, (b) 10,469 elements with discretization of the first order for all the conservation equations with 60,000 degrees of freedom.

In order to specify zero velocity in solid cells, Gartling’s method [15] has been used. The dynamic viscosity in solid cells should acquire high values in comparison with other terms in the momentum equation, which is achieved following this method [14] by using a large computational constant B. The values of the constant B were examined in the range from 10 to 10,000, and B = 1000 kg/ ms was found to be the lowest order of magnitude which diminishes the velocity magnitude correctly in the entire solid region.

l ¼ ll þ B½1  F l ðTÞ:

ð9Þ

The enthalpy of gallium was calculated by:

H ¼ cp T þ DH;

ð10Þ

where cp is the average specific heat of the solid and liquid metal according to the volume fraction value Fl, and DH is the modified latent heat added to the enthalpy in cells with melting metal, which is expressed by a normal Gaussian formula:

 2 L ðT  T m Þ DH ¼ pffiffiffiffi exp  : ðT l  T s Þ pðT l  T s Þ

ð11Þ

Here L is the latent heat of gallium, Tm is the mean temperature value, and Tl  Ts is temperature variance. The boundary condition for the momentum equation is a non-slip condition at the container walls. At the upper boundary, a pressure outlet to the atmosphere is taken into account. In the energy equation, THot and TInt correspond to the two vertical walls, and the remaining walls are insulated.

The comparison between the results was examined (see Fig. 2), and only minor differences can be seen at the solid–liquid interface profile and liquid fraction volume. Therefore, method (b) was chosen for decreasing calculation time and computer resources. For method (b) with Tl  Ts = 0.5 °C, energy balance was calculated in order to verify the numerical scheme. Fig. 3 shows a comparison between the heat required for heating and melting the gallium, Qmelt, and the heat difference between the inlet and outlet, Qin  Qout. The agreement is very good, the mean error is about 3%, and the energy balance in the numerical scheme is observed. Time discretization was computed by backward differentiation formula (BDF) algorithm [30] which estimates the next time step size for the convergence. BDF algorithm always aims at an increase in the next time step size for shortening the calculation time without damaging the convergence of the solution. Segregation was computed by the conservation equations, the flow equations (~ u; P) were solved by PARDISO [31] solver (parallel sparse direct linear solver), and the energy equation (T) was solved by MUMPS [32] solver (multifrontal massively sparse direct solver). The iterative calculation was performed with an under-relaxation factor allowing the convergence in high-gradient problems. For smaller under-relaxation factors, the convergence becomes slower, and the number of iterations grows. After examining a number of under-relaxation factors, optimal values were found. For the flow equation, 0.5 factors were used, and for the energy equation, an automatic algorithm embedded in the COMSOL Multiphysics software is found to be optimal.

3.2. Accomplishment of the computations The computations were performed for a rectangular container of 0.06  0.06  0.09 m. Three cases were examined (Table 1), for gallium with the properties listed in Table 2. The numerical solution was built using COMSOL Multiphysics software which uses the finite element method (FEM). Both 3D and 2D models were numerically explored to describe the liquid– solid interface by using the volume-of-fluid technique [20]. The grid was a built-in rectangular mesh. The grid size was chosen after careful examination of the solution convergence, the total cells number being 8932 for the 2D. The discretization of the conservation equations is done by Galerkin’s method [29] with Lagrange’s elements of the first order for the pressure and heat transfer equations and of the second order for the flow equation in the 2D model. The number of degrees of freedom (DOF) in the

4. Experimental setup In order to validate our theoretical model and numerical simulation, an experimental system was designed and constructed. Melting experiments were performed in a Perspex 0.06  0.06  0.09 m rectangular container (Fig. 4). Two vertical walls serving as a multipass heat exchanger were made of brass (5 in Fig. 4). The rest of the container walls were made of 0.01 m thick Perspex. Additional Perspex walls (3 in Fig. 4) were placed for better insulation. 99% metallic gallium was used, whose properties allow measurements in the temperature region close to room temperature. Temperature regimes were specified by Haake F3 laboratory thermostats (1 in Fig. 4) with the precision to ±0.5 °C. For temper-

Table 2 Gallium properties [33,34].

Gallium (solid) 298K < T < 303K Gallium (liquid) 303K < T < 323K

Property

Value

qs

5905 kg/m3 (60.66  0.183T + 6.03  104T2  7.136  107T3)W/m  K (9858.85 + 64.57T  0.1014T2)J/kg  K 303K (6273.98  0.605T + 4.82  105T2)kg/m3 (15.70 + 0.031T + 2.97  105T2)W/m  K (828.32  3.056T + 0.0079T2  9.13  106T3)J/kg  K (0.0156  1.053  104T + 3.198  107T2  5.005  1010T3)kg/m  s 80,200 J/kg  K

ks cps Tm

ql kl cpl

ll L

264

O. Ben-David et al. / International Journal of Heat and Mass Transfer 67 (2013) 260–271

Fig. 2. Comparison of solid–liquid interface profiles at different times and liquid fraction volume with time calculated by methods (a) and (b).

ature distribution measurements in the liquid and solid phases of the metal and for determining the liquid–solid interface position, eighteen calibrated T-type thermocouples (1/16-inch diameter stainless steel tubes) with ±0.3 °C accuracy were uniformly spaced along the length and height of the working container. Thermocouples readings were registered by the Agilent 34970A Data Logger

connected to PC. LabVIEW 2012 software of National Instruments was used for data processing. The interface position, its displacement with time, and one of the melt velocity components between thermally stabilized opposite vertical walls were recorded by ultrasonic Doppler velocimeter (UDV) DOP2000, Model 2150, Signal Processing, Lausanne, Switzerland [7]. Ultrasound Doppler velocimetry is based on the Doppler shift theory. The device sends an ultrasonic pulse into the fluid and receives an echo from particles of the fluid. The delay of the echo provides the distance between the transducer and a particle, and the Doppler shift reveals the particle velocity. It takes a short time to obtain a fluid velocity profile along the ultrasonic beam (Fig. 5). In spite of specific problems connected with such measurements (acoustic coupling, wetting conditions, medium reflecting capabilities etc.), they usually provide sufficient information on the mean flow velocity [25]. Five ultrasonic transducers TR0405RS (4 MHz) and TR1005RS (10 MHz) were placed into special slots in a hot plate 7 (Fig. 4) with a possibility to switch the transducers according to a specified order by a multi-channel device multiplexer. It was taken into account that transducers with higher emission frequency allowed to obtain a higher resolution, but with a concurrent growth of signal attenuation. Ultrasonic Doppler velocimetry provides a full profile of the flow velocity along the beam line (internal plot in Fig. 5). In our experimental situation, the divergence of the ultrasonic beam (1.8–2.5°) allows to estimate the lateral size of the measuring region close to the far cold wall depending on the distance from the transducers 2.8 mm to 3.8 mm. Before the experiment, the solid gallium is melted and poured into the experimental container through tube 8 (Fig. 4) at the top of the box. A tube 9 serves to let the air out. Both tubes also allow changing the volume during melting. Gallium solidifies by keeping a uniform initial temperature below melting temperature for several hours at the container using two heat exchangers. Melting is initiated by switching the left bath to a proper temperature exceeding the melting point. The left container wall has a constant temperature during the experiment. The temperature of the opposite wall is kept below the melting point, while other walls are thermally insulated. Reading and recording of the thermocouples and UDV transducers starts in parallel. The mentioned temporal advancement of the experimental liquid–solid interface is illustrated in Fig. 4. For its plotting, both the registered time-varied velocities profiles and peak positions of UDV power display were used. A great number of distinct points

Fig. 3. Comparison between the heat required for heating and melting the gallium, Qmelt, and the heat difference between the inlet and outlet, Qin  Qout, for two DT values (normalized to Qin  Qout maximal value).

O. Ben-David et al. / International Journal of Heat and Mass Transfer 67 (2013) 260–271

265

Fig. 4. Schematic diagram of the experimental system: 1 – constant temperature bath, 2 – plexiglas walls, 3 – plexiglas walls for insulated condition, 4 – thermocouples inside the container, 5 – heat exchanger, 6 – heat exchanger thermocouples, 7 – ultrasonic transducers, 8,9 – pouring and compensating tubes.

Fig. 5. Liquid–solid interface temporal progress for specific UDV transducer at the level y = 53 mm. Inset: Typical UDV display where the horizontal axis is the flow velocity [mm/s], and the vertical axis is the length of the container [mm].

266

O. Ben-David et al. / International Journal of Heat and Mass Transfer 67 (2013) 260–271

Table 3 Dimensionless parameters at the temperature TH for the melting process (K = 0.67 is aspect ratio, qm is the density at melting temperature 303 K, and ql, ml, al are liquid gallium properties at TH). Case

DT(oC)

PrT H

SteT H

RaT H

1 2 3

15 18 25

0.0247 0.0244 0.0242

0.07444 0.08933 0.12406

3.712  105 6.017  105 8.376  105

Fig. 6. Comparison of liquid–solid interface position at different times in GV experiment and in numerical 2D and 3D solutions at Thot = 38 °C and DT = 10 °C.

of liquid–solid interface were sampled, and a polynomial fitting was used to illustrate its temporal motion and profile. The maximal error of the points processing is described by R2 = 0.99.

5. Results and discussion The numerical modeling of the melting problem is extremely complicated, since besides the presence of two phases, there is also a moving liquid–solid interface and a mushy zone. Moreover, metal properties change with temperature during the process. The same reasons cause many impediments to experimental measurements in low-temperature metal. Note that our problem is solved in a dimensional form, but to make it possible to unify the obtained results and make general conclusions, the values of main dimensionless numbers are given in Table 3.

5.1. Comparison with experimental data 5.1.1. Solid–liquid interface position and profile Fig. 6 shows a comparison of liquid–solid interface position at different times in Gau and Viskanta experiments [2] and in our numerical models results for mid vertical x–y section. For the 2D model, the agreement is reasonable at the moment 120 s, but since the moment 360 s, the profile shape and its variation with time is different from the 3D case and the experimental results. For the subsequent process development, the error is larger in the upper part of the container (for example, at t = 1020 s the mean error is about 13% and the maximal error is about 23%). For the 3D model, the results agreement is better (the mean error is about 9% and the maximal error is about 20%).

Fig. 7. Comparison of liquid–solid interface position at different moments of time obtained from: – UDV experimental data, 2D and 3D computations for three cases: (a) DT = 15 °C, (b) DT = 18 °C, (c) DT = 25 °C.

A set of experiments was carried out to validate the developed numerical model. As noted above, liquid melt dynamics was studied by means of ultrasonic measurements. It is established that the lower boundary of the measured velocity is within the range down to 1 mm/s [7] and [9], and the estimation of measurements accuracy varies within the range of about

O. Ben-David et al. / International Journal of Heat and Mass Transfer 67 (2013) 260–271

267

Fig. 8. Temporal development of melting metal velocities at the point (45, 53), DT = 15 °C o – experimental points, – fitting of experimental data, 2D and 3D computations.

2–8% and even more. At the same time, it is noteworthy that the measurements accuracy depends on the experimental conditions. For its unbiased assessment, one should have the results of alternative measurements performed in a different way, or a standard case for which the velocity distribution is known (in time, along the coordinates, etc.). In the absence of such comparison tests, we have to be based on the described findings. We have also to note some problems with UDV measurement caused by melting peculiarities – transient character of the process under study, including thermal variations of medium properties, possible presence of the mushy zone with inconstant sound speed, gallium oxides appearance [25,33]. Moreover, varying conditions may be a cause of ultrasound pulses advancing failure and their breakdown, etc. Therefore, velocity profiles obtained using UDV require some kind of processing, and they were averaged, for the sake of comparison, by several approximation methods. The most convenient and simple among them is the polynomial fitting for spatial averaging. Time-averaged profiles were obtained from long-term recorded velocity profiles by choosing close-in-time values and accomplishing standard statistical procedures. For our experimental situation, most results were obtained using 128 emissions per profile. The sensitivity of 10 MHz transducers was acceptable at pulse repetition frequency (PRF) in the range from 160 Hz to 1 kHz to obtain velocities in the range from 0.012 m/s to 0.072 m/s. For 4 MHz transducers, the same PRF interval allows to get velocities from 0.028 m/s to 0.18 m/s. Fig. 7 displays a comparison of liquid–solid interface position at different moments of time between the ultrasonic experimental results and 2D and 3D numerical solutions for a vertical x–y section at a distance of 20 mm from the center. Experimental velocity values were processed using OriginPro 7.5 statistics tools, where the standard error of the mean value is characterized by R2 = 0.99. One can see a good agreement for the 3D solution and a less exact one for the 2D solution, meaning that the 3-rd dimension effects cannot be neglected, and a flow through the depth of the container has an essential effect on the heat transfer and the melting rate. The mean error for the 2D solution is 15–20.5% depending on the temperature difference DT, whereas the maximal one is about 40%. For the 3D solution the mean error is from 8.8% to 11.3%, and it also depends on the temperature difference, while the maximal error is about 14%.

Fig. 9. Longitudinal velocity profile temporal increase at the height y = 53 mm (600, 700, 800 s one after another); o – experimental points, parabolic approximation, 3D computations.

5.1.2. Liquid metal velocity Temporal behavior of the melting metal longitudinal velocity component at a fixed point (x = 45 mm, y = 53 mm) is illustrated in Fig. 8. As indicated, 3D method reveals a better matching with the velocity development. Here the experimental points were approximated using OriginPro 7.5 software by an exponential dependence of the 2-nd order with a reduced chi-squared value 1.58 and R2 = 0.82.

268

O. Ben-David et al. / International Journal of Heat and Mass Transfer 67 (2013) 260–271

Fig. 11. Temporal development of liquid–solid interface position, 2D model. Case 1: DT = 15 °C, maximal velocities vary from 8 to 15 mm/s.

Fig. 10. Temporal temperature change at the distance x = 40 mm for different heights: 15 mm, 30 mm, 40 mm; DT = 15 °C o – experimental points, 2D and 3D computations.

Fig. 9 shows velocity profiles comparison at the height y = 53 mm along the length of the container at different times in a vertical x–y section 20 mm from the center. Parabolic fitting of experimental data with R2 = 0.99 and SD = 0.25 demonstrates a good agreement with 3D computations.

5.1.3. Temperature distribution Temperature distribution can be a direct indication of the liquid–solid interface position. Fig. 10 shows a comparison of the temporal temperature change between the numerical models and the experiment in mid x–y section. A considerable temperature increase occurs both in computations and measurements when the liquid–solid interface passes over the point of interest. As expected, because of weak convective flow, one can see that temperature values fall with the level decrease. Here a better agreement exists between the measurements and 3D computations, too, which proves that 3D effects have a great influence on the melting process.

O. Ben-David et al. / International Journal of Heat and Mass Transfer 67 (2013) 260–271

269

Figs. 11 and 12 illustrate the solid–liquid interface temporal extension according to 2D numerical results in different cases, where solid phase is shown in blue, the liquid in red and the mushy zone in a transition range between blue and red. As expected, a convective flow is developing in time. First, the flow is negligibly weak, and the main heat transfer mechanism is thermal conduction, and thus the solid–liquid interface profile is nearly a straight line. Later, the convective flow becomes more significant [17]. Convection becomes the main heat transfer mechanism and has a great influence on the liquid–solid interface morphology. The liquid phase heats up at the hot wall, changes its density and flows upward. Then it flows toward the liquid–solid interface, cools and goes down. A comparison between Figs. 11 and 12 shows that the convection becomes more significant much earlier, and the volume flow is higher in case 2 because of a higher temperature difference. Accordingly, the liquid–solid interface profile becomes deformed and moves much faster. Melting conditions determine an interesting flow regime effect related to small rolls nucleation at the early stages of liquid phase flow (Fig. 13). These rolls are observed in 2D flows, and their dynamics can be analyzed within the framework of this model, while in 3D model (Fig. 14) a single roll dominates through the entire melting process. One can see in Fig. 14 the liquid–solid interface position at different times according to 3D numerical results, where the third dimension influence is considered. On the right side of Fig. 13 a top view of the melting process at the height y = 45 mm can be seen. The liquid–solid profile is curved due to the flow regime developing along the third dimension. The presence of side walls

Fig. 12. Temporal development of liquid–solid interface position, 2D model. Case 3: DT = 25 °C. Maximal velocities vary from 15 to 22 mm/s.

Comparative analysis of the obtained computational data with experimental ones (previous and described herein) reveals their acceptable fit. The developed numerical 3D model is suitable for the prediction of both the temporal dynamics of liquid–solid interface motion inside the container during melting, and melts flow in liquid phase including early stages of its formation. 5.2. 2D and 3D numerical simulation-flow analysis Computations give a possibility to distinguish the melt flow details from the early stages including rolls birth and its further development.

Fig. 13. Rolls nucleation and transformation at early flow stages, time from 15 to 40 s, case 3: DT = 25 °C.

270

O. Ben-David et al. / International Journal of Heat and Mass Transfer 67 (2013) 260–271

Fig. 14. Temporal development of liquid–solid interface position, 3D model. Case 3: DT = 25 °C. Left side: midplane of vertical cross-section with maximal velocities from12 to 14 mm/s; right side: midplane of horizontal cross-section with maximal velocities from 1.2 to 6.1 mm/s.

and non-slip boundary conditions influences the flow regime along the third dimension, which is absent in the 2D model. The conditions of our computations, temperature regimes and container dimensions were chosen in compliance with those described by Gau and Viskanta [2]. But even in this particular case, the influence of the side walls on the process is confirmed by calculations and indicates that 3D models are necessary for similar problems. A comparison between Figs. 12 and 13 displays that in the 3D numerical model liquid–solid interface advance is less dynamic, and its volume flow is less intensive. Another difference appears at the early stages of flow formation. In the 2D model, flow origin

is accompanied with small rolls development. Their quantity definitely depends on the process parameters and changes during the melting progress. On the contrary, in the 3D model a single roll is formed in the beginning, enlarging during the melting process. Alexiades et al. [21] noted that several studies using finite element approach in order to predict gallium melting in a rectangular container got different solutions. One of them found that during a rather long time, the flow structure consists of a number of small rolls, while another found only a single roll. It seems that in the present study the 2D and 3D numerical models provide both solutions, but a physical explanation has to be developed.

O. Ben-David et al. / International Journal of Heat and Mass Transfer 67 (2013) 260–271

For example, in the 2D model case 1 (DT = 15 °C) one can see four convective rolls in the liquid phase at the early periods (until 100 s).Throughout 300 s, the rolls quantity is reduced to three, later (at 500 s) there are two rolls left, and after that (1020 s) a single big roll remains. In case 3 (DT = 25 °C) the rolls confluence process is much faster. At the beginning only three rolls are present in the flow, after that a reduction to a pair of rolls can be seen (140 s), and one big roll remains during the further process development (about 500 s). Beyond any doubt, the temperature difference between the hot and cold walls has a great influence on the flow regime and its intensity, including the rolls quantity. Detailed analysis of this dependence is not executed here, but we emphasize the numerical model capabilities for the study of the effect of geometrical ratios, boundary conditions and temperature gradient on the rolls quantity. A certain transition region (mushy zone) can be also seen in both 2D and 3D computations. The width of this region at the lower side of the container is much wider due to a lower temperature gradient. 6. Conclusions A numerical and experimental study of solid–liquid phase transition in a rectangle container is performed. An advanced numerical model was developed for the computations of metal melting in a rectangular container with a heated vertical wall. The thermal dependences of all material properties were taken into account at the calculation. The Carman–Kozeny constant values were examined, and their optimal value was found for the prediction of the liquid–solid interface behavior. The error of solid liquid interface position determination grows with the temperature difference, which can be explained by the way we calibrate the numerical models according to experimental results. Several different values of the constant C were examined for such calibration of the first case, DT = 15 °C, and C = 1011 was found to be the most accurate. In other cases, the same constant was used. The growing error may indicate to a certain dependence between the constant C and the temperature difference. According to 2D computations, the melt flow origin is accompanied by small rolls nucleation. Their quantity depends on the process parameters and reduces during the melting progress in such a way that one big roll remains in the further process development. Calculations reveal the mushy zone appearance close to the liquid–solid interface. As expected, this zone is wider in the lower part of the container because of smaller temperature gradients. A set of experimental data in liquid phase was obtained using the ultrasonic Doppler velocimeter. These results, together with collateral temperature measurements, describe liquid–solid interface dynamics and temporal behavior of longitudinal mean flow velocity component. The local velocities obtained by ultrasonic measurements were compared with those predicted by the numerical simulations, both 2D and 3D. The comparison demonstrates a better agreement between the experimental data and the results of the 3D numerical simulation, which confirms a great influence of the third dimension. It was shown that 3D computations reveal an acceptable fit with experimental data and can be exploited as an effective instrument for the prediction of the liquid–solid interface dynamics and melt flow features, starting from the early stages of flow development. The obtained data require a certain expansion. Nevertheless, they obviously demonstrate the potentialities of both the

271

described numerical 3D model and measuring method as effective tools for the study of similar problems. References [1] C. Gau, R. Viskanta, Effect of natural convection on solidification from above and melting from below of a pure metal, Int. J. Heat Mass Transfer 28 (1984). [2] C. Gau, R. Viskanta, Melting and solidification of pure metal on a vertical wall, Heat Transfer 108 (1986) 174–181. [3] T.A. Campbel, J.N. Koster, Visualization of liquid–solid interface morphologies in gallium subject to natural convection, Cryst. Growth 140 (1994) 414–425. [4] A.S. Menon, M.E. Weber, A.S. Mujumdar, The dynamics of energy storage for paraffin wax in cylindrical containers, Can. J. Chem. Eng. 61 (1983) 643–653. [5] D. Pal, Y.K. Joshi, Melting in a side heated tall enclosure by a uniformly dissipating heat source, Int. J. Heat Transfer 44 (2001) 375–387. [6] L. Katsman, V. Dubovsky, G. Ziskind, R. Letan, Experimental investigation of solid–liquid change in cylindrical geometry, in: Thermal Engineering Summer Heat Transfer Conference, Vancouver, Canada, 2007. [7] http://signal-processing.com. [8] P. Oborin, I. Kolesnichenko, Study of liquid metal flow and crystallization in a plane later with stirring, in: Conference on Theoretical and Applied Magnetohydrodynamics, Perm, Russia, 2012, p. 75. [9] A. Cramer et al., Liquid metal model experiments on casting and solidification processes, J. Mater. Sci. 39 (2004) 7285–7294. [10] N. Ramachandran, J.R. Gupta, Y. Jalunu, Thermal and fluid flow effects during solidification in a rectangular cavity, Int. J. Heat Mass Transfer 25 (1982) 187– 194. [11] A. Gadgil, D. Gobin, Analysis of two dimensional melting in rectangular enclosures in the presence of convection, J. Heat Transfer 106 (1984) 20–26. [12] M.R. Albert, K. O’Neill, Transient two-dimensional phase change with convection using deforming finite elements, Int. Comput. Tech. Heat Transfer 1 (1985). [13] M. Vynnycky, S. Kimura, An analytical and numerical study of coupled transient natural convection and solidification in a rectangular enclosure, Int. J. Heat Mass Transfer 50 (2007) 5204–5214. [14] K. Morgan, A numerical analysis of freezing and melting with convection, Comput. Meth. Appl. Eng. 28 (1981) 275–284. [15] D.K. Gartling, Finite element analysis of convective heat transfer problems with change of phase, in: K. Morgan et al. (Eds.), Computer Methods in Fluids, Pentech, London, 1980, pp. 257–284. [16] V.R. Voller, M. Cross, N.C. Markatos, An enthalpy method for convection/ diffusion phase change, Int. J. Numer. Methods Eng. 24 (1987) 271–284. [17] C. Prakash, V.R. Voller, A fixed grid numerical modelling methdology for convection–diffusion mushy region phase-change problem, Int. J. Heat Mass Transfer 30 (1987) 1709–1719. [18] A.D. Brent, V.R. Voller, K.J. Reid, The enthalpy-porosity technique for modeling convection–diffusion phase change: application to the melting of a pure metal, Numer. Heat Transfer 13 (1988) 297–318. [19] W.D. Carrier III, Goodbye, Hazen: Hello, Kozeny–Carman, J. Geotech. Geoenviron. Eng. (2003) 1054–1056. [20] H. Shmueli, G. Ziskind, R. Letan, Melting in a vertical cylindrical tube: numerical investigation and comparison with experiments, Int. J. Heat Mass Transfer 53 (2010) 4082–4091. [21] V. Alexiades, N. Hannoun, T.Z. Mai, Resolving the controversy over Tin and Gallium melting, Numer. Heat Transfer 44 (Part B) (2003) 253–276. [22] J.A. Dantzig, Modeling liquid–solid phase changes with melt convection, Int. J. Numer. Methods Eng. 28 (1989) 1769–1785. [23] F. Stella, M. Ginagi, Melting of a pure metal on a vertical wall: numerical simulation, Numer. Heat Transfer Appl. 38 (Part A) (2000) 193–208. [24] Y. Takeda, Ultrasonic Doppler method for velocity profile measurement in fluid dynamics and fluid engineering, Exp. Fluids 26 (1999) 177–178. [25] D. Brito, H.C. Nataf, P. Cardin, J. Aubert, J.P. Masson, Ultrasonic Doppler velocimetry in liquid gallium, Exp. Fluid 31 (2001) 653–663. [26] C. Beckermann, R. Viskanta, Natural convection solid/liquid phase change in porous media, Int. J. Heat Transfer 31 (1987) 35–46. [27] B.W. Webb, R. Viskanta, Analysis of heat transfer during melting of a pure metal from an isothermal vertical wall, Numer. Heat Transfer 9 (1986) 539– 558. [28] N.S. Cheng, Application of Ergun equation to computation of critical shear velocity subject to seepage, J. Irr. Drain. Eng. ASCE 129 (4) (2003) 278–283. [29] O.C. Zienliewicz, P. Nithiarasu, R.L. Taylor, Finite Elem. Methods Fluid Dyn. 6 (2005). [30] P.N. Brown, K.E. Grant, S.L. Lee, R. Serban, D. EShmaker, A.C. Hindmarsh, C.S. Woodward, Sundials: suite of nonlinear and differential/algebraic equation solver, ACNT Math. Softw. 31 (2005) 363. [31] O. Schenk, K. Gartner, PARDISO User Guide, 2011, http://www.pardisoproject.org. [32] P. Amestoy et al., 2011. http://www.graal.ens-lyon.fr/MUMPS. [33] Y. Tasaka, Y. Takeda, T. Yanagisawa, Ultrasonic visualization of thermal convective motion in a liquid gallium layer, Flow Meas. Instrum. 19 (2008) 131–137. [34] D.R. Lide (Ed.), CRC Handbook of Chemistry and Physics, 84th ed., CRC Press, 2003.