CFD Model of Shell-and-Tube Latent Heat Thermal Storage Unit Using Paraffin as a PCM

Abstract

This chapter validates the capability of CFD modelling technique to accurately describe processes in the thermal storage system with the PCM. For validation purposes, CFD modelling using FLUENT ANSYS was conducted and the predicted results were compared with the experimental and numerical data from the literature. The comparison between experimental and numerical results was carried out in terms of the temperature distributions and average volume of the PCM liquid fraction. Additionally, the detailed parametric study of the storage system with the PCM was performed and results obtained were discussed with dimensional correlations for the Nusselt number being proposed to be used in the designing process. Finally, a correlation was developed to estimate the total melting time at the thermal storage system.

1. Introduction

The thermal energy storage systems can be classified into several main groups, namely thermochemical storage, sensible heat storage and latent heat storage, or combination of these both [1]. The energy is stored in the latent heat storage systems regarding the phase change of materials as a constant or nearly constant temperature. It is critically reviewed in several recent publications that the PCM thermal storage system is one of the most efficient heat storage methods because it provides a considerable amount of energy during the charging and discharging process compared to that of the conventional sensible heat energy storage.

This means that the latent heat storage system using a PCM requires a much smaller volume of materials to store a certain amount of energy. Recently, several studies have been carried out to study the thermal behaviour of PCM thermal storage systems using FLUENT. Al-Maghalseh [2] conducted a comprehensive review of the heat transfer enhancement methods in the thermal storage systems using PCM. Several techniques of heat transfer intensification methods were discussed in details, including both experimental and numerical studies for fins, fins materials and geometry, filling materials, nano-fluid, nano-particles, microencapsulation and thermal conductivity enhancement. Another study by the same author [3] numerically studied the effects of natural convection on the heat transfer of the PCM thermal storage system. He found that the natural convection has a considerable effect on the heat transfer inside the storage system, and therefore reducing the total melting time of the PCM. A detailed review of energy storage using PCMs has been performed in [4, 5, 6, 7, 8] as well as Perez-Raya [9], Joybari [10], Riahi [11], Kozak [12], Liu [13], and Mario [14].

Jian et al. [15] developed a numerical model to predict the transient thermal behaviour during charging and discharging processes of a latent thermal storage unit involving a triplex concentric tube with the PCM filling the middle channel Freeman et al. [16] carried out some investigations into a small-scale solar organic ranking cycle (ORC) with integrated PCM thermal energy storage(TES) unit. The system was examined for selected months in the contrasting climates of Cyprus and the UK. The performance indicator of the ORC engine and the required TES volume with and without the PCM are compared and discussed. It was found that the system with evacuated flat-plate collectors has a better performance compared with using low-cost evacuated-tube heat-pipe collectors. Furthermore, using PCMs for the TES shown better performance and a smaller equivalent storage volume than water.

Dal Magro et al. [17] used the PCM to improve the efficiency of the ORC system operating under thermal power fluctuations. He found that using the PCM allows the capacity factor to increase from 38–52% and the average thermal efficiency to increase from 15.5% to 16.4%. Sagar et al. [18] developed a numerical model for ORC based solar thermal power plant integrated with latent heat thermal energy storage system. Shell and tube latent heat thermal storage system was designed to generate 200 kW during the discharging process. However, the overall performance of the solar thermal power plant was evaluated for ten days of operation. Another study by Manfrida et al. [19] developed a simulation model for PCM thermal storage system coupled with solar-powered ORC. The study examined the thermal performance of the system over several conditions. Further, a case study for the operating of the system during one week was numerically examined. The results clearly showed that the system we able to provide power in 78.5% of the time, with weekly averaged efficiencies of 13.4% for the ORC unit, and 3.9% for the whole plant. A simple numerical method, called, the temperature and thermal resistance iterations, was used in the numerical calculation. The data from the numerical model was then compared with experimental results, and a good agreement was observed. Ho and Chen [20] also developed a numerical model for the melting of ice around a horizontal isothermal cylinder. The model’s results were compared with experimental data published by White in [21], and a good agreement was found. It was concluded that the melting process of ice is strongly affected by the changing recirculation occurring in the molten water. Another numerical model of melting around a horizontal pipe was developed by Rieger et al. [22]. The numerical solution was obtained for Rayleigh numbers (Ra) up to 1.5 × 10⁵, Stefan numbers in the range of 0.005 ≤ Ste ≤0.08, and for Pr=50. It was found that the natural convection is the dominant process in the heat transfer mechanism throughout almost the entire melting process.

