Journal of Fluid Flow, Heat and Mass Transfer (JFFHMT)

ISSN: 2368-6111

Nomenclature:

Latin Symbols
	Mushy zone constant	kg m−3 s−1
	Specific heat capacity	J kg–1 K–1
	Gravity vector	m s−2
	Specific enthalpy	J kg–1
	Latent heat	kJ kg–1
	Pressure	Pa
	Buoyancy source term	kg m−2 s−2
	Darcy type source term	kg m−2 s−2
	Surface tension force term	kg m−2 s−2
	Surface area	m
	Time	s
	Temperature	K
	Velocity vector	m s–1
Greek Symbols
	Volume fraction	-
	Thermal diffusivity	m2 s−2
	Coeff. Of thermal expansion	K–1
	Melt fraction	-
	Thermal conductivity	W m–1 K–1
	Eigenvalue	-
	Density	kg m–3
	Dynamic viscosity	kg m–1 s–1
Subscripts
c	Cold
Ref	Reference
liq	Liquidus
sol	Solidus
m	Melt
h	Hot
Abbreviations
CFD	Computational Fluid Dynamics
DSLR	Digital Single Lens Reflex
LED	Light Emitting Diode
PCM	Phase Change Material
RT	Rubitherm®
VoF	Volume of Fluid

1. Introduction

The solid-liquid phase change in the presence of third phase (liquid or gas) is of the utmost importance in many industrial processes such as casting, welding, thermal storage systems, powder atomization, etc. Therefore, the understanding of solidification and melting, its initial and thermal conditions plays a crucial role to make the process more efficient.  Over the past few decades numerical simulations of solid-liquid phase change is gaining popularity in the industry not only to improve the quality of their products, but also to make cost effective manufacturing decisions [1]. A common approach for modeling such multi physical phenomena is coupling of Volume of Fluid (VoF) method with a solid-liquid phase change model in Eulerian frame.

Since the introduction of the enthalpy method by Voller et al. [2], several numerical methods, formulations and models have been found in the literature to tackle the solid-liquid phase change phenomena and especially in handling of the velocity field in the solid phase. Hu and Argyropoulos [3] reviewed various methods and stated the pros and cons of the discussed methods according to the nature of the problem and computational expenses. Rösler and Brüggemann [4] developed a numerical model to study the melting process of phase change material (PCM) in latent heat storage systems and incorporated an error function to describe the melt fraction. Miehe [5] implemented and compared the viscosity and Darcy model to investigate the twin-roll casting process of the magnesium alloy. Kasibhatla et al. [6] adopted the variable viscosity model to study the melting and settling of PCM in a capsule. More recently König-Haagen et al. [7] performed a comprehensive benchmarking of fixed-grid methods and presented a qualitative and quantitative analysis of most common energy formulations for the modeling of melting.

The VoF approach proposed by Hirt and Nichols [8] is the well-known interface-capturing technique, which uses an indicator function to mark different fluids. This approach has also been previously used for the solid-liquid phase change simulations in a multiphase system. Assis et al. [9] performed parametric CFD analysis of melting process of PCM (Rubitherm® RT27) in a partially filled spherical shell. The VoF method was employed to account for PCM-air system and the enthalpy-porosity method was adopted for the solid liquid phase change. A temperature dependent density function was used for air phase, whereas a variable density based on Boussinesq approximation was implemented in the liquid PCM.

