Natural convection in melting and solidification

From Thermal-FluidsPedia

Interest in utilizing clean energy is growing due to environmental considerations. Solar energy is one of the most promising clean energies, but due to its periodic nature, a heat storage device is needed. A thermal energy storage system using a Phase Change Material (PCM) is an attractive method because a large amount of latent heat can be stored or released during the melting or solidification process of the PCM. During the thermal energy storage process, the hot working fluid from the solar energy collector flows into the thermal energy storage device to melt the PCM and store thermal energy in the form of latent heat. During the nighttime, the PCM solidifies and the latent heat is released. Natural convection plays a significant role in both melting and solidification processes due to the temperature difference in the liquid phase. While melting is promoted by the presence of natural convection, solidification is slowed down by its presence.

Natural convection during solidification around a horizontal tube and during melting in an enclosure heated from the side have been investigated by (Zhang et al., 1997) and (Zhang, 2002), respectively.

Solidification around Horizontal Cylinder

Figure 1 Solidification around a horizontal tube.

The theoretical model to be employed in this study is shown in Fig. 1. At the very beginning of the process (t = 0), the tube, which has a radius of R0, is surrounded by a liquid phase change material with uniform temperature T_f > T_m. The temperature of the working fluid inside the tube is T_i, and the convective heat transfer coefficient between the working fluid and the internal tube wall is h_i. Both h_i and T_i are kept constant throughout the process. The thickness of the tube is assumed to be very thin, so the thermal resistance of the tube wall can be neglected. The phase change material can be treated as if it were directly in contact with the coolant inside the tube. The liquid adjacent to the cooled surface will be solidified, and the temperature difference between the solid-liquid interface and the otherwise quiescent liquid will drive natural convection in the liquid region. The liquid region is assumed to be sufficiently large so that the convection can be treated as natural convection between the surface and an extensive fluid medium. Conduction is the only heat transfer mode in the solidified region. The analysis is based on several additional assumptions:

1. The liquid is Newtonian and Boussinesq as well as incompressible. The Prandtl number of the liquid phase change material is greater than unity.

2. The solid is homogeneous and isotropic.

3. The liquid motion induced by volumetric variation during solidification is neglected, i.e., the density of the solid is equal to the density of the liquid. In addition, the properties of the phase change material are constant in the liquid and solid regions.

4. The solid-liquid interface is assumed to be a smooth cylinder concentric with the cooled tube.

5. The effects of natural convection are restricted to within the boundary layer and the bulk liquid has a uniform temperature, T_f.

Compared to the thickness of the solidified layer, the thickness of the natural convective boundary layer on the solid-liquid interface is very thin, except at the onset of freezing (the boundary layer thickness in Fig 1 is significantly exaggerated so as to be clearly visible). However, since the onset of freezing is a very short period compared to the whole solidification process, it is reasonable to neglect the curvature effect in the equations of the boundary layer. The boundary layer equations of the problem can be written as follows:

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$	(1)

$\nu \frac{{{\partial }^{2}}u}{\partial {{y}^{2}}}+g\beta (T-{{T}_{m}})\sin \theta =0$	(2)

$\frac{\partial T}{\partial t}+\frac{\partial (uT)}{\partial x}+\frac{\partial (vT)}{\partial y}={{\alpha }_{f}}\frac{{{\partial }^{2}}T}{\partial {{y}^{2}}}$	(3)

Since the Prandtl number of the liquid is much larger than unity, the inertia terms in the momentum equation have been neglected. These equations were solved by an integral method. Assuming a polynomial profile, the temperature and velocity profiles inside the boundary layer of thickness $δ$ are expressed as

$T={{T}_{f}}-({{T}_{f}}-{{T}_{m}}){{\left( 1-\frac{y}{\delta } \right)}^{2}}$	(4)

$u=U\frac{y}{\delta }{{\left( 1-\frac{y}{\delta } \right)}^{2}}$	(4)

where U is a characteristic velocity, $δ$ is the boundary layer thickness, and both of them are functions of time t and angle θ. Integrating both sides of eq. (2) with respect to y within the interval (0, $δ$ ), the expression for U in eq. (4) should be

$U=\frac{g\beta \sin \theta }{3\nu }({{T}_{f}}-{{T}_{m}}){{\delta }^{2}}$	(5)

Integrating eq. (3) in the same manner, one obtains

$-\frac{1}{3}\frac{\partial \delta }{\partial t}+\frac{dR}{dt}-\frac{1}{90}({{T}_{f}}-{{T}_{m}})\frac{g\beta \sin \theta }{\nu R}\frac{{{\partial }^{3}}\delta }{\partial {{\theta }^{3}}}+\frac{2{{\alpha }_{f}}}{\delta }=0$	(6)