Trp studied the transient heat transfer in the shell-and-tube thermal storage system in an experimental and numerical study [23, 24]. He developed a mathematical model based on the non-isothermal phase transition, and it was implemented as a FORTRAN computer code. The numerical results were validated with experimented data, and it was concluded that heat transfer from the HTF to the PCM was low due to the large Prandtl numbers of the HTF. Therefore, a large amount of heat was carried downstream with the HTF, whilst a small amount of heat was transferred to the PCM upstream. The same author [25] numerically investigated the effects of several geometrical parameters and different HTF operational conditions on heat transfer during both melting and solidification processes by measuring the transient temperature distribution of the HTF, PCM and tube wall.

This chapter presents the results of the 3-D CFD modelling of the PCM in the shell-and-tube thermal storage system. Then, the numerical results obtained by the CFD were compared with the experimental and numerical data from the literature. Finally, a detailed parametric study of heat transfer processes in the melting PCM was carried out and results were discussed.

2. Validation of the FLUENT model with the experimental case study by Lacroix

Lacroix [26], conducted a series of experiments to study the heat transfer performance of the shell-and-tube thermal storage unit using PCM. The PCM installed on the shell side, while the Heat Transfer Fluid (HTF) flowing inside the tube. The effects of several thermal and geometric parameters on the heat process were investigated. The schematic diagram of the model is presented in Figure 1. The PCM fills the shell with the diameter of D_e, whereas the HTF flows through the tube with a diameter of Di. The PCM is commercially available material n-Octadecane. The thermophysical properties of the n-Octadecane are presented in Table 1.

Figure 1.

Schematic representation of the test unit.

Properties	Value
Liquid density, ρl	814 kg/m^3
Solid density, ρs	814 kg/m^3
Liquid Thermal conductivity, kl	0.148 W/(m^oC)
Solid Thermal conductivity, ks	0.358 W/(m^oC)
Liquid specific heat, Cpl	2200 J/Kg ^oC
Solid specific heat, Cps	1900 J/Kg ^oC
Latent heat, L	243.5 KJ/Kg
Viscosity, μ	$3.878\times {10}^{-3}$Pa s
Thermal expansion coefficient, β	0.00091 1/K
Melting Temperature, Tm	300.7 K

Table 1.

Thermo-physical properties of the n-Octadecane [27, 28, 29].

Figure 2 illustrates the test unit scheme. Two concentric tubes were used. The inner tube (D_i = 12.7 mm, D_o = 15.8 mm, and L = 1 m) is made of copper, and outside tube (D_i = 25.8 mm, and L = 1 m) is made of Plexiglas. Thick pipe insulation (Rubates Armstrong Armflex II) was used to isolate the system. The space between two tubes was filled with N-Octadecane as a PCM, while water was used as an HTF. An electrical heater inside a tank was used to maintain the inlet temperature of the HTF. Then, the HTF was circulating inside the copper tube with the mass flow ranging from 0.03 to 0.07 kg/s. Three thermocouples were used to record the temperature inside the PCM at different locations. Further, two thermocouples being used to record the inlet and outlet temperature of the HTF. A data acquisition unit was used to record the thermocouples signals into a PC. Finally, the storage unit was positioned vertically to depress the natural convection effects on the heat transfer inside the system.