Shmueli et al. [10] investigated the influence of the large mushy zone constant  in vertical cylindrical tube, partially filled with PCM (RT27) with a flat surface. A temperature dependent viscosity function and variable density  were used in the liquid PCM. Two sets of thermophysical properties were examined, at first uniform values of specific heat capacity ( and thermal conductivityfor solid and liquid phase was assigned, later a small variation of thermal conductivity between phases was considered (i.e.). Rösler [11] investigated close contact melting of a PCM with air above it in a capsule. The PCM was assumed to be an incompressible fluid, whereas the gas phase (air) was modeled as a compressible fluid. Similarly, Kasibhatla et al. [12] also investigated the close contact melting in a capsule using a multiphase approach and variable viscosity model to enforce the velocity in the solid region to zero.

Furthermore, Tan et al. [13] performed simulations in small scale geometries to investigate the effects of thermocapillary convection and the dynamics of the gas phase during the solidification of pure bismuth with an argon gas on the top using a multiphase approach. Saldi [14] studied the Marangoni driven free surface hydrodynamics in liquid steel weld pools during solid-liquid phase change with free surface deformation. Finally, Richter et al. [15] performed mold filling and solidification simulations; the gas phase was treated as an incompressible fluid with constant properties. A special treatment was employed to handle the Darcy-type source term at the interface between PCM and air.

As evident from the above literature research, the coupling of VoF Method with the enthalpy porosity method is an effective way to tackle such multi-physical phenomena numerically. Therefore, in this study a fundamental numerical model of solid-liquid phase change in the presence of gas phase is presented. The PCM and gas phase are considered as two-phase immiscible incompressible fluid and are modeled using the VoF [8] approach based on interface capturing method. A fixed grid enthalpy-porosity method [16, 17] is employed for the solid-liquid phase change. Thermophysical properties of both the fluids are constant and independent of temperature. The Boussinesq approximation is used for both the fluids (liquid PCM and air).

2. Numerical Model

2. 1. Governing Equations

The solid-liquid phase change with a free surface is handled by solving the governing equations discussed in this section, using a single domain formulation. The conservation equations accounts for both the differences in fluid properties as well as surface tension forces acting at the interface. The mass, momentum and energy conservation equations are given by [13, 14, 16, 17]:

where   is the velocity vector, t the time, p the pressure,  the density,  the dynamic viscosity,  the specific enthalpy, T the temperature,  the thermal conductivity,  the buoyancy source term, the surface tension force term,  the Darcy type source term.

2. 2. Volume of Fluid (VoF) method

The VoF uses an indicator function  to mark different fluids:

in the present work fluid A is the PCM and fluid B is the gas phase. The volume fraction α transport equation can be written as:

where, is the relative velocity between two fluids. On the contrary to the classical VoF method originally proposed by Hirt and Nichols [8], the VoF method in OpenFOAM® uses an additional convective term also referred to as the artificial compression term (the third term in Eq.(5)) to minimize the numerical diffusion at the interface [18]. This term is only active within the cells containing interface due to its dependence on where  (0, 1). The relative velocity also known as the compression velocity is computed by:

where  is the fluid flux, surface area and  compression coefficient. In the present work is set to 1. The thermophysical properties of each fluid are expressed in terms of volume fraction:

where  is either.

The surface tension force term  in Eq. (2) is the continuum surface force (CSF) model proposed by Brackbill et al. [19] and is expressed as:

where  is the surface tension and the surface curvature. The surface curvature is defined as:

2. 3. Solid-liquid phase change

The specific enthalpy  in Eq. (3) can be defined as the sum of sensible heat and latent heat, such that:

where is the specific heat,  the reference temperature,  restricts the term in fluid A (PCM),  the melt fraction as a function of temperature and the latent heat of fusion.

By combining Eq. (3) and Eq. (10), we get the energy conservation equation in terms of temperature, as shown below:

The melt fraction  [0, 1] is defined as a linear function of temperature:

where   is the solidus temperature and   the liquidus temperature.

 

The Darcy type source term  in Eq. (2) is modeled using the Carman-Kozeny equation for flows in porous media [16] and is defined as:

 

where  is the large constant that describes the morphology of the melting front in the mushy zone and  = 1  10^-3 is the small computational constant to avoid "divide by zero" error. In Eq. (13), when  = 0 the source term dominates in the cells containing solid phase and enforces the velocity in that cell to zero, and when  = 1 the source term vanishes in the cells containing liquid phase.

 

The thermophysical properties of each phase are expressed in terms of melt fraction^[1], as shown in Eq. (14), where  is either  or.

The buoyancy source term in Eq. (2) is described by the Boussinesq approximation, such that:

where  is the gravity vector.  

Figure 1 illustrates an overview of the complete model and is implemented in the open source FVM based CFD software OpenFOAM®. For the pressure-velocity coupling, the PIMPLE (PISO and SIMPLE) algorithm is employed. The VoF method is solved explicitly using a modified FCT (Flux Corrected Transport) method called MULES (Multidimensional Universal Limiter with Explicit Solution) and is based on the work of Zalesak [20] with an iterative calculation of the weighting factors [21]. A first order implicit Euler scheme is employed for the temporal discretization. For the gradient terms a second order Gauss linear scheme is applied. A variety of HRS (High Resolution Schemes) discretization schemes are used to solve the divergence terms of the equations.

3. Validation

3. 1. Stefan problem

The one-dimensional one-phase Stefan problem is a well-known transient heat transfer problem for directional solidification [1]. Figure 2 shows the schematic diagram of the problem.

The fluid with a melting point of  = 882.15 K is used. The left and right boundaries are kept at constant temperatures  = 822.15 K () and  = 882.7 K () respectively. The initial temperature of the system is  = 882.7 K and the thermophysical properties of the fluid are tabulated in Table 1. The analytical solution for the determination of interface location  at time t is given by:

where  is the thermal diffusivity, and  is the non-linear eigenvalue and is determined by solving the below transcendental equation using iterative methods.

Table 1. Thermophysical properties of the fluid.

Properties	Values
[kg m–3]	1780
[kg m–1 s–1]	1 x 10–3
[W m–1 K–1]	80
[J kg–1 K–1]	1050
[kJ kg–1]	339

In the simulations a computational domain of size 150 mm x 2.5 mm (L  B) is selected and discretized with equidistant cell size of x = 0.5 mm, resulting in 1500 cells. The top and bottom boundaries are considered as adiabatic; an empty^[2] boundary condition is applied on the front and back boundaries of the domain.

At time t = 0, solidification commences from the left boundary and continues to propagate towards the right boundary.  Figure 3 shows a comparison of numerical simulation and analytical solution for the one-dimensional Stefan problem. The interface locations predicted by the simulation are in excellent agreement with the analytical solution.

3. 2. Melting of gallium

The melting of gallium in a rectangular enclosure is a well-documented solid-liquid phase change exp-eriment performed by Gau and Viskanta [22]. The computational domain is shown in Figure 4, the size of the domain is 88.9 mm x 63.5 mm and is divided into 42 x 32 hexahedral numerical cells.

The left and right walls are maintained at constant temperatures  = 311.15 K,  = 301.45 K respectively. The top and bottom walls are assumed to be adiabatic. The initial temperature of the system is  = 301.45 K. The thermophysical properties of gallium used in the simulations are given in Table 2 Furthermore, the mushy zone constant C and the computational constant ε in Eq. (13) are set to 1.6  10⁶ kg m⁻³ s⁻¹ and 1  10⁻³ respectively, see Brent et al [16].

Table 2. Thermophysical properties of gallium.

Properties	Values
[kg m–3]	6093
[kg m–1 s–1]	1.81 x 10–3
[W m–1 K–1]	32
[J kg–1 K–1]	381.5
[kJ kg–1]	80.160
[K-1]	1.2 x 10-4

The temporal evolution of the melting front obtained from simulation exhibits good agreement with experiments of Gau and Viskanta [22] as shown in Figure 5.

The slight deviations between the simulation and experiments can be explained as follows:

a.       the melting fronts are obtained from two-dimensional simulation, whereas the experiments performed by Gau and Viskanta [22] were three-dimensional (for further information refer [5, 23]),

b.       temperature independent thermophysical properties are used in the simulation,

c.       the experiment was stopped multiple times and the melt was drained in order to track and record the melting front, which may have influenced the results.

3. 3. Free surface thermocapillary flow with phase change

A free surface thermocapillary flow with phase change in a two-dimensional rectangular container under microgravity conditions investigated by Tan et al. [13] is used for validation. For the sake of validation, the source term  in Eq. (8) is resolved into the normal  and the tangential  components [13, 14, 24] and is expressed by:

Figure 6 shows the computational domain of size 15 mm x 5 mm and is discretized into 140 x 62 hexahedral cells [13, 14]. The bismuth melt is filled up to height  = 4 mm and the rest is filled with argon gas. The left and right walls are kept at constant temperatures  = 552.55 K and  = 540.55 K respectively. The top and bottom walls are enforced with a linear temperature distribution from  to   , such that, the isoline of melting temperature  = 544.55 K is at a distance  from the left wall. The dimensionless numbers corresponding to bismuth and argon gas are listed in Table 3[13]

Table 3. Dimensionless numbers.

Symbol	Values
	244
	0.031
	0.0022
	0.033
	1.88 x 10–4
	0.019

where  is the Marangoni number,  Rayleigh number,  Capillary number,  Stefan number,  Biot number and  Prandtl number.  The comparison of solid-liquid interface of the present work, the simulations by Tan et al. [13], and the simulations by Saldi [14] is shown in Figure 7. An excellent agreement is achieved in comparison with literature.

4. Experimental setup

Experiments are conducted to further investigate the introduced model. Figure 8 illustrates the schematic of the experimental setup adapted from [11, 22, 26]. The domain consists of a square enclosure with inside dimensions of 50 mm x 50 mm x 50 mm made of Plexiglas. The left and right walls are made of aluminum plates which serve as heat source and heat sink respectively. A couple of multi-purpose riser pipes are screwed on the top wall of the enclosure which allows free flow of air to compensate for volume expansion/contraction during phase change. The enclosure is initially filled with paraffin wax (RT31) up to a height of 40 mm in liquid form through one of the riser pipes and is allowed to solidify for a certain period of time. All the walls are insulated except for the front and back walls. A DSLR camera is placed at an appropriate distance from the front wall to capture the sharp images of the temporal evolution of the melting front. A LED diffuser is placed at a suitable distance from the back wall in order to differentiate the solid and liquid phase during the melting process of the wax. Finally, all the images are processed and automated using ImageJ [25] for the evaluation of melting fronts and melt volume fraction of the paraffin wax.

5. Results and discussions

The thermophysical properties of RT26 to RT35 provided by the manufacturer [27] are comparable, the specific heat capacity and the thermal conductivity of the above mentioned PCMs are same and have uniform values for each (solid/liquid) phase. Shmueli et al. [10] compared their numerical results with experiments using the thermophysical properties provided by the manufacturer with a slightly varied thermal conductivities for solid and liquid phase. The thermophysical properties used in this study are adapted from the above literature and are tabulated in Table 4.

Table 4. Thermophysical properties used in CFD Simulations

Properties	Rubitherm® RT31		Air
	Solid	Liquid
[kg m–3]	880	760	1.1885
[kg m–1 s–1]	-	3.08 x 10–3	1.81 x 10–5
[W m–1 K–1]	0.24	0.15	0.0258
[J kg–1 K–1]	2000	2000	1006
[kJ kg–1]		160	-
[K-1]		7.6 x 10-4	3.38 x 10-3

The grid convergence study is performed to investigate the influence of cell size on the melt fraction. Different grid resolutions of 2D computational domain have been investigated, from coarse (200 x 200 x 1) cells to fine (400 x 400 x 1) cells increased in steps of 25 cells in each x and y direction. A grid resolution of (275 x 275 x 1) is selected in this study. The relative error of the melt fraction between the selected grid to the finest grid resolution is less than 2%. An adaptive time stepping with Courant Number () of 0.75 is adopted.

Figure 9 shows the experimental and numerical, temporal evolution of the melting front at 15, 30, 45 and 60 min respectively, for the temperature at hot and cold walls  = 313.15 K and  = 295.15 K respectively. The liquid paraffin close to the hot wall rises upwards and flow horizontally towards the cold melting front, where it releases heat at the upper region of the solid paraffin and descends along the cold melting front, thereby advancing the melting front rapidly in the upper region than in the lower region. The streamlines in the liquid paraffin shows that there is no movement approximately at the center of the melt, this is because of the upwards and downwards flow of liquid paraffin close to hot wall and cold melting front (indicated by solid arrows in the bottom row of Figure 9). It is observed that the air flow above liquid paraffin also progresses along with the melting front in a circular motion^[3]. A secondary air flow is also observed in the region above the solid paraffin. However, the strength of this secondary air flow is minuscule. In the experiment it is notable that the liquid paraffin-air interface is thicker than in simulations. This is because of the contact angle between liquid paraffin-air interface and the Plexiglas wall, which leads to the formation of meniscus near the Plexiglas wall. Furthermore, the melting fronts in the simulation appear to be sharper than in the experiment near the paraffin-air interface. This is because in experiments the melting fronts are obtained by the evaluation of photographs taken through front Plexiglas wall. Nevertheless, it is noteworthy to mention that the surface tension effects are neglected during simulations.

Figure 10 shows the comparison of melting fronts between experiments and simulation. The height of the liquid paraffin in the experimental domain varies due to the volume expansion during phase change. Therefore, for simplicity the temporal evolution of melting fronts has been normalized to the width of the experimental domain in x direction and to the height of the liquid paraffin in y direction. The melting fronts obtained by the numerical model are in good agreement with experiments.

 However, the melting fronts at the lower region advances faster in experiments than in simulations, this is because the thermal conductivity of the Plexiglas walls also has an effect on the melting of the paraffin, however the influence of the Plexiglas walls is not considered during simulations.

6. Conclusions

In this study, a fundamental numerical model to tackle the solid-liquid phase change phenomena in the presence of gas phase has been investigated by means of well documented case studies in the literature and experimental validation. The model is based on the use of VoF approach to account for PCM and gas phase and enthalpy-porosity method to account for solid-liquid phase change in the PCM. The model has been validated using multiple references cases. The simulated results are in good agreement with the results described in the literature.

   To further investigate the limitations of the model, experiments with paraffin wax (RT31) are conducted. Experiments are performed in a square enclosure partially filled with paraffin wax in the presence of air. The left and right walls of the enclosure are maintained at temperatures  = 313.15 K and   = 295.15 K respectively, whereas the remaining walls are assumed to be adiabatic. The temporal evolution of the melting front obtained from simulation show a good agreement with the experiments. The small discrepancies in the results can be accounted for, as follows:

• In experiments it is observed that the side Plexiglas walls also influence the melting phenomena within the enclosure, whereas in simulations these walls were not considered and are modeled as perfectly insulated.
• Temperature independent thermophysical properties are used in the simulations for PCM (solid/liquid) and gas phase.
• Only 2 dimensional simulations are carried out and the influence of the third dimension is not yet investigated.
• Finally, the errors in evaluation of photographs obtained from experiments due to volume expansion and surface tension effects.

Acknowledgements

The authors wish to acknowledge the Federal Ministry of Economic Affairs and Energy (BMWi), Germany for the financial support of MOTEVAS bearing project number 03ET1192B. Special thanks to the Computer Centre (URZ) of TU Bergakademie Freiberg for the HPC-Cluster support.