Heat transfer in the solidified layer is dominated by conduction. The governing equation of the solid layer and the corresponding boundary conditions are as follows:

$\frac{1}{r}\frac{\partial }{\partial r}\left( r\frac{\partial T}{\partial r} \right)=\frac{1}{{{\alpha }_{s}}}\frac{\partial T}{\partial t}\quad {{R}_{0}}<r<R\quad t>0$	(7)

${{k}_{s}}\frac{\partial T}{\partial r}={{h}_{i}}(T-{{T}_{i}})\quad r={{R}_{0}}\quad t>0$	(8)

${{T}_{s}}(r,t)={{T}_{m}}\quad r=R\quad t>0$	(9)

${{\left. {{k}_{s}}\frac{\partial {{T}_{s}}}{\partial r} \right\|}_{r=R}}-{{\left. {{k}_{f}}\frac{\partial T}{\partial y} \right\|}_{y=0}}=\rho {{h}_{s\ell }}\frac{dR}{dt}\quad r=R\quad t>0$	(10)

Equations (6) – (10) can be nondimensionalized and solved using an integral approximation method. When $\text{Bi}\to \infty$ , it corresponds to boundary conditions of the first kind, i.e. the tube wall temperature is $T i$ and is kept steady throughout the process. Wang et al. (1991) experimentally investigated the solidification process around a horizontal cooled tube. Zhang et al., 1997 compared the predicted solidification rate, $V/{{V}_{0}}={{(R/{{R}_{0}})}^{2}}$ with the experimental results. When Ra = 0, i.e., no superheat exists in the liquid region or the solidification process is dominated by conduction, the predicted value is 18% lower than the experimental data. During the conduction-dominated freezing process, the front of the freezing layer is a dendritic layer. Therefore, it is believed that the solid-liquid interface is extended by the dendritic layer. As the Rayleigh number increases, natural convection occurs in the liquid region, and the solid-liquid interface becomes smooth because of the natural motion of the liquid. When $\text{Ra}=1.8\times {{10}^{5}}$ , the predicted value is only 8% lower than the experimental data, so the agreement is satisfactory.

Melting in a Rectangular Enclosure Heated from the Side

Figure 2 Melting in an enclosure heated from the side.

Natural convection controlled melting in a rectangular cavity under boundary condition of the first and second kinds has been investigated extensively (Ho and Viskanta, 1984; Gadgil and Gobin, 1984; Bejan, 1989; Cao and Faghri, 1990; Zhang and Bejan, 1989; Zhang and Chen, 1994; Wang et al., 2010). In a solar energy utilization system, the working fluid absorbs heat in the solar energy collector and transfers this heat to the latent heat thermal energy storage system via convection. Therefore, melting in the latent heat thermal energy storage system is under boundary condition of neither the first kind nor the second kind. It is instead under boundary condition of the third kind. In order to thoroughly understand the mechanism of the latent heat thermal energy storage unit, it is necessary to study melting in an enclosure under the boundary condition of the third kind.

The physical model under consideration is shown in Fig. 2. The enclosure with height of H and width of L is filled with a solid PCM at its melting point, Tm. The top, bottom and the left side of the enclosure are adiabatic. At time t > 0, the right side of the enclosure is exposed to the working fluid at temperature of Te. The convective heat transfer coefficient between the working fluid and the right wall is constantly equal to he. It is assumed the thickness of the right-side wall is very small and therefore, the conduction thermal resistance of the wall is negligible. The melting process can be divided into two stages: conduction and convection stages. During the conduction stage, conduction is the dominant mode of heat transfer. The problem can be simplified as a 1-D melting problem under boundary condition of the third kind and the analytical solution is available. The melting process enters the convection stage when the thickness of the melting layer increases and natural convection becomes the dominant mode of heat transfer. Although the process is still unsteady, the transient terms in the governing equations describing the melting process have negligible effect and the process can be considered to be quasi-steady (Bejan, 1989).

When the melting is in the early stage of the convection regime, the solid-liquid interface can be assumed to be sufficiently straight and vertical. In addition, the following assumptions are made:

1. The properties of liquid PCM are constants except for density, which is a function of temperature (Boussinesq assumption). The liquid is Newtonian and incompressible.

2. The Prandtl number of the liquid phase change material is much greater than 1. This assumption is valid for most PCMs (Zhang and Bejan, 1989).

3. The liquid motion induced by volumetric variation during melting is neglected, i.e. the density of the solid is equal to that of the liquid.

The liquid phase can be divided into three regions: (1) a cold boundary layer near the solid-liquid interface, (2) a warm boundary layer near the heated vertical wall, and (3) the core region between these two boundary layers (Bejan, 1989).