Figure 2.

Experimental test unit.

For the validation purpose, Lacroix’s experiments were numerically restudied using the ANSYS FLUENT software. In the preliminary simulations, different grid sizes and time steps were carefully examined to obtain computational grid convergence. The computational grid was constructed using 282504 hexahedral elements and boundary layers were used surrounding the pipe. Transient simulations were run using the k-epsilon turbulence model and the time step used in calculations was set to 0.1 s. To study the phase change phenomena in the PCM, the solidification/melting model was enabled. The first-order upwind spatial discretization and the pressure solver with the PRESTO algorithm for pressure–velocity coupling were selected to obtain a converged solution. Convergence criteria were established by setting the absolute residual values to 10⁻⁶ for energy and 10⁻³ for all other variables. Zero heat flux boundary conditions were set on all sides of the shell. The mass flow rate and temperature of the HTF were specified at the inlet of the copper pipe. The mathematical formulations for solving PCM related problems have been categorized [30] as fixed grid, variable grid, front-fixing, adaptive grid generation, and enthalpy methods. Two methods are used to analyse the heat transfer in solid-liquid PCMs. These are the temperature-based and enthalpy-based methods. In the former, temperature is considered to be a single dependent variable. The energy equations for both solid and liquid are formulated separately; and thus the solid-liquid interface positions can be tracked easily to achieve an accurate solution for the problem [31, 32].

An enthalpy-porosity method is used for modelling the solidification/melting process [33]. This technique is described in detail by Voller and Prakash [32, 34].

The energy conservation equation for this case is written as:

$$\frac{\partial }{\partial t}\left(\mathit{\rho H}\right)+\nabla .\left(\rho \overrightarrow{v}H\right)=\nabla .\left(k\nabla T\right)+S$$E1

The enthalpy of the material is calculated as the sum of the sensible heat, h, and latent heat, ΔH:

$$H=h+\mathrm{\Delta }H$$E2

The sensible heat is calculated as:

$$h=h_{\mathit{ref}}+\underset{T_{\mathit{ref}}}{\overset{T}{\int }}{c}_p\mathit{dT}$$E3

The latent heat is also calculated as:

$$\mathrm{\Delta }H={\beta }_{l}L$$E4

The liquid fraction, ${\beta }_l$, can be calculated as:

$$\begin{array}{l}{\beta }_l=0,\mathit{\text{when}}T<T_{\mathit{\text{solid}}}\\ {\beta }_l=1,\mathit{\text{when}}T>T_{\mathit{\text{solid}}}\\ {\beta }_l=\frac{T-T_{\mathit{\text{solid}}}}{T_{\mathit{\text{liquid}}}-T_{\mathit{\text{solid}}}}\mathit{\text{if}}{T}_{\mathit{\text{solid}}}<T<T_{\mathit{\text{liquid}}}\end{array}$$E5

The solid and liquid temperatures are also calculated as

$$T_{\mathit{\text{solid}}}=T_{\mathit{\text{melt}}}+\sum _{\mathit{\text{solutes}}}{K}_i{m}_i{Y}_i$$E6

$$T_{\mathit{\text{liquid}}}=T_{\mathit{\text{melt}}}+\sum _{\mathit{\text{solutes}}}{m}_i{Y}_i$$E7

The source term in the momentum equation can be written as [33]:

$$S=\frac{\left(1-\beta \right)}{\left({{\beta }_l}^3+ϵ\right)}{A}_{\mathit{\text{mush}}}\left(\overrightarrow{v}-{\overrightarrow{v}}_p\right)$$E8

Due to Darcy’s law damping terms as a source term are added to the momentum equation because of the effect of phase change on convection, whereas ε is a small constant number (0.001) used to prevent division by zero and $A_{\mathit{\text{mush}}}$ is the mushy zone constant. Values between 10⁴ and 10⁷ are recommended for most computations [33]. In the present study, the mushy zone was set to 10⁵. ${\overrightarrow{v}}_p$ is the solid velocity due to pulling solidification materials out of the domain; and in the present study, pull velocities are not included in the solution and so ${\overrightarrow{v}}_p$ is set to zero. More details about the numerical model can be found at [32].

