- Email: [email protected]
CFD simulation of melting process of phase change materials (PCMs) in a spherical capsule
Accepted Manuscript Title: CFD simulation of melting process of phase change materials (PCMs) in a spherical capsule Author: H. Sattari, A. Mohebbi, M.M. Afsahi, A. Azimi Yancheshme PII: DOI: Reference:
S0140-7007(16)30287-0 http://dx.doi.org/doi: 10.1016/j.ijrefrig.2016.09.007 JIJR 3421
To appear in:
International Journal of Refrigeration
Received date: Revised date: Accepted date:
15-2-2016 16-7-2016 5-9-2016
Please cite this article as: H. Sattari, A. Mohebbi, M.M. Afsahi, A. Azimi Yancheshme, CFD simulation of melting process of phase change materials (PCMs) in a spherical capsule, International Journal of Refrigeration (2016), http://dx.doi.org/doi: 10.1016/j.ijrefrig.2016.09.007. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
CFD simulation of melting process of phase change materials (PCMs) in a spherical capsule H. Sattaria, A. Mohebbia,*, M. M. Afsahia, A. Azimi Yancheshmeb a
Department of Chemical Engineering, Faculty of Engineering, Shahid Bahonar University of Kerman, Kerman, Iran b
Chemical and Petroleum Engineering Department, Sharif University of Technology, Tehran, Iran *Corresponding author: Tel & Fax: +983432118298 E-mail: [email protected], [email protected]
Highlights The effect of container’s surface temperature on melting rate is more than other parameters. At the beginning of the melting process of PCM, the conduction mode is governing. Natural convection is the ruling mode at the top half of container. The molten liquid at the top half of the container is naturally stable. Liquid layer, which is formed at the bottom of the container is unstable.
Abstract The present study is focused on CFD simulation of constrained melting of Phase Change Materials (PCMs) in a spherical container. To investigate the melting process of the PCM, its melting fraction was analyzed at different 1 Page 1 of 31
times. The results indicated the existence of thermally stable layers on the top of the sphere. Moreover, inspection of the calculated temperatures at different points along the vertical axis indicates the existence of some disturbances at the bottom of the sphere due to the natural convection. After the validation of the results, the effects of different parameters such as the surface temperature of the capsule, the initial temperature and the size of the spherical capsules, on the melting process were investigated. The initial temperature did not affect the melting rate, whereas melting rate increased by increasing the surface temperature of the capsule and also decreasing the diameter of the sphere. The results showed that the surface temperature of capsule compared to geometrical parameters and other operational conditions can have a greater influence on the melting rate and the heat flux. Keywords: Constrained melting; Phase change materials; Computational fluid dynamics (CFD); Spherical container
1. Introduction Recently, latent heat thermal energy storage systems have been widely applied in cool storage for central air-conditioning, hot storage for air heating, solar energy, energy efficient buildings and waste heat recovery. Storage 2 Page 2 of 31
systems allow for the efficient and rational utilization of the available resources or renewable energies, by using the time lag between production, as well as the availability of the energy and its consumption in the demanding systems. Phase change materials (PCM) are mainly used to provide higher storage densities. As a serious problem PCMs have relatively low thermal conductivities. Therefore, increasing the surface/volume ratio of PCMs in order to increase the heattransfer rate is very appealing. This can be done by packing a volume with a great number of PCM capsules. The spherical geometry is one interesting case for heat storage applications, since spheres are much employed in packed beds. Due to the complexity of such systems, it is often more effective to first model the behavior of an individual sphere and then describe it with a simple parametric model in the packed bed modeling. Roy and Sengupta [1] examined the melting process with the solid phase initially uniformly super cooled. In order to include the effects of a temperature gradient in the solid core, they modified the heat transfer equation. At each time step, the unsteady conduction equation has been solved numerically using a toroidal coordinate system with a suitable immobilization of the moving boundary to transform the infinite domain into a finite one. In another paper [2], they investigated the outcome of natural convection on the melting process. In order to reduce the computational effort and the time, they made some simplifications. They obtained the non-dimensional melting time and the heat
3 Page 3 of 31
transfer coefficient as functions of the PCM property values, the operating temperature and the physical size. Ettouney et al. [3] experimentally evaluated the heat transfer (during energy storage) and release for the phase change of paraffin wax in the spherical shells. They showed that an increment in the Nusselt number of the sphere with a larger diameter is attributed to the increase in the natural convection cells in the PCM inside the sphere. The natural convection role is enhanced upon an increase in the sphere diameter and the air temperature. On the other hand, during solidification the wax layers are formed inward, nucleating on the sphere walls. As solidification progresses, the melt volume becomes smaller and the role of natural convection diminishes rapidly. Khodadadi and Zhang [4] also considered the effects of buoyancy-driven convection on the constrained melting of phase change materials within spherical containers. They reported that during the initial melting process, the conduction mode of the heat transfer is dominant. As the buoyancy-driven convection is strengthened due to the increase in the melting zone, the process at top region of the sphere is much faster than at the bottom due to the increment of the conduction mode of heat transfer. They found that buoyancydriven convection speeds up the melting process compared to the diffusioncontrolled melting. Assis et al. [5] numerically and experimentally explored the melting process of a PCM in a spherical geometry. The results of their experimental 4 Page 4 of 31
investigation, including visualization, were favorably comparable with the numerical results and therefore suitable for validation of the mathematical approach. Their computational results showed how the transient phase-change process depends on the thermal and geometrical parameters of the arrangement. Veerappan et al. [6] investigated the phase change behavior of 65 mol% capric acid and 35%mollauric acid, calcium chloride hexa hydrate, n-octadecane, nhexadecane, and n-eicosane inside the spherical containers to identify a suitable heat storage material. They created analytical models for solidification and melting of PCM in a spherical shell with conduction, natural convection, and heat generation and found a good agreement between the analytical predictions and the experimental data. Both models were validated with the experimental work of Eames and Adref [7]. Regin et al. [8] examined the heat transfer performance of a spherical capsule using paraffin wax as PCM placed in a convective environment during the melting process. The model results were in a good agreement with the experimental data. Tan [9] investigated the melting of the phase change material (PCM) inside a sphere using n-octadecane for both constrained and unconstrained melting processes. For constrained melting, paraffin wax (n-octadecane) was immobilized, through the use of thermocouples when melting is done inside a transparent glass sphere. Their experimental setup provided a detailed temperature data that were gathered along the vertical diameter of the sphere 5 Page 5 of 31
during the melting process. Tan and Khodadadi [10] also experimentally and numerically investigated the constrained melting of PCMs inside a spherical capsule to understand the role of the buoyancy-driven convection. In this study, a computational fluid dynamics (CFD) modeling on the constrained melting of phase change material (PCM) in a spherical container was performed. The effects of different parameters like the capsule size, the surface temperature, the initial temperature and the Stefan number, on the PCM melting process were investigated.
2. CFD simulation CFD simulation was carried out using the commercial Fluent 6.3 software. The physical model and computational procedure are discussed below.
2.1. Physical model Consider a spherical shell with an internal radius (Ri) of 50.83 mm and a wall thickness of 1.5 mm initially filled with a solid PCM at an initial temperature (Ti) lower than the melting temperature (Tm).The schematic diagram of the physical model is shown in Fig. 1. At t (time) > 0, a constant temperature (T0) greater than the melting temperature is imposed on the surface of the sphere, i.e. T0>Tm. Melting starts at the surface, moving the solid-liquid interface towards the middle of the field. Since the present work is focused on
6 Page 6 of 31
the analysis of constrained melting, both the solid and liquid phases have the same density. The container was made of glass with a thermal conductivity of 1.14 Wm-1K-1. The properties of the PCM, based on n-octadecane as a commercially available material, are given in Table 1.
For the experimental data that we used [9], the initial temperature was 1°C below the melting temperature and the surface temperature was 40°C. The Prandtl, Stefan and Rayleigh numbers for this problem were 59.6, 0.113, and 2.63×108, respectively.
2.2. Computational procedure The numerical approach makes it possible to predict specifications of the melting process that take place inside the sphere. The flow was considered unsteady, laminar, incompressible and two-dimensional. In order to simulate the constrained melting, it is supposed that both the liquid and the solid phases are
7 Page 7 of 31
isotropic, homogeneous and remain in thermal equilibrium at the interface. For the phase-change region inside the PCM, enthalpy-porosity approach [11-12] was used. In this technique, the melt interface is not tracked explicitly. Alternatively, a quantity called the liquid fraction, indicating the fraction of the cell volume that is in liquid form, is associated with each cell in the domain. In this way, the governing equations for the PCM are as follows: Continuity: (1)
Momentum: (2)
Energy: (3)
where V is the fluid velocity vector, ρ is the density, μ is the dynamic viscosity, P is the pressure, g is the gravitational acceleration, k is the thermal conductivity and H is the specific enthalpy, which is defined as the sum of the sensible enthalpy (h) and the latent heat (ΔH):
(4)
8 Page 8 of 31
where (5)
The content of the latent heat can vary between zero (for a solid) and L (for a liquid):
(6)
The liquid fraction (β) can be written as:
(7)
By definition, solidus and liquidus temperatures are the maximum temperature at which the material is at solid state and minimum temperature at which the material is at liquid state, respectively. However, for melting of a pure substance, phase change occurs at a distinct melting temperature. Therefore, solidus and liquidus temperatures are the same. In order to determine the liquid fraction at mushy zone, according to the method was given by Voller and Prakash[11], a tiny temperature interval
near the melting point temperature
(Tm) was considered and assumed that liquefaction occurs at a certain range of temperature (
and
), as shown in Fig. 2. Consequently, using
9 Page 9 of 31
Equation (7) and by substituting
and
for
and
respectively, liquid fraction was calculated.
The source term S in the momentum equation is the Darcy's law damping term that is added to the momentum equation due to the phase change effect on the convection. That is given by: (8)
The coefficient C is a mushy zone constant. This constant measures the amplitude of the damping. The higher this value, the steeper is the transition of the velocity of the material to zero as it solidifies. Values between 10 4 and 107 are the recommended values for most calculations [12] and values exceeding this range may cause the solution to oscillate. In the present research, different values of the constant C were investigated to examine the effect of the mushy zone constant. In order to handle pressure-velocity coupling, SIMPLE (semi-implicit method for pressure-linked equations) algorithm of Patankar [13] was employed. The PRESTO scheme was used for the pressure correction equation. This scheme uses the discrete continuity balance for a “staggered” control volume about the face to compute the “staggered” pressure [13]. For flows with natural convection, rotating flows, and flows in strongly curved domain, using 10 Page 10 of 31
PRESTO scheme is recommended [12]. Other equations (i.e. momentum and energy) were discretized using power law scheme, which interpolates the value of a variable using the exact solution to a one-dimensional convection-diffusion equation. The under-relaxation factors for the velocity components, thermal energy, pressure correction and liquid fraction were 0.1, 0.3, 1 and 0.9, respectively. Different grid densities were tried to make sure that the solution is independent from the adopted grid density (see Fig. 3). As a result, 7200 grids were found sufficient for the present study. Adopting a fine grid distribution in the radial direction allows for the use of longer time steps. The time required for achieving the full melting is a good measure of time step dependence. By comparing the selected quantities obtained from simulations using 0.05, 0.1, 0. 5 and 1 s, the time step for integrating the time derivatives was set to 1 s.
3. Results and discussion The effect of different values of the mushy zone constant on the melt fraction with a Stefan number of 0.132 and a shell diameter of 101.66 mm, is shown in Fig. 4. As one can see from this figure, the value of 107 is suitable for the present simulation where the numerical result is in most accordance with the experimental data. Fig. 4 also confirms the validation of the result of CFD work.
11 Page 11 of 31
The colored contours of solid-liquid front of the PCM at various time instants for constrained melting process at a wall temperature of 40 °C and an initial temperature of 27.2 °C are shown in Fig.5. The blue and red colors represent the solid and melted parts, respectively. In the beginning of the process, the outer surface of the solid PCM is in contact with the inner wall of the sphere so that heat conduction dominates between the solid PCM and the wall. This makes the formation of a thin layer of liquid between the solid PCM and the wall of the sphere. As time progresses, the molten zone expands and the liquid layer grows. The liquid film is then forced up on the interior surface of the sphere, occupying the region above the solid PCM. By moving towards the top part of the sphere, the warm liquid is replaced with the cold liquid. Therefore, natural convection in combination with a hot rising curved wall jet are dominant at the top and side regions of the sphere. After 40 and 60 minutes, as can be seen in Fig. 5, more intense melting in the upper zone compared to the side and bottom sections is observed. This phenomenon is attributed to the strong role of natural convection in the upper zone of the sphere. After 40, 60 and 80 minutes, a waviness at the bottom part of the PCM is observed. This can be attributed to the existence of stable and unstable thermal layers at the bottom part of the sphere. This result as well as the observed shape of the solid, are in good agreement with the experimental data in the literature [9].
12 Page 12 of 31
The temperature contours of the right side of the sphere as well as the streamlines of the left side of the sphere, at a wall temperature of 40 °C for a sphere of 101.66 mm diameter, are shown in Fig. 6. At 20 and 40 minutes, formation of recirculating vortexes between the solid PCM and the inner wall of the sphere as a result of the replacement of the warm liquid with cold liquid at the bottom part of the sphere is obvious. It can also be seen that more melting occurs at the upper part of the sphere due to the contribution of natural convection (see Fig. 5). The whole PCM is melted after 120 minutes.
13 Page 13 of 31
Comparison of experimental temperatures with the results of CFD simulation for three points A, B, and C along the vertical diameter of the sphere (see Fig.1) is shown in Fig. 7.
The disordered temperature data at points A and B, show the chaotic convection motion in the unstable liquid layer at the bottom of the sphere. The computed temperatures at point C (located on the upper side of the sphere) are in a good agreement with the experimental data. This is due to the stable nature of the liquid layer at the upper part of the sphere. To prove the existence of liquid layer where the point A is located, the variation of liquid fraction at the point A during melting time was calculated (see Fig. 8). The results revealed that the PCM was at a transition state from solid to liquid during 7.5 to 9.5 minutes and after 9.5 minutes it was completely at liquid state. Therefore, the high frequency oscillation in the temperature versus time may be found at this point.
The effect of initial temperature is showed in Fig. 9. It is observed that different initial temperatures do not have a substantial effect on the melting process. The melting rates for initial temperatures of 18.2 °C and 8.2 °C are the same but the initial temperature of 27.2 °C has a few faster melting rate.
14 Page 14 of 31
Fig. 10 shows variations in the liquid fraction with time at different Stefan numbers. The Stefan number (Ste) is the ratio of sensible to latent heat of the PCM. The surface temperature (TS) directly affects the value of the Stefan number: (9)
here Cpl is the specific heat of the liquid PCM and L and Tm are the latent heat and melting temperature of the PCM, respectively. The present simulation was performed at three different Stefan numbers, namely Ste= 0.065 (Ts= 35 °C), Ste= 0.113 (Ts= 40 °C) and Ste= 0.161 (Ts=45°C). By increasing the surface temperature, the Stefan number will be increased whereas the time required for completing the melting will be decreased.
The effects of different shell diameters and Stefan numbers on the changes in the PCM liquid fraction with time are shown in Fig.11.The slopes of these graphs represent the melting rates. The highest melting rates belong to the beginning of the melting process because of the heat conduction between the inner wall of the shell and the solid PCM. As expected, for each shell diameter the liquid fraction increases more rapidly when the Stefan number (or surface temperature, Ts) is larger. Also, for each Stefan number, full melting is accomplished more rapidly in a shell with a smaller diameter. However, the effect of these parameters (i.e. diameter and surface temperature or Stefan 15 Page 15 of 31
number) are not the same on the melting time and rate. With respect to the data of Fig. 11, it can be extracted that, 28% increase in surface temperature (from to
) results in about 170% to 198% decrease in melting
time for each constant diameter. While 27% increase in diameter of container (from
to
) results in about 80% increase in melting
time. Therefore, the effect of surface temperature on melting time as well as melting rate is more than the effect of container’s diameter on these parameters.
16 Page 16 of 31
Fig.12 indicates variations in the heat flux versus time at different shell diameters and Stefan numbers. In all plots, a sharp drop in the heat flux is observed at the beginning of the melting process but after a few minutes, it decreases smoothly. When the whole PCM is melted, the heat flux almost reaches zero. The results show that by increasing the Stefan number, the heat
transfer rate increases. Integration of the area under the heat flux graphs in Fig.12 gives the total amount of heat collection. As expected, by increasing the diameter of the sphere and due to the increase in the amount of the PCM, the heat collection is enhanced.
4. Conclusions The numerical study of the constrained melting process of n-octadecane paraffin wax as a PCM in a spherical container at different operating conditions and various diameters of the container was successfully done. Detailed phase and flow fields were provided inside the system by simulations and the results were discussed and compared with the experimental data. At the beginning, the heat conduction to the PCM is dominant but after formation of a thin layer of melted liquid between the inner wall and the solid PCM, natural convection becomes dominant at the top half of the sphere. The comparison between the 17 Page 17 of 31
calculated and measured temperatures at the top half of the sphere clearly shows the stable nature of the molten liquid layer in that zone. The unstable liquid layer at the bottom of the sphere is responsible for the waviness of the PCM at this zone. Based on the simulation results, change in surface temperature has the most deterministic effect on melting rate, melting time and the rate of heat transfer in comparison with other parameters such as diameter of container and initial temperature of PCM.
Nomenclature Cpl
Specific heat of liquid PCM, Jkg-1 K-1
H
Specific enthalpy, Jkg-1
k
Thermal conductivity, Wm-1 K-1
L
Latent heat, Jkg-1
P
Pressure, Pa
18 Page 18 of 31
Ri
Internal radius, m
t
Time, s
Ti
Initial temperature, K
Tm
Melting temperature, K
Ts
Surface temperature, K
V
Velocity, ms-1
Greek symbols
Liquid fraction, dimensionless
Density, kgm-3
µ
Dynamic viscosity, kgm-1s-1
References [1] Roy SK, Sengupta S. Melting of a free solid in a spherical enclosure: effects of subcooling. J Solar Energy Eng. 1989; 111:32–6. [2] Roy SK, Sengupta S. Gravity-assisted melting in a spherical enclosure: effects of natural convection. Int J Heat Mass Transfer 1990; 33:1135–47.
19 Page 19 of 31
[3] Ettouney H, El-Dessouky H, Al-Ali A. Heat transfer during phase change of paraffin wax stored in spherical shells. ASME J Solar Energy Eng. 2005; 127: 357–65. [4] Khodadadi JM, Zhang Y. Effects of buoyancy-driven convection on melting within spherical containers. Int J Heat Mass Transfer 2001; 44:1605–18. [5] Assis E, Katsman L, Ziskind G, Letan R. Numerical and experimental study of melting in a spherical shell. Int J Heat Mass Transfer 2007; 50 (9-10): 1790– 804. [6] Veerappan M, Kalaiselvam S, Iniyan S, Goic R. Phase change characteristic study of spherical PCMs in solar energy storage. Solar Energy 2009; 83: 1245– 52. [7] Eames IW, Adref KT. Freezing and melting of water in spherical enclosures of the type used in thermal (ice) storage systems. ApplTherm Eng. 2002; 22: 733–45. [8] Regin AF, Solanki SC, Saini JS. Experimental and numerical analysis of melting of PCM inside a spherical capsule. In: Proceedings of the 9th AIAA/ASME joint thermos physics and heat transfer conference (CD ROM). Paper AIAA 2006-3618; 2006. p. 12. [9] Tan FL. Constrained and unconstrained melting inside a sphere. Int Commun Heat Mass Transfer 2008; 35:466–75.
20 Page 20 of 31
[10] Tan FL, Hosseinizadeh SF, Khodadadi JM, Fan L. Experimental and computational study of constrained melting of phase change materials (PCM) inside a spherical capsule. Int J Heat Mass Transfer 2009; 52:3464–72. [11] Voller V, Prakash C. A fixed grid numerical modeling methodology for convection-diffusion mushy region phase-change problems. Int J Heat Mass Transfer 1987; 30:1709–19. [12] Fluent 6.3.22 Users’ Guide, Fluent, Inc., 2006. [13] Patankar S. V. Numerical Heat Transfer and Fluid Flow. Hemisphere, Washington, DC. 1980.
21 Page 21 of 31
C
B A
Fig.1.The computational domain
Fig. 2. The relation between enthalpy and temperature of a pure substance.
22 Page 22 of 31
1 0.9
Liquid fraction
0.8 0.7 0.6 0.5 0.4 0.3
12800 node
0.2
5000 node
0.1
7200 node
0 0
20
40
60
80
100
120
140
Time(min)
Fig. 3. Variations of liquid fraction vs. time at different grid sizes.
Fig. 4. Variations in the liquid fraction vs. time at different mushy zone
23 Page 23 of 31
constants (Ste= 0.113 and D=101.66 mm).
0 min
20 min
40 min
60 min
80 min
100 min
120 min
24 Page 24 of 31
Fig. 5. Constrained melting phase front at wall temperature of 40 °C and initial temperature of 27.2 °C at different times.
20 min
40 min
60 min
80 min
25 Page 25 of 31
100 min
120 min
Fig. 6. Detailed temperature contours (T in Kelvin) and streamlines for a wall temperature of 40 °C at different times.
42
40
40
38 Temperature(°C)
Temperature (°c)
38 36 34 32 30 28
present work
26
36 34 32 30 28 present work
26
Experimental data[9]
24
Experimental data[9]
24 0
20
40
60 80 100 120 140 Time (min)
0
20
40
60
80
100 120 140
Time(min)
(a)
(b)
26 Page 26 of 31
40 38
Temperature(°C)
36 34 32 30 28 present work
26
Experimental data[9]
24 0
20
40
60 80 100 Time(min)
120
140
(c)
Fig. 7. Comparison of the computed and measured temperatures at (a)point A at a distance of 47.8 mm below the center (b) point B at a distance of 37.5 mm below the center (c) point C at a distance of 25 mm above the center.
Liquid fraction at point A
1 0.8 0.6 0.4 0.2 0 7
7.5
8
8.5
9
9.5
10
Time (min)
Fig. 8. Changes in liquid fraction versus time at point A at distance of 47.8 mm
27 Page 27 of 31
below the center.
1
Liquid fraction
0.8 0.6 0.4 initial temperature at 8.2 °C initial temperature at 27.2 °C
0.2
initial temperature at 18.2 °C 0 0
20
40
60 80 Time(min)
100
120
140
Fig. 9. Variations in the liquid fraction versus time at different initial temperatures.
1
Liquid fraction
0.8 0.6 0.4
ste=0.065(Ts=35°C) ste=0.113(Ts=40°C)
0.2
ste=0.161(Ts=45°C) 0 0
30
60
90
120 150 Time(min)
180
210
240
28 Page 28 of 31
Fig.10. Variations in the liquid fraction versus time at different surface temperatures.
0.8
0.8 Liquid fraction
1
Liquid fraction
1
0.6 0.4 Ste=0.065
0.2
Ste=0.113
0.6 0.4 Ste=0.065 Ste=0.113
0.2
Ste=0.161
Ste=0.161 0
0 0
20
40
60
80
100 120 140
0
20 40 60 80 100 120 140 160 180 200
Time(min)
Time(min)
(a) D=60 mm
(b) D=80 mm
1
Liquid fraction
0.8 0.6 0.4 Ste=0.065 Ste=0.113
0.2
Ste=0.161 0 0
30
60
90 120 150 180 210 240 Time(min)
(c) D=101.66 mm
Fig.11. Changes in the liquid fraction versus time at different Stefan numbers
29 Page 29 of 31
8000
6000
7000
5000
6000
Heat flux(Wm-2)
7000
4000 3000 ste=0.065 ste=0.113 ste=0.161
2000 1000 0
50 100 Time(min)
5000 4000 3000 ste=0.065 ste=0.113 ste=0.161
2000 1000
0
0
150
0
50
(a) D=60 mm
100 Time(min)
150
200
(b) D=80 mm
7000 6000 Heat flux(Wm-2)
Heat flux(Wm-2)
and different shell diameters: (a) D=60 mm (b) D=80 mm (c) D=101. 66 mm.
5000 4000 3000 ste=0.065
2000
ste=0.113
1000
ste=0.161
0 0
50
100 150 Time(min)
200
250
(c) D=101.66 mm
Fig.12. Variations in the heat flux at different Stefan numbers and shell diameters: (a) 60 mm (b) 80 mm (c) 101.66 mm.
30 Page 30 of 31
Table 1. Thermophysical properties of n-octadecane. Melting temperature
28.2 oC
Density
772 kgm-3
Kinematic viscosity
5 x 10-6 m2s-1
Specific heat
2330 Jkg-1K-1
Thermal conductivity
0.1505 Wm-1K-1
Latent heat of fusion
243.5 kJkg-1
Thermal expansion
0.00091 K-1
coefficient
31 Page 31 of 31
S0140-7007(16)30287-0 http://dx.doi.org/doi: 10.1016/j.ijrefrig.2016.09.007 JIJR 3421
To appear in:
International Journal of Refrigeration
Received date: Revised date: Accepted date:
15-2-2016 16-7-2016 5-9-2016
Please cite this article as: H. Sattari, A. Mohebbi, M.M. Afsahi, A. Azimi Yancheshme, CFD simulation of melting process of phase change materials (PCMs) in a spherical capsule, International Journal of Refrigeration (2016), http://dx.doi.org/doi: 10.1016/j.ijrefrig.2016.09.007. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
CFD simulation of melting process of phase change materials (PCMs) in a spherical capsule H. Sattaria, A. Mohebbia,*, M. M. Afsahia, A. Azimi Yancheshmeb a
Department of Chemical Engineering, Faculty of Engineering, Shahid Bahonar University of Kerman, Kerman, Iran b
Chemical and Petroleum Engineering Department, Sharif University of Technology, Tehran, Iran *Corresponding author: Tel & Fax: +983432118298 E-mail: [email protected], [email protected]
Highlights The effect of container’s surface temperature on melting rate is more than other parameters. At the beginning of the melting process of PCM, the conduction mode is governing. Natural convection is the ruling mode at the top half of container. The molten liquid at the top half of the container is naturally stable. Liquid layer, which is formed at the bottom of the container is unstable.
Abstract The present study is focused on CFD simulation of constrained melting of Phase Change Materials (PCMs) in a spherical container. To investigate the melting process of the PCM, its melting fraction was analyzed at different 1 Page 1 of 31
times. The results indicated the existence of thermally stable layers on the top of the sphere. Moreover, inspection of the calculated temperatures at different points along the vertical axis indicates the existence of some disturbances at the bottom of the sphere due to the natural convection. After the validation of the results, the effects of different parameters such as the surface temperature of the capsule, the initial temperature and the size of the spherical capsules, on the melting process were investigated. The initial temperature did not affect the melting rate, whereas melting rate increased by increasing the surface temperature of the capsule and also decreasing the diameter of the sphere. The results showed that the surface temperature of capsule compared to geometrical parameters and other operational conditions can have a greater influence on the melting rate and the heat flux. Keywords: Constrained melting; Phase change materials; Computational fluid dynamics (CFD); Spherical container
1. Introduction Recently, latent heat thermal energy storage systems have been widely applied in cool storage for central air-conditioning, hot storage for air heating, solar energy, energy efficient buildings and waste heat recovery. Storage 2 Page 2 of 31
systems allow for the efficient and rational utilization of the available resources or renewable energies, by using the time lag between production, as well as the availability of the energy and its consumption in the demanding systems. Phase change materials (PCM) are mainly used to provide higher storage densities. As a serious problem PCMs have relatively low thermal conductivities. Therefore, increasing the surface/volume ratio of PCMs in order to increase the heattransfer rate is very appealing. This can be done by packing a volume with a great number of PCM capsules. The spherical geometry is one interesting case for heat storage applications, since spheres are much employed in packed beds. Due to the complexity of such systems, it is often more effective to first model the behavior of an individual sphere and then describe it with a simple parametric model in the packed bed modeling. Roy and Sengupta [1] examined the melting process with the solid phase initially uniformly super cooled. In order to include the effects of a temperature gradient in the solid core, they modified the heat transfer equation. At each time step, the unsteady conduction equation has been solved numerically using a toroidal coordinate system with a suitable immobilization of the moving boundary to transform the infinite domain into a finite one. In another paper [2], they investigated the outcome of natural convection on the melting process. In order to reduce the computational effort and the time, they made some simplifications. They obtained the non-dimensional melting time and the heat
3 Page 3 of 31
transfer coefficient as functions of the PCM property values, the operating temperature and the physical size. Ettouney et al. [3] experimentally evaluated the heat transfer (during energy storage) and release for the phase change of paraffin wax in the spherical shells. They showed that an increment in the Nusselt number of the sphere with a larger diameter is attributed to the increase in the natural convection cells in the PCM inside the sphere. The natural convection role is enhanced upon an increase in the sphere diameter and the air temperature. On the other hand, during solidification the wax layers are formed inward, nucleating on the sphere walls. As solidification progresses, the melt volume becomes smaller and the role of natural convection diminishes rapidly. Khodadadi and Zhang [4] also considered the effects of buoyancy-driven convection on the constrained melting of phase change materials within spherical containers. They reported that during the initial melting process, the conduction mode of the heat transfer is dominant. As the buoyancy-driven convection is strengthened due to the increase in the melting zone, the process at top region of the sphere is much faster than at the bottom due to the increment of the conduction mode of heat transfer. They found that buoyancydriven convection speeds up the melting process compared to the diffusioncontrolled melting. Assis et al. [5] numerically and experimentally explored the melting process of a PCM in a spherical geometry. The results of their experimental 4 Page 4 of 31
investigation, including visualization, were favorably comparable with the numerical results and therefore suitable for validation of the mathematical approach. Their computational results showed how the transient phase-change process depends on the thermal and geometrical parameters of the arrangement. Veerappan et al. [6] investigated the phase change behavior of 65 mol% capric acid and 35%mollauric acid, calcium chloride hexa hydrate, n-octadecane, nhexadecane, and n-eicosane inside the spherical containers to identify a suitable heat storage material. They created analytical models for solidification and melting of PCM in a spherical shell with conduction, natural convection, and heat generation and found a good agreement between the analytical predictions and the experimental data. Both models were validated with the experimental work of Eames and Adref [7]. Regin et al. [8] examined the heat transfer performance of a spherical capsule using paraffin wax as PCM placed in a convective environment during the melting process. The model results were in a good agreement with the experimental data. Tan [9] investigated the melting of the phase change material (PCM) inside a sphere using n-octadecane for both constrained and unconstrained melting processes. For constrained melting, paraffin wax (n-octadecane) was immobilized, through the use of thermocouples when melting is done inside a transparent glass sphere. Their experimental setup provided a detailed temperature data that were gathered along the vertical diameter of the sphere 5 Page 5 of 31
during the melting process. Tan and Khodadadi [10] also experimentally and numerically investigated the constrained melting of PCMs inside a spherical capsule to understand the role of the buoyancy-driven convection. In this study, a computational fluid dynamics (CFD) modeling on the constrained melting of phase change material (PCM) in a spherical container was performed. The effects of different parameters like the capsule size, the surface temperature, the initial temperature and the Stefan number, on the PCM melting process were investigated.
2. CFD simulation CFD simulation was carried out using the commercial Fluent 6.3 software. The physical model and computational procedure are discussed below.
2.1. Physical model Consider a spherical shell with an internal radius (Ri) of 50.83 mm and a wall thickness of 1.5 mm initially filled with a solid PCM at an initial temperature (Ti) lower than the melting temperature (Tm).The schematic diagram of the physical model is shown in Fig. 1. At t (time) > 0, a constant temperature (T0) greater than the melting temperature is imposed on the surface of the sphere, i.e. T0>Tm. Melting starts at the surface, moving the solid-liquid interface towards the middle of the field. Since the present work is focused on
6 Page 6 of 31
the analysis of constrained melting, both the solid and liquid phases have the same density. The container was made of glass with a thermal conductivity of 1.14 Wm-1K-1. The properties of the PCM, based on n-octadecane as a commercially available material, are given in Table 1.
For the experimental data that we used [9], the initial temperature was 1°C below the melting temperature and the surface temperature was 40°C. The Prandtl, Stefan and Rayleigh numbers for this problem were 59.6, 0.113, and 2.63×108, respectively.
2.2. Computational procedure The numerical approach makes it possible to predict specifications of the melting process that take place inside the sphere. The flow was considered unsteady, laminar, incompressible and two-dimensional. In order to simulate the constrained melting, it is supposed that both the liquid and the solid phases are
7 Page 7 of 31
isotropic, homogeneous and remain in thermal equilibrium at the interface. For the phase-change region inside the PCM, enthalpy-porosity approach [11-12] was used. In this technique, the melt interface is not tracked explicitly. Alternatively, a quantity called the liquid fraction, indicating the fraction of the cell volume that is in liquid form, is associated with each cell in the domain. In this way, the governing equations for the PCM are as follows: Continuity: (1)
Momentum: (2)
Energy: (3)
where V is the fluid velocity vector, ρ is the density, μ is the dynamic viscosity, P is the pressure, g is the gravitational acceleration, k is the thermal conductivity and H is the specific enthalpy, which is defined as the sum of the sensible enthalpy (h) and the latent heat (ΔH):
(4)
8 Page 8 of 31
where (5)
The content of the latent heat can vary between zero (for a solid) and L (for a liquid):
(6)
The liquid fraction (β) can be written as:
(7)
By definition, solidus and liquidus temperatures are the maximum temperature at which the material is at solid state and minimum temperature at which the material is at liquid state, respectively. However, for melting of a pure substance, phase change occurs at a distinct melting temperature. Therefore, solidus and liquidus temperatures are the same. In order to determine the liquid fraction at mushy zone, according to the method was given by Voller and Prakash[11], a tiny temperature interval
near the melting point temperature
(Tm) was considered and assumed that liquefaction occurs at a certain range of temperature (
and
), as shown in Fig. 2. Consequently, using
9 Page 9 of 31
Equation (7) and by substituting
and
for
and
respectively, liquid fraction was calculated.
The source term S in the momentum equation is the Darcy's law damping term that is added to the momentum equation due to the phase change effect on the convection. That is given by: (8)
The coefficient C is a mushy zone constant. This constant measures the amplitude of the damping. The higher this value, the steeper is the transition of the velocity of the material to zero as it solidifies. Values between 10 4 and 107 are the recommended values for most calculations [12] and values exceeding this range may cause the solution to oscillate. In the present research, different values of the constant C were investigated to examine the effect of the mushy zone constant. In order to handle pressure-velocity coupling, SIMPLE (semi-implicit method for pressure-linked equations) algorithm of Patankar [13] was employed. The PRESTO scheme was used for the pressure correction equation. This scheme uses the discrete continuity balance for a “staggered” control volume about the face to compute the “staggered” pressure [13]. For flows with natural convection, rotating flows, and flows in strongly curved domain, using 10 Page 10 of 31
PRESTO scheme is recommended [12]. Other equations (i.e. momentum and energy) were discretized using power law scheme, which interpolates the value of a variable using the exact solution to a one-dimensional convection-diffusion equation. The under-relaxation factors for the velocity components, thermal energy, pressure correction and liquid fraction were 0.1, 0.3, 1 and 0.9, respectively. Different grid densities were tried to make sure that the solution is independent from the adopted grid density (see Fig. 3). As a result, 7200 grids were found sufficient for the present study. Adopting a fine grid distribution in the radial direction allows for the use of longer time steps. The time required for achieving the full melting is a good measure of time step dependence. By comparing the selected quantities obtained from simulations using 0.05, 0.1, 0. 5 and 1 s, the time step for integrating the time derivatives was set to 1 s.
3. Results and discussion The effect of different values of the mushy zone constant on the melt fraction with a Stefan number of 0.132 and a shell diameter of 101.66 mm, is shown in Fig. 4. As one can see from this figure, the value of 107 is suitable for the present simulation where the numerical result is in most accordance with the experimental data. Fig. 4 also confirms the validation of the result of CFD work.
11 Page 11 of 31
The colored contours of solid-liquid front of the PCM at various time instants for constrained melting process at a wall temperature of 40 °C and an initial temperature of 27.2 °C are shown in Fig.5. The blue and red colors represent the solid and melted parts, respectively. In the beginning of the process, the outer surface of the solid PCM is in contact with the inner wall of the sphere so that heat conduction dominates between the solid PCM and the wall. This makes the formation of a thin layer of liquid between the solid PCM and the wall of the sphere. As time progresses, the molten zone expands and the liquid layer grows. The liquid film is then forced up on the interior surface of the sphere, occupying the region above the solid PCM. By moving towards the top part of the sphere, the warm liquid is replaced with the cold liquid. Therefore, natural convection in combination with a hot rising curved wall jet are dominant at the top and side regions of the sphere. After 40 and 60 minutes, as can be seen in Fig. 5, more intense melting in the upper zone compared to the side and bottom sections is observed. This phenomenon is attributed to the strong role of natural convection in the upper zone of the sphere. After 40, 60 and 80 minutes, a waviness at the bottom part of the PCM is observed. This can be attributed to the existence of stable and unstable thermal layers at the bottom part of the sphere. This result as well as the observed shape of the solid, are in good agreement with the experimental data in the literature [9].
12 Page 12 of 31
The temperature contours of the right side of the sphere as well as the streamlines of the left side of the sphere, at a wall temperature of 40 °C for a sphere of 101.66 mm diameter, are shown in Fig. 6. At 20 and 40 minutes, formation of recirculating vortexes between the solid PCM and the inner wall of the sphere as a result of the replacement of the warm liquid with cold liquid at the bottom part of the sphere is obvious. It can also be seen that more melting occurs at the upper part of the sphere due to the contribution of natural convection (see Fig. 5). The whole PCM is melted after 120 minutes.
13 Page 13 of 31
Comparison of experimental temperatures with the results of CFD simulation for three points A, B, and C along the vertical diameter of the sphere (see Fig.1) is shown in Fig. 7.
The disordered temperature data at points A and B, show the chaotic convection motion in the unstable liquid layer at the bottom of the sphere. The computed temperatures at point C (located on the upper side of the sphere) are in a good agreement with the experimental data. This is due to the stable nature of the liquid layer at the upper part of the sphere. To prove the existence of liquid layer where the point A is located, the variation of liquid fraction at the point A during melting time was calculated (see Fig. 8). The results revealed that the PCM was at a transition state from solid to liquid during 7.5 to 9.5 minutes and after 9.5 minutes it was completely at liquid state. Therefore, the high frequency oscillation in the temperature versus time may be found at this point.
The effect of initial temperature is showed in Fig. 9. It is observed that different initial temperatures do not have a substantial effect on the melting process. The melting rates for initial temperatures of 18.2 °C and 8.2 °C are the same but the initial temperature of 27.2 °C has a few faster melting rate.
14 Page 14 of 31
Fig. 10 shows variations in the liquid fraction with time at different Stefan numbers. The Stefan number (Ste) is the ratio of sensible to latent heat of the PCM. The surface temperature (TS) directly affects the value of the Stefan number: (9)
here Cpl is the specific heat of the liquid PCM and L and Tm are the latent heat and melting temperature of the PCM, respectively. The present simulation was performed at three different Stefan numbers, namely Ste= 0.065 (Ts= 35 °C), Ste= 0.113 (Ts= 40 °C) and Ste= 0.161 (Ts=45°C). By increasing the surface temperature, the Stefan number will be increased whereas the time required for completing the melting will be decreased.
The effects of different shell diameters and Stefan numbers on the changes in the PCM liquid fraction with time are shown in Fig.11.The slopes of these graphs represent the melting rates. The highest melting rates belong to the beginning of the melting process because of the heat conduction between the inner wall of the shell and the solid PCM. As expected, for each shell diameter the liquid fraction increases more rapidly when the Stefan number (or surface temperature, Ts) is larger. Also, for each Stefan number, full melting is accomplished more rapidly in a shell with a smaller diameter. However, the effect of these parameters (i.e. diameter and surface temperature or Stefan 15 Page 15 of 31
number) are not the same on the melting time and rate. With respect to the data of Fig. 11, it can be extracted that, 28% increase in surface temperature (from to
) results in about 170% to 198% decrease in melting
time for each constant diameter. While 27% increase in diameter of container (from
to
) results in about 80% increase in melting
time. Therefore, the effect of surface temperature on melting time as well as melting rate is more than the effect of container’s diameter on these parameters.
16 Page 16 of 31
Fig.12 indicates variations in the heat flux versus time at different shell diameters and Stefan numbers. In all plots, a sharp drop in the heat flux is observed at the beginning of the melting process but after a few minutes, it decreases smoothly. When the whole PCM is melted, the heat flux almost reaches zero. The results show that by increasing the Stefan number, the heat
transfer rate increases. Integration of the area under the heat flux graphs in Fig.12 gives the total amount of heat collection. As expected, by increasing the diameter of the sphere and due to the increase in the amount of the PCM, the heat collection is enhanced.
4. Conclusions The numerical study of the constrained melting process of n-octadecane paraffin wax as a PCM in a spherical container at different operating conditions and various diameters of the container was successfully done. Detailed phase and flow fields were provided inside the system by simulations and the results were discussed and compared with the experimental data. At the beginning, the heat conduction to the PCM is dominant but after formation of a thin layer of melted liquid between the inner wall and the solid PCM, natural convection becomes dominant at the top half of the sphere. The comparison between the 17 Page 17 of 31
calculated and measured temperatures at the top half of the sphere clearly shows the stable nature of the molten liquid layer in that zone. The unstable liquid layer at the bottom of the sphere is responsible for the waviness of the PCM at this zone. Based on the simulation results, change in surface temperature has the most deterministic effect on melting rate, melting time and the rate of heat transfer in comparison with other parameters such as diameter of container and initial temperature of PCM.
Nomenclature Cpl
Specific heat of liquid PCM, Jkg-1 K-1
H
Specific enthalpy, Jkg-1
k
Thermal conductivity, Wm-1 K-1
L
Latent heat, Jkg-1
P
Pressure, Pa
18 Page 18 of 31
Ri
Internal radius, m
t
Time, s
Ti
Initial temperature, K
Tm
Melting temperature, K
Ts
Surface temperature, K
V
Velocity, ms-1
Greek symbols
Liquid fraction, dimensionless
Density, kgm-3
µ
Dynamic viscosity, kgm-1s-1
References [1] Roy SK, Sengupta S. Melting of a free solid in a spherical enclosure: effects of subcooling. J Solar Energy Eng. 1989; 111:32–6. [2] Roy SK, Sengupta S. Gravity-assisted melting in a spherical enclosure: effects of natural convection. Int J Heat Mass Transfer 1990; 33:1135–47.
19 Page 19 of 31
[3] Ettouney H, El-Dessouky H, Al-Ali A. Heat transfer during phase change of paraffin wax stored in spherical shells. ASME J Solar Energy Eng. 2005; 127: 357–65. [4] Khodadadi JM, Zhang Y. Effects of buoyancy-driven convection on melting within spherical containers. Int J Heat Mass Transfer 2001; 44:1605–18. [5] Assis E, Katsman L, Ziskind G, Letan R. Numerical and experimental study of melting in a spherical shell. Int J Heat Mass Transfer 2007; 50 (9-10): 1790– 804. [6] Veerappan M, Kalaiselvam S, Iniyan S, Goic R. Phase change characteristic study of spherical PCMs in solar energy storage. Solar Energy 2009; 83: 1245– 52. [7] Eames IW, Adref KT. Freezing and melting of water in spherical enclosures of the type used in thermal (ice) storage systems. ApplTherm Eng. 2002; 22: 733–45. [8] Regin AF, Solanki SC, Saini JS. Experimental and numerical analysis of melting of PCM inside a spherical capsule. In: Proceedings of the 9th AIAA/ASME joint thermos physics and heat transfer conference (CD ROM). Paper AIAA 2006-3618; 2006. p. 12. [9] Tan FL. Constrained and unconstrained melting inside a sphere. Int Commun Heat Mass Transfer 2008; 35:466–75.
20 Page 20 of 31
[10] Tan FL, Hosseinizadeh SF, Khodadadi JM, Fan L. Experimental and computational study of constrained melting of phase change materials (PCM) inside a spherical capsule. Int J Heat Mass Transfer 2009; 52:3464–72. [11] Voller V, Prakash C. A fixed grid numerical modeling methodology for convection-diffusion mushy region phase-change problems. Int J Heat Mass Transfer 1987; 30:1709–19. [12] Fluent 6.3.22 Users’ Guide, Fluent, Inc., 2006. [13] Patankar S. V. Numerical Heat Transfer and Fluid Flow. Hemisphere, Washington, DC. 1980.
21 Page 21 of 31
C
B A
Fig.1.The computational domain
Fig. 2. The relation between enthalpy and temperature of a pure substance.
22 Page 22 of 31
1 0.9
Liquid fraction
0.8 0.7 0.6 0.5 0.4 0.3
12800 node
0.2
5000 node
0.1
7200 node
0 0
20
40
60
80
100
120
140
Time(min)
Fig. 3. Variations of liquid fraction vs. time at different grid sizes.
Fig. 4. Variations in the liquid fraction vs. time at different mushy zone
23 Page 23 of 31
constants (Ste= 0.113 and D=101.66 mm).
0 min
20 min
40 min
60 min
80 min
100 min
120 min
24 Page 24 of 31
Fig. 5. Constrained melting phase front at wall temperature of 40 °C and initial temperature of 27.2 °C at different times.
20 min
40 min
60 min
80 min
25 Page 25 of 31
100 min
120 min
Fig. 6. Detailed temperature contours (T in Kelvin) and streamlines for a wall temperature of 40 °C at different times.
42
40
40
38 Temperature(°C)
Temperature (°c)
38 36 34 32 30 28
present work
26
36 34 32 30 28 present work
26
Experimental data[9]
24
Experimental data[9]
24 0
20
40
60 80 100 120 140 Time (min)
0
20
40
60
80
100 120 140
Time(min)
(a)
(b)
26 Page 26 of 31
40 38
Temperature(°C)
36 34 32 30 28 present work
26
Experimental data[9]
24 0
20
40
60 80 100 Time(min)
120
140
(c)
Fig. 7. Comparison of the computed and measured temperatures at (a)point A at a distance of 47.8 mm below the center (b) point B at a distance of 37.5 mm below the center (c) point C at a distance of 25 mm above the center.
Liquid fraction at point A
1 0.8 0.6 0.4 0.2 0 7
7.5
8
8.5
9
9.5
10
Time (min)
Fig. 8. Changes in liquid fraction versus time at point A at distance of 47.8 mm
27 Page 27 of 31
below the center.
1
Liquid fraction
0.8 0.6 0.4 initial temperature at 8.2 °C initial temperature at 27.2 °C
0.2
initial temperature at 18.2 °C 0 0
20
40
60 80 Time(min)
100
120
140
Fig. 9. Variations in the liquid fraction versus time at different initial temperatures.
1
Liquid fraction
0.8 0.6 0.4
ste=0.065(Ts=35°C) ste=0.113(Ts=40°C)
0.2
ste=0.161(Ts=45°C) 0 0
30
60
90
120 150 Time(min)
180
210
240
28 Page 28 of 31
Fig.10. Variations in the liquid fraction versus time at different surface temperatures.
0.8
0.8 Liquid fraction
1
Liquid fraction
1
0.6 0.4 Ste=0.065
0.2
Ste=0.113
0.6 0.4 Ste=0.065 Ste=0.113
0.2
Ste=0.161
Ste=0.161 0
0 0
20
40
60
80
100 120 140
0
20 40 60 80 100 120 140 160 180 200
Time(min)
Time(min)
(a) D=60 mm
(b) D=80 mm
1
Liquid fraction
0.8 0.6 0.4 Ste=0.065 Ste=0.113
0.2
Ste=0.161 0 0
30
60
90 120 150 180 210 240 Time(min)
(c) D=101.66 mm
Fig.11. Changes in the liquid fraction versus time at different Stefan numbers
29 Page 29 of 31
8000
6000
7000
5000
6000
Heat flux(Wm-2)
7000
4000 3000 ste=0.065 ste=0.113 ste=0.161
2000 1000 0
50 100 Time(min)
5000 4000 3000 ste=0.065 ste=0.113 ste=0.161
2000 1000
0
0
150
0
50
(a) D=60 mm
100 Time(min)
150
200
(b) D=80 mm
7000 6000 Heat flux(Wm-2)
Heat flux(Wm-2)
and different shell diameters: (a) D=60 mm (b) D=80 mm (c) D=101. 66 mm.
5000 4000 3000 ste=0.065
2000
ste=0.113
1000
ste=0.161
0 0
50
100 150 Time(min)
200
250
(c) D=101.66 mm
Fig.12. Variations in the heat flux at different Stefan numbers and shell diameters: (a) 60 mm (b) 80 mm (c) 101.66 mm.
30 Page 30 of 31
Table 1. Thermophysical properties of n-octadecane. Melting temperature
28.2 oC
Density
772 kgm-3
Kinematic viscosity
5 x 10-6 m2s-1
Specific heat
2330 Jkg-1K-1
Thermal conductivity
0.1505 Wm-1K-1
Latent heat of fusion
243.5 kJkg-1
Thermal expansion
0.00091 K-1
coefficient
31 Page 31 of 31