The analysis of the cold boundary layer is performed in a frame (x-y system) that is attached to the solid-liquid interface. Since the Prandtl number of the liquid is much larger than unity, the inertia terms in the momentum equation can be neglected. The advection induced by the solid- liquid interface motion in the boundary layer equations can also be neglected since the solid-liquid interface velocity is typically several orders of magnitude smaller than the natural convective velocity in the boundary layer (Gadgil and Gobin 1984). The momentum and energy equations for the cold boundary layer are:

$\frac{{{\partial }^{3}}v}{\partial {{x}^{3}}}+\frac{g\beta }{\nu }\frac{\partial T}{\partial x}=0$	(11)

$u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\alpha \frac{{{\partial }^{2}}T}{\partial {{x}^{2}}}$	(12)

which can be solved by an integral method. By following the procedure similar to Bejan (1989) and Zhang and Bejan (1989), one obtains an ordinary differential equation for cold boundary layer thickness, δ:

$\frac{2k{{({{T}_{c}}-{{T}_{m}})}^{2}}}{\rho {{h}_{s\ell }}\delta }+\frac{g\beta ({{T}_{c}}-{{T}_{m}})}{36\nu }\frac{d}{dy}\left[ {{\delta }^{3}}({{T}_{c}}-{{T}_{m}}) \right]-\frac{g\beta }{60\nu }\frac{d}{dy}\left[ {{\delta }^{3}}{{({{T}_{c}}-{{T}_{m}})}^{2}} \right]=-2\alpha \frac{{{T}_{c}}-{{T}_{m}}}{\delta }$	(13)

where T_c is the core temperature. Since the cold boundary layer starts from the top of the enclosure, the boundary condition of eq. (13) is

$δ(H) = 0$	(14)

The solid-liquid interface velocity, u₀, can be obtained by energy balance at the interface. Considering the assumed temperature profile in the cold boundary layer used in the integral solution, the solid-liquid interface velocity is

${{u}_{0}}=\frac{2k({{T}_{c}}-{{T}_{m}})}{\rho {{h}_{s\ell }}\delta }$	(15)

The warm boundary layer is attached to the heated wall, i.e. on to a stationary frame (xr-y). The analysis of the warm boundary layer is very similar to that of the cold boundary layer. The momentum and the energy equations are the same as eqs. (11) and (12) except (x, v, T) are replaced with (xr, vr, Tr). After performing the integral solution, the equation of the warm boundary layer thickness becomes:

$\frac{g\beta }{90\nu }\frac{d}{dy}\left[ {{\lambda }^{3}}{{({{T}_{w}}-{{T}_{c}})}^{2}} \right]+\frac{g\beta ({{T}_{w}}-{{T}_{c}}){{\lambda }^{3}}}{36\nu }\frac{d{{T}_{c}}}{dy}=-2\alpha \frac{{{T}_{w}}-{{T}_{c}}}{\lambda }$	(16)

which is subject to the following boundary condition:

$λ(0) = 0$	(17)

The wall temperature, T_w, can be obtained by considering energy balance at the right side wall of the enclosure:

$\frac{2k({{T}_{w}}-{{T}_{c}})}{\lambda }={{h}_{e}}({{T}_{e}}-{{T}_{w}})$	(18)

The temperature of the core region, Tc, can be obtained by considering the fact that the horizontal entrainment velocities of the two boundary layers represent the same core flow (Bejan, 1989; Zhang and Bejan, 1989).

$δ 3 (T c - T m) = λ 3 (T w - T c)$	(19)

The core temperature is the lowest at the bottom and highest at the top of the enclosure:

$T c (0) = T m$	(20)

$T c (H) = T w (H)$	(21)

Defining the following dimensionless variables:

$\text{Ste}=\frac{{{c}_{p}}({{T}_{e}}-{{T}_{m}})}{{{h}_{s\ell }}},\text{ R}{{\text{a}}_{*}}=\frac{g\beta ({{T}_{e}}-{{T}_{m}}){{H}^{3}}}{\nu \alpha },\text{ }Y=\frac{y}{H},\text{ }\Delta =\frac{\delta }{H}\text{R}{{\text{a}}_{*}}^{\frac{1}{4}},\text{ }\Lambda =\frac{\lambda }{H}\text{R}{{\text{a}}_{*}}^{\frac{1}{4}}$

$\theta =\frac{T-{{T}_{m}}}{{{T}_{e}}-{{T}_{m}}},\text{ Ra}=\overline{{{\theta }_{w}}}\text{R}{{\text{a}}_{}},\text{ Bi}=\frac{{{h}_{e}}H}{k},\text{ }{{U}_{0}}=\frac{{{u}_{0}}H\rho {{h}_{s\ell }}}{k({{T}_{e}}-{{T}_{m}})R{{a}_{}}^{1/4}}$	(22)