The liquid velocity can be calculated by the following Equation [33]:

$${\overrightarrow{v}}_{\mathit{liq}}=\frac{\left(\overrightarrow{v}-{\overrightarrow{v}}_p\left(1-{\beta }_l\right)\right)}{{\beta }_l}$$E9

The validation of the CFD model was carried out by comparing numerical results from ANSYS FLUENT to experimental data obtained by Lacroix [26]. The comparison was carried out for three different cases during the melting process. These are for three different HTF inlet temperatures above the melting temperature of n-Octadecane by 5, 10 and 20 K. The HTF mass flow rate was maintained at a constant value of 0.0315 Kg/s. Figures 3-5 show the temporal temperature variations in the experiment and CFD data at locations T1 (h = 0.51 m, r = 0.002 m) and T2 (h = 0.95 m, r = 0.001 m) inside the PCM. It is clearly shown that the predicted numerical results of the temperature follow the experimental trend in Lacroix [26]. The main discrepancies between the numerical and experimental results can be attributed to the measurements uncertainties and the difference in the PCM physical properties in the solid and liquid phase. Furthermore, CFD numerical results on the evolution of the liquid fraction in the PCM were compared with calculations of Lacroix [26]. Figure 6 shows the variation of molten volume fraction of the PCM in the test unit as a function of time. The outer diameter (D_e) of the storage unit is 22 mm, the inside tube diameter (Di) is 12.7 mm, and the storage length (L) is 1 m. The HTF mass flow rate in simulations ranges from 1.5 × 10⁻⁴ to 1.5 × 10⁻² and the HTF inlet temperature was 20 K above the PCM melting temperature. It can be seen in Figures 7-9, the current CFD results are in a very good agreement with calculations of Lacroix [26]. In general, comparison of CFD results with results presented in Lacroix [26] demonstrated that the developed CFD model accurately describes processes taking place in the experimental test rig and therefore can be used with confidence for further transient heat transfer simulations in the shell-and-tube latent thermal storage unit.

Figure 3.

The variation of the predicted and experimental temperature at locations T1 and T2 (T_in = T_m + 20 K), the mass flow rate = 0.0315 kg/s.

Figure 4.

The variation of the predicted and experimental temperature at locations T1 and T2 (T_in = T_m + 10 K), the mass flow rate = 0.0315 kg/s.

Figure 5.

The variation of the predicted and experimental temperature at locations T1 and T2 (T_in = T_m + 5 K), the mass flow rate = 0.0315 kg/s.

Figure 6.

The CFD liquid fraction variation in time against the numerical liquid fraction from [26] at several mass flow rate.

Figure 7.

Temperature distribution in the PCM for (A) bottom section plane, (B) middle section plane, and (C) top section plane (T_in = 320 K, HTF mass flow rate = 0.0315Kg/s., the elapsed time is 350 sec).

Figure 8.

PCM temperature versus time at the different radial positions: R = 0.001, 0.002, and 0.004 m from the axis of the computational domain at y = 0.5 m, T_in = 320 K, HTF mass flow rate = 0.0315Kg/s.

Figure 9.

PCM temperature versus time along the axis (y = 0.05, 0.49, and 0.95 m) for the radial distances z = 0.001 and 0.004 m, T_in = 320 K, HTF mass flow rate = 0.0315Kg/s.

3. Heat transfer performance

The heat transfer performance of the storage unit was numerically examined during the charging process. Figure 7 shows the numerical results of the temperature and liquid fractions along the storage unit axis during the charging process when the elapsed time is 350 seconds (the inlet temperature of the HTF is 320 K, the mass flow rate is 0.0315 Kg/s).

It can be seen that the highest temperature of the PCM can be observed at the domain’s top region close to the inlet of the HTF and can rise gradually in the regions closed to at the vicinity of the tube walls. Therefore, the top part of the domain converted into a liquid first and later on, melting expands to lower regions on the domain. Figure 8, illustrates the temperature variation in the PCM for three different radial locations. As expected, the higher temperatures are noticed in the regions close to the surface of the wall with the HTF, where the melting process takes place first. The temperature variations along the axis at the outer surface of the HTF tube are shown in Figure 9. It can be seen that higher temperatures exist at the domain’s top. This is mainly because the HTF flows from the top to the bottom and thus the temperature rises faster at the vicinity of the tube walls close to the inlet. At this period, the effect of natural convection is not profound yet.

The effect of HTF inlet temperature on the melting process is demonstrated in Figure 10. It can be seen that the inlet temperature of the HTF considerably affects the rate of the melting and the PCM temperature distribution. The increase in the inlet temperature of the HTF leads to a rise in the temperature difference between the tube walls and the bulk of the PCM and thus enhance the heat transfer rate. It results in the faster rise of the liquid fraction and decreases the total melting time. Figure 11 shows that the time for completion of melting for the HTF inlet temperature of 305 K is 2974 s, for 310 K this time is reduced to 1617s and finally, for 320 K, the time of melting is only 963 s. Therefore, the total melting time is reduced approximately by 68% when the inlet temperature is increased from 305 to 320 K and by 45.6% when the inlet temperature is increased from 305 to 310 K.

Figure 10.

Effect of the HTF inlet temperature on the melting process: (A) location T1 (h = 0.51 m, r = 0.002 m) and (B) location T2 (h = 0.95 m, r = 0.001 m), T_intial = 282.5 K, HTF mass flow rate = 0.0315Kg/s.

Figure 11.

The effect of the inlet temperature of the HTF on the formation of the liquid fraction of the PCM and on the melting time.

Figure 12 demonstrates the effect of the HTF flow rate on the meting process. It can be seen that the flow rate accelerates the melting process due to the increased heat transfer rate. It can be seen from Figure 13 that when the flow rate increases from 0.000315 to 0.00315 kg/s, the PCM melting time is reduced from 2781 to 1173 s (reduction by 57%). Also, the melting time is reduced by 17.8% when the mass flow rate increases from 0.00315 to 0.0315 kg/s. The effect of the HTF mass flow rate rise is less profound when compared with the effect of the rise in the HTF inlet temperature. This is demonstrated in Figure 14 using temperature distribution counters and the PCM fluid fraction evolution diagrams.

Figure 12.

Effect of the HTF flow rate on the melting process. The PCM temperature at locations (A) T1 (h = 0.51 m, r = 0.002 m) and (B) T2 (h = 0.95 m, r = 0.001 m), T_intial = 282.5 K, HTF mass flow rate = 0.0315Kg/s.

Figure 13.

The effect of the HTF flow rate on the liquid fraction formation and the melting time, T_intial = 282.5 K, T_inlet = 320 K.

Figure 14.

Temperature distribution and melting process at the bottom of the computational domain for the different HTF mass flow rates and inlet temperatures (elapsed time is 850 s).

4. Effect of natural convection

To study the effect of the natural convection, the system needs to be installed in the horizontal position. Several points were created inside the computational domain to monitor the variation of the temperature inside the PCM during numerical CFD modelling. These monitoring points are placed in three planes, which are perpendicular to the axis of the system and located at distances of 0.07, 0.51 and 0.95 m from the front of the system (see Figure 1). Figure 15 indicates the locations of monitoring points in the plane at a distance of 0.07 m from the front of the system.

Figure 15.

Locations of the monitoring points around the pipe.