the following dimensionless equations are obtained:

$2\text{Ste}\frac{{{\theta }_{c}}^{2}}{\Delta }+\frac{{{\theta }_{c}}}{36}\frac{d}{dY}({{\Delta }^{3}}{{\theta }_{c}})-\frac{1}{60}\frac{d}{dY}({{\Delta }^{3}}{{\theta }^{2}}_{c})=-\frac{2{{\theta }_{c}}}{\Delta }$	(23)

$\frac{1}{90}\frac{d}{dY}\left[ {{\Lambda }^{3}}{{({{\theta }_{w}}-{{\theta }_{c}})}^{2}} \right]+\frac{{{\Lambda }^{3}}({{\theta }_{w}}-{{\theta }_{c}})}{36}\frac{d{{\theta }_{c}}}{dY}=-\frac{2({{\theta }_{w}}-{{\theta }_{c}})}{\Lambda }$	(24)

$Δ 3 θ c = Λ 3 (θ w - θ c)$	(25)

$\frac{2({{\theta }_{w}}-{{\theta }_{c}})}{\Lambda }=BiR{{a}_{*}}^{-1/4}(1-{{\theta }_{w}})$	(26)

and their corresponding boundary conditions are:

$Y=0\begin{matrix} {} & \Lambda =0\begin{matrix} {} & {} \\ \end{matrix} \\ \end{matrix}{{\theta }_{c}}=0$	(27)

$Y=1\begin{matrix} {} & \Delta =0\begin{matrix} {} & {} \\ \end{matrix} \\ \end{matrix}{{\theta }_{c}}={{\theta }_{w}}$	(28)

The dimensionless solid-liquid interface velocity is:

$U 0 = 2θ c / Δ$	(29)

The set of ordinary differential equations (23) and (24), with boundary conditions specified by eq. (27) and (28) is a boundary value problem. The boundary value problem was solved by using a shooting method in the interval of Y = 0 and Y = 1. The objective of the shooting method was to satisfy boundary condition Δ(1) = 0 by properly adjusting Δ(0). Equations (23) and (24) were discretized using an implicit scheme. The discretized equations and eqs. (25) and (26) are solved using an iteration method. The converged solution for the entire problem was obtained when shooting error for Δ(1) was less than 10^-9. Zhang (2002) numerically solved the problem and investigated the effect of Biot number on the wall and core temperatures, the velocity of the solid-liquid interface, and heat transfer characteristics.

References

Bejan, A., 1989, “Analysis of Melting by Natural Convection in an Enclosure,” Int. J. Heat and Fluid Flow, Vol. 10, pp. 245-252.

Cao, Y., and Faghri, A., 1990, “A Numerical Analysis of Phase Change Problem Including Natural Convection,” ASME Journal of Heat Transfer, Vol. 112, pp. 812-815.

Faghri, A., and Zhang, Y., 2006, Transport Phenomena in Multiphase Systems, Elsevier, Burlington, MA

Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.

Gadgil, A., and Gobin, D., 1984, “Analysis of Two Dimensional Melting in Rectangular Enclosures in Presence of Convection,” AMSE J. Heat Transfer, Vol. 106, pp. 20-26.

Ho, C.J. and Viskanta, R., 1984, “Heat Transfer During Melting from an Isothermal Vertical Wall,” ASME J. Heat Transfer, Vol. 106, pp. 12-19.

Wang, S., Faghri, A., and Bergman, T.L., 2010, “A Comprehensive Numerical Model for Melting with Natural Convection,” Int. J. Heat Mass Transfer, Vol. 53, pp. 1986-2000.

Zhang, Y., 2002, “Analysis of Melting in an Enclosure under Boundary Condition of the Third Kind,” Proceedings of the 12th International Heat Transfer Conference, Vol. 3, pp. 315-320, Grenoble, France, August 18-23, 2002.

Zhang, Y., Chen, Z.Q., and Faghri, A., 1997, “Heat Transfer During Solidification around a Horizontal Tube with Internal Convective Cooling,” ASME Journal of Solar Energy Engineering, Vol. 119, pp. 44-47.

Zhang, Z., and Bejan, A., 1989, “Melting in an Enclosure Heated at Constant Rate,” Int. J. Heat Mass Transfer, Vol. 32, pp. 1063-1076.

Retrieved from "http://www.thermalfluidscentral.org/encyclopedia/index.php/Natural_convection_in_melting_and_solidification"

Natural convection in melting and solidification

From Thermal-FluidsPedia

Solidification around Horizontal Cylinder

Melting in a Rectangular Enclosure Heated from the Side

References

Views

Personal tools

Navigation

Search