Figure 16 shows the temperature variation at some of the monitoring points at the front plane of the system with the PCM for the case when the inlet temperature of the HTF is 320.7 K and the mass flow rate is 0.0315 kg/s. It can be seen that initially, the temperature increases rapidly from 280 to 299 K due to the heat transfer from the pipe walls to the solid PCM by conduction. The temperatures at monitoring points u1, r1, and b1 rise more rapidly due to their proximity to the pipe. Temperatures in monitoring points u3, r3, and b3 rise considerable slower since these points are located on the edge of the storage unit. During the melting process the temperature increases from 299 to 300.7 K and equalise at all monitoring points. Initially, during the melting process, a thin liquid layer is formed between the pipe and the solid PCM. Gradually, the solid–liquid interface expands in the axial and radial directions and the melting then is dominated by natural convection in the PCM’s liquid regions. The melting process is intensified in the upper regions of the container, resulting in higher temperature recordings at the top of the computational domain (point u1).

Figure 16.

Variation of the base PCM temperature with time.

The velocity vectors in the liquid pure PCM are shown in Figure 17 for the elapsed time of 300, 400 and 550 s. It can be seen that the molten PCM ascends upwards from the top regions at the centre of the unit and then after cooling flows downwards to complete the natural convection circle. The convection is intensified as the liquid fraction volume increases. The velocity magnitude gradually decreases in time due to a reduction in the temperature difference in the molten PCM. These results are in good agreement with the results of several experimental and numerical investigations [35, 36, 37].

Figure 17.

Velocity vectors in pure PCM.

5. Average heat transfer coefficient

The local heat transfer coefficient could not be estimated accurately for the present thermal storage system as there is a temperature difference between the outer surface of the HTF pipe and the PCM along with both axial and radial directions. Consequently, the average heat transfer coefficient for the melting process is calculated instead using the following Equations [38] for the temporal heat transfer coefficient

$$h_p=\frac{q}{A\times \mathrm{\Delta }T\left(\mathit{\text{LMTD}}\right)}$$E10

where the surface area of HTFP is calculated using the following equation:

$$A=\pi D_{o}l$$E11

The heat transfer rate (q) in the thermal storage can be calculated through the HTF’s enthalpy reduction rate in the HTF [39]. This enthalpy reduction can be calculated through the following equation:

$$q=\dot{m}{c}_p\left(T_i-T_o\right)$$E12

where, $\dot{m}$ is the HTF mass flow, T_i and T_o are the HTF’s inlet and outlet temperatures respectively.

The average heat transfer coefficient is:

$$\overline{h_p}=\frac{Q_{\mathit{\text{total}}}}{A\times \mathrm{\Delta }T\left(\mathit{\text{LMTD}}\right)t_n}$$E13

$$\mathrm{\Delta }T\left(\mathit{\text{LMTD}}\right)=\frac{\left(T_{w1}-T_{u1}\right)-\left(T_{w2}-T_{u2}\right)}{\left[ln\left(T_{w1}-T_{u1}\right)/\left(T_{w2}-T_{u2}\right)\right]}$$E14

where T_u1 and T_u2 are the PCM temperature at points u1 (first measurement plane) and u2 (last measurement plane); T_w1 and T_w2 are the pipe wall temperatures at the same corresponding measurement planes.

Finally, the temporal Nusselt number (Nu) is used to quantify heat transfer and this can be calculated using the following equation:

$$\mathit{Nu}=\frac{h_p{r}_{\mathit{eq}}}{k_{\mathit{pcm}}}$$E15

The time-averaged Nusselt number $\overline{\mathit{Nu}}$ is defined as

$$\overline{\mathit{Nu}}=\frac{\overline{h_p}{r}_{\mathit{eq}}}{k_{\mathit{pcm}}}$$E16

Different cases were analysed with various system’s geometrical and thermophysical parameters. To generalise results it is vital to characterise them in the dimensionless form. More details about the dimensional parameters calculation and analysis can be found at [32].

The heat transfer coefficient values were calculated for the top, side and bottom regions of the storage unit for several cases during the melting processes.. Figure 18 present results on the heat transfer coefficient changing as a function of time for the PCM during the case when the inlet temperature of the HTF is 320.7 K and the mass flow rate is 0.0315 kg/s. It can be seen that the heat transfer coefficients for both top and side regions increase with time. This agrees well with the initial solid PCM’s temperature rise followed by the melting process. The higher heat transfer coefficient values can be observed in the top regions of the system. This can be attributed to the effect of natural convection. For the longer elapsed times, the heat transfer value stabilises. This is because the temperature distribution becomes more established within the PCM body when the system reaches the steady-state operation. At this stage, the majority of the PCM is melted and a very small part of it, which is close to the bottom regions of the unit, maybe in the solid-state.

Figure 18.

The heat transfer coefficient variation at the top, side, and bottom regions of the system.

The average heat transfer coefficient values for various HTF inlet temperature and flow rates were obtained and used to calculate the Nusselt number for the pure PCM. The flow rates of the HTF considered were 0.000315, 0.00315, 0.005, and 0.0063 Kg/s. The HTF inlet temperatures used in the modelling are 305, 310, and 320 K. To derive generic heat transfer correlations, the Nusselt numbers were calculated at the top, side and bottom regions of the storage unit. Thereafter, the average Nusselt number value for all regions is calculated. Figures 19-21 show the results of these CFD modelling. As it can be seen in Figure 19, the Nusselt number increases with the rise in the Stefan number (see Eq. 18) which is proportional to the difference between the inlet temperature of the HTF and the melting temperature of the PCM (the increase in the inlet temperature of HTF increases the Stefan number). Similar observations can be made concerning the Rayleigh number (see Eq. 19), see Figure 20. The Rayleigh number is also proportional to the difference between the inlet temperature of the HTF and the melting temperature of the PCM.

Figure 19.

Variation in the Nusselt number as a function of Stefan number for (A) upper regions (B) side regions (C) bottom regions (D) the average Nusselt number (in the system).

Figure 20.

Variation in the Nusselt number as a function of Rayleigh number for (A) upper regions (B) side regions (C) bottom regions (D) the average Nusselt number (in the system).

Figure 21.

Variation in the Nusselt number as a function of Fourier number for (A) upper regions (B) side regions (C) bottom regions (D) the average Nusselt number (in the system).

$$\mathit{Ste}=\frac{{\mathit{Cp}}_{\mathit{pcm}}\mathrm{\Delta }{T}_H}{{\lambda }_{\mathit{pcm}}}$$E17

$$\mathit{Ra}=\frac{\mathit{g\beta }{\mathit{Cp}}_{\mathit{pcm}}{\rho }_{\mathit{pcm}}\mathrm{\Delta }{T}_H{r}_{\mathit{eq}}^3}{{\mu }_{\mathit{pcm}}{k}_{\mathit{pcm}}}$$E18

$$\mathit{Fo}=\frac{\left({\alpha }_{\mathit{pcm}}{t}_m\right)}{{\left(r_{\mathit{eq}}\right)}^2}$$E19

Figure 21 show that the Nusselt number decreases with the rise in the Fourier number (see Eq. 20). The Fourier number is proportional to the melting time. However, an increase in the inlet temperature and flow rate of the HTF will lead to a decrease in the total melting time, and thus, reduce the Fourier number.

The average Nusselt numbers in the system are presented in Figure 19D,20D, and 21D. This data was used to derive the Nusselt number correlations for the pure PCM. The correlations are derived as a function of the Stefan, Fourier, and Rayleigh numbers.

The Nusselt number for the pure PCM is ($R^2=0.8524$)

$$\mathit{Nu}=2.9883{\left(\mathit{Ste}\right)}^{0.0758}{\left(\mathit{Fo}\right)}^{-0.095}{\left(\mathit{Ra}\right)}^{0.0759}$$E20

The correlation between the numerically obtained and calculated using Eqs. (20) is shown in Figure 22. It can be seen that the Nusselt number varies between 2 and 4.3 for the system under investigations. The lower Nusselt numbers are for the low HTF inlet temperature, and the high Nusselt numbers values are obtained at the high HTF inlet temperatures. The total melting time of the PCM can then be calculated using the following formula (R² = 0.966)

Figure 22.

The correlation between numerically obtained and calculated using Eq. (20) Nusselt numbers - the pure PCM.

$$t_m=\frac{\mathit{\rho D}}{1.056\overline{h}}\left[\frac{c_p\left(T_m-T_{int\mathit{ial}}\right)+{\lambda }_{\mathit{pcm}}}{\left(T_{\mathit{\text{inlet}}}-T_m\right)}\right]$$E21

6. Conclusions

The 3D CFD simulation model was developed for shell-and-tube thermal storage system using FLUENT/ANSYS. For the validation purpose, the results from the numerical model were compared with the experimental results of Lacroix [26]. The results show that a considerable agreement between the numerical and experimental results. Therefore, the results demonstrate that the developed CFD model accurately describes processes taking place in the experimental test rig and therefore can be used with confidence for further transient heat transfer simulations in the shell-and-tube latent thermal storage unit. 3D CFD simulations were performed for a range of the HTF inlet temperature values and its mass flow rate. Thereby, the results were used to derive Nusselt number correlations as a function of Stefan, Rayleigh and Fourier numbers to take into account the effect of all design and operational conditions. The Nusselt number for the system with the pure PCM was found to be ($R^2=0.8524$)

$$\mathit{Nu}=2.9883{\left(\mathit{Ste}\right)}^{0.0758}{\left(\mathit{Fo}\right)}^{-0.095}{\left(\mathit{Ra}\right)}^{0.0759}$$

The total melting time of the PCM can be estimated by (R² = 0.966)

$$t_m=\frac{\mathit{\rho D}}{1.056\overline{h}}\left[\frac{c_p\left(T_m-T_{int\mathit{ial}}\right)+{\lambda }_{\mathit{pcm}}}{\left(T_{\mathit{\text{inlet}}}-T_m\right)}\right]$$

Subscripts

L

Liquid

eff

Effective

f

Fluid

i

Inlet

m

Melting

o

Outlet

s

Solid

T

Total

w

Wall

Greek symbols

ρ

Density ($\mathrm{Kg}/{\mathrm{m}}^3$)

$\beta$

Coefficient of thermal expansion (1/K)

μ

Viscosity ($\mathrm{Kg}/\mathrm{ms}$)

${\beta }_t$

Liquid fraction

α

Permeability

ϕ

Volumetric fraction

Nomenclature

A	Area (m2)
AMush	Mushy zone constant
CP	Specific heat (J/KgoC)
D	Diameter
g	Gravitational acceleration
Gr	Grashof number
h	Sensible enthalpy (J)
H	Enthalpy
H	Height
href	Reference enthalpy
k	Thermal conductivity (W/moC)
keff	Effective thermal conductivity
L	Latent heat (Kj/Kg)
L	Length
LMTD	The logarithmic mean temperature difference
m	Mass in each space (Kg)
ṁ	Mass flow rate
Nu	Nusselt number
Pr	Prandtl number
q	Heat transfer rate
r	Radial coordinate
Ra	Rayleigh number
Re	Reynolds number
S	Source term
Sh	The heat from any volumetric sources
Sm	Mass source
Ste.	Stefan number
t	Time (s)
T	Temperature
Tref	Reference temperature
u	X velocity components (m/s)
v	Y velocity components (m/s)
vn	Normal components of the velocity of the interface
vr	Radial velocity component
vx	Axial velocity component
w	Z velocity components (m/s)
W	Width
x	Axial coordinate
YM	Contribution of the fluctuating dilatation incompressible turbulence to The overall dissipation rate