1 Introduction

In the two last decades, urban growth created some over-populated regions where the demanded electricity could not be delivered by the network. For this reason, enhancing the thermal performance of devices [1] as well as using renewable resources such as solar energy [2, 3] became a field of interest in lots of researches. Changes in electrical energy demand during the day and night have forced the government to set different prices for on-peak and off-peak periods [4]. Thus shifting the energy demand from on-peak to off-peak periods can be economical and it can prevent a power outage. In fact, peak demand management can save 10–15 billion dollars annually in the US market [5]. Shifting the energy demand load from on-peak to off-peak periods of time can be achieved by using different thermal energy facilities such as building thermal mass (BTM), phase change material (PCM) and thermal energy storage system (TES) [6, 7]. By implementing these systems the energy demand can be more balanced, they can also minimize the size and price of the required thermal equipment because while using these TES, the energy will be produced and stored in a lower rate over a longer time and these equipment does not have to be so powerful that it can handle a narrow peak demand [8]. One of the most advanced technologies in the energy field is thermal energy storage (TES) which has drawn much attention recently. Three main types of storing energy in TES are sensible, latent, and thermochemical. In selecting or designing a TES, operating condition, price and the required storage period should be considered. [9]. Khan et al. [10] reviewed studies conducted on the performance of TES-integrated solar collectors comprehensively. They revealed paraffin wax is extensively employed in flat plate solar collectors, and in some cases, nanoparticles were used by researches for enhancing the collector performance. In a comprehensive review article conducted by Qureshi et al. [11] different approaches for enhancing the PCM thermal conductivity were reviewed. Their study mainly focused on additives for improving PCM thermal conductivity and performance.

The applications of PCMs are not limited in building TES units since they have been used in different areas to enhance device efficiencies. Ali et al. [12] experimentally analyzed the thermal performance of different kind of PCM embedded heat sinks, including rectangular, triangular, and circular pin–fin. Besides different heat sink geometries, various types of PCM were evaluated in their study. In the experimental research done by Arshad et al. [13] the effects of different parameters on paraffin wax-based aluminum pin–fin heat sinks enhancements were studied. They compared the results of the traditional pin–fin heat sink and PCM-based pin–fin heat sink. The impacts of PCM volume fraction and pin diameters were evaluated. In the experimental research conducted by Rehman et al. [14], PCM-based copper foam heat sinks were studied. They used different PCMs in both charge and discharge processes. It was revealed that with higher PCM volume fraction leads to a more effective temperature reduction in the base. Ali et al. [15] investigated the ability n-eicosane based circular pin heat sink for cooling portable electrical devices. The investigated parameters in their research included PCM volume, fin configuration, and different heat fluxes. It was shown that the existence of PCM reduces the base temperature and it can be beneficial for electronic cooling. Besides these researches, some other studies [16,17,18,19,20] were also aimed the applications of PCM in heat sinks for electronic system cooling.

An ice storage system is a latent thermal storage system which can be very effective in demand-side management [21]. In recent years, researchers have focused on enhancing the operation and efficiency of the ice storage systems. Yang et al. [22] numerically investigated the thermal effects of coolant inlet temperature on the charging process of an ice storage system. Their results showed that a lower coolant temperature enhances the thermal efficiency as well as ice formation rate. In the experimental and numerical research conducted by Erek and Ezan [23] the influence of different coolant flowrate and temperatures were studied during the charging process of an ice-on-coil ice storage system. Their investigations showed that for the ice formation, the change of coolant temperature is much more influential compared to changing the flow rate. Kim et al. [24] numerically studied the melting process in a cold storage system containing a PCM. They examined the effects of some parameters including length and number of fins on the thermal performance of the storage. Their observations indicated that with larger fins the area that was influenced by these fins were enlarged and the discharge process occurred with a higher speed. Jannesari and Abdollahi [25] studied the effects of annular fins and thin rings on improving the ice formation rate in an ice-on-coil ice storage system. The results indicated that using annular fins or thin rings could improve the formation of ice up to 21 and 34% respectively, compared to when a bare tube was used. Also, the rate of ice formation was increased by 15%. Sheikholeslami [26] investigated the process of water solidification in a latent TES in existence of Copper oxide nanoparticles. His study showed that these additives enhance the PCM thermal conductivity and the rate of solidification. In another investigation conducted by Sheikholeslami [27] the effects of metallic fins length were evaluated for further improving the freezing rate of a nano-enhanced PCM (NEPCM). The transient PCM solidification process near annular fins in a cylindrical shell was studied by Mosaffa et al. [28]. They used Green’s function for acquiring the solution. Zheng et al. [29] examined the influence of thermal resistance, tube diameter and material on an ice storage system. Their study showed that larger tube diameter or utilizing a tube material which has a higher thermal conductivity can enhance the heat transfer and thermal efficiency. Bondareva and Sheremet [30], used numerical methods to simulate the natural convection during the process of melting in PCM-filled storage in both two-dimensional (2D) and two-dimensional (3D) forms. The results indicated that the natural convection was more intense during the 3D simulation. Korti [31] studied the effects of some parameters on enhancing the performance of a latent heat storage system. In his numerical study, the amount of PCM, and water flowrate as the HTF was examined and some guidelines were presented for improving the design of the storage. Mousavi et al. [32] numerically studied the melting process in a shell and coil ice storage system. They evaluated the impacts of different operational and geometrical parameters on the discharge process. Their results showed the dominant heat transfer mechanism in the early stages of the melting process is heat conduction, but the natural convection effects appear when the time passes.

There are two types of ice storage discharge (melt) process: internal and external. In the internal melt process, the coolant flows through coils in both charge and discharge processes, while in an external melt process warm water flows over the ice-covered coils in the storage itself. Since the internal melt process has a more complex structure and the flow the cannot be directly used as it could for the external melt process, the external melt process is more preferable [33].

The heat exchanger is a device which can be utilized in order to transfer heat from one fluid to another one. Heat exchangers are used in different areas such as HVAC, cars, petroleum and some other industries [34]. These devices are so important that a tremendous number of studies have been done in order to enhance them [35]. Different forms of heat exchangers can be found in the literature. Shell and tube, plate, double pipe, heat pipes, and helically coiled heat exchanger are some of them [36]. Due to the secondary vortex flow which is generated inside the helical coils because of their curvature, and also thinner thermal boundary layer in the flow inside them, helical coils tend to provide a better transfer coefficient [37]. Because of compact structure and high heat transfer rate in these heat exchangers [38], they have been widely utilized in different industries like power generation, nuclear technology, heat recovery and refrigeration [39].

Kong et al. [36] examined some geometrical parameters in a helical-coil heat exchanger which was designed for heat recovery applications. These geometrical parameters included pitch length and tube and coil diameters. Their results showed that in low Reynolds numbers, smaller tube diameters enhances the heat transfer, on the contrary, for Reynolds numbers more than 3500, larger tube diameters improve it. Mirgolbabaei [40] studied the performance of a vertical helical coil and evaluated the influence of different flowrate and geometrical parameters such as the dimensionless ratio of coil diameter to tube diameter on the heat transfer. In an experimental examination conducted by Izadpanah et al. [41], the effects of natural convection in a cold storage system equipped with a helical-coil heat exchanger were studied. In their study, R-134a was flowing in the helical coil as a coolant and the cold storage was containing water. They discussed the Nusselt number and Coefficient of Performance (COP) variations due to changes in the number of helical coil turns and coolant flow rate. Zhao et al. [42] studied the heat transfer for membrane helical coil and membrane serpentine tube heat exchangers. In their study the flow was turbulent and the heat transfer was examined under different circumstances such as flow inlet velocities, operating pressures, and pitch lengths. It was shown that the overall performance of the heat exchanger with the membrane helical coil was better compared with the heat exchanger with serpentine tubes. In another numerical study conducted by Zhao and Che [43], the effects of different tube arrangements on hydraulic and thermal aspects of the performance of a helical membrane heat exchanger were studied. They examined how Reynolds number, pitch ratio, and tube arrangements can affect the flow and heat transfer. The results indicated that the helical coil heat exchanger has a substantial effect on heat transfer enhancement.

In the present article, an ice storage system with double helical coil heat exchanger is three-dimensionally modeled. To the author’s knowledge and based on the reviewed literature, there is not a study which investigated the performance of a double coil heat exchanger in an ice storage system, thus using this heat exchanger in the ice storage system is the novelty of this study. The main goal here is to investigate the effects of some geometrical parameters including helical coil pitch length (P) and the distance between inner and outer helical coils (H) on the charging process (solidification) in an ice-on-coil ice storage system. For this purpose, eight different geometries of this heat exchanger with the same heat transfer area are simulated in an ice storage system using the finite volume numerical method. The transient heat transfer and phase change process are modeled for four hours of flow time.

2 Physical model and assumptions

2.1 Physical model

In this study, a 15 L ice storage system has been modeled. The ice storage had the shape of a cylinder with diameter and height of 135 and 262 mm respectively, and it contained pure water. It was tried to consider the storage large enough to prevent the ice from reaching the storage wall because the low temperature can cause deformations in the walls. As it can be observed in Fig. 1, a cold double helical coil heat exchanger was present in the ice storage system for cooling down the water and transforming it to ice. As time passes the water contained in the storage becomes colder and ice layers are formed on the coils. It should be noted this research mainly focuses on the freezing process in the shell. Figure 2 lists the geometrical parameters which can be altered in the double helical coil heat exchanger.

Fig. 1
figure 1

General schematic of an ice storage system with double helical coils

Full size image
Fig. 2
figure 2

The considered double helical coil effective parameters

Full size image

Among these parameters, the effects of two of them, which are coil pitch length (P) and inner and outer coils distance (H) on the phase change process of the pure water in the ice storage system has been studied.

To do so, eight geometries for this heat exchanger was modeled in which the P and H parameters were changed, even though the heat transfer area was considered to be fixed. Figure 3 shows these examined geometries in the ice storage system, and the geometrical details can be found in Table 1. The helical coil height (S) was selected as the free parameter to keep the heat transfer area constant. While keeping the pitch length fixed, the helical coil height was changed by altering the number of helical coil turns (N). The helical coil length was calculated by:

$$S = N*P$$
(1)

Also, the heat transfer area for one of the coils was obtained using:

$$A = N\pi D\sqrt {P^{2} + (\pi D_{c} )^{2} }$$
(2)

It should be noted that the pitch length for the inner and outer helical coils was the same. Also, in the geometries for the H parameter, the diameter of the outer coil was kept fixed while the inner one was changed to study the effects of the distance between inner and outer coils.

Fig. 3
figure 3

Schematics of the models

Full size image
Table 1 Detailed geometrical specifications for different cases
Full size table

The ice storage contained water and the material used for the tubes was copper. It was tried to use standard copper tube diameters. The physical properties for the materials used in this study have been presented in Table 2. These data have been extracted from published articles [25, 44].

Table 2 Thermal and physical properties of used materials [25, 44]
Full size table

2.2 Assumptions

The enthalpy porosity method has been used for simulating the phase change process in the storage [45, 46].

The following assumptions were applied to the problem:

  1. 1.

    The flow is considered to be incompressible, laminar, and transient.

  2. 2.

    The ice storage outer wall is considered to be completely insulated.

  3. 3.

    Although the effects of buoyancy on natural convection were included by the Boussinesq approximation, the effects of density variation in the fluid and solid phases were ignored and an average value of 958.4 kg/m3 was applied to both solid and liquid phases of the PCM. Thus, the effects of heat conduction and natural convection were both considered during this simulation.

  4. 4.

    Even though the effects of temperature on the physical properties such as the thermal conductivity and specific heat capacity were neglected for liquid and solid phases, two different values were applied, i.e. the effect of phase change on the physical properties including thermal conductivity and specific heat capacity was considered.

  5. 5.

    In order to reduce the computation cost and the number of grids, the difference between inlet and outlet coolant temperatures was neglected and a constant wall temperature boundary was applied on the tube inner walls. This assumption is proved to be valid for high coolant flowrates and it has been utilized in both phase change process [47, 48] and shell and coil heat exchanger [49] published studies.

3 Mathematical formulation and simulation procedure

In this section, the governing equations and the procedure of the numerical simulation will be discussed.

3.1 Governing equations

In order to include the effects of buoyancy on natural convection in the storage, the Boussinesq approximation was applied:

$$\rho = \rho_{0} \left[ {\beta \left( {T - T_{liquid} } \right) + 1} \right]^{ - 1}$$
(3)

Continuity:

$$\nabla \cdot \overrightarrow {V} = 0$$
(4)

Momentum:

$$\frac{\partial V}{\partial t} + \overrightarrow {V} \cdot \nabla \overrightarrow {V} = \frac{1}{\rho }( - \nabla \rho + \mu \nabla^{2} \overrightarrow {V} + \rho \beta (T - T_{ref} )) + \overrightarrow {S}$$
(5)

Energy:

$$\frac{{\partial h_{sens} }}{\partial t} + \frac{{\partial h_{lat} }}{\partial t} + \nabla \cdot (\overrightarrow {V} h_{sens} ) = \nabla \cdot \left( {\frac{k}{{\rho c_{p} }}\nabla h_{sens} } \right)$$
(6)

The material total enthalpy which can be computed by the sum of sensible and latent heat:

$$h_{tot} = h_{sens} + h_{lat}$$
(7)

In which:

$$h_{sens} = h_{ref} + \int\limits_{{T_{ref} }}^{T} {c_{p} dT = h_{ref} + c_{p} } \int\limits_{{T_{ref} }}^{T} {dT}$$
(8)
$$h_{lat} = \sum\limits_{i = 1}^{n} {\lambda_{i} h_{sf} }$$
(9)

In which hlat can vary from zero (for solid-phase) to hsf (for the liquid phase) and λ is the liquid fraction. According to the correlation, the total latent heat can be computed by the summation of latent heat from every cell in the computational domain in each time step. After that, the sensible heat can be obtained from:

$$h_{sens\,} = h_{tot} - h_{lat}$$
(10)

And the liquid fraction, λ, can be defined as [50]:

$$\lambda = \left\{ {\begin{array}{*{20}l} {\frac{{h_{lat} }}{{h_{sf} }} = 0} \hfill & {\quad if\quad T \le T} \hfill \\ {\frac{{h_{lat} }}{{h_{sf} }} = 1} \hfill & {\quad if\quad T \ge T_{liq} } \hfill \\ {\frac{{h_{lat} }}{{h_{sf} }} = \frac{{T - T_{s} }}{{T_{liq} - T_{s} }}} \hfill & {\quad if\quad T_{s} < T < T_{liq} } \hfill \\ \end{array} } \right\}$$
(11)

Darcy’s law damping term is usually added to the momentum equation for considering the influence of the phase change in convective heat transfer. The term \(\overrightarrow {S}\) in Eq. 5 is the source term presenting Darcy’s law damping term, and it can be defined as:

$$\overrightarrow {S} = - \frac{{\left( {1 - \lambda } \right)^{2} }}{{\lambda^{3} }}C_{mush} \overrightarrow {V}$$
(12)

where Cmush is the mushy zone constant. This constant is usually varying in the range of 104–107 kg/m3 s. In this study, Cmush is assumed to be constant and its value was set to 105 kg/m3 s.

3.2 Boundary and initial conditions

The ice storage outer walls are considered to be insulated. Besides, a coupled boundary condition has been applied between the tube outer walls and the ice storage contact region with the tubes. A constant temperature of 263.15 K was applied on the tube inner walls. The initial temperature for both tube and ice storage is 274.15 K, one degree over the solidus temperature of the pure water contained in the ice storage.

3.3 Numerical and solution method

The numerical simulation of the phase change process was done by using the Finite Volume Method (FVM). For doing that ANSYS Fluent 19.2 commercial CFD software was employed.

In order to study the phases, change process in the problem, the enthalpy-porosity strategy was utilized. Also, for keeping the simulation accurate, the modeling was done in double precision mode, and the residuals for the continuity, x, y and z velocities were considered to be 10−3 while for the energy equation it was 10−6. In the solution for pressure–velocity coupling SIMPLE algorithm was used while for gradient spatial discretization and transient formulation, least-square cell-based and first-order implicit were applied to the simulation, respectively. Also, for discretizing both energy and momentum equations QUICK scheme was applied.

It should be noted the values for under-relaxation factors were kept at the default values of this software which are 0.3, 1, 1, 0.7, 0.9 and 1 for pressure, density body forces, momentum, liquid fraction and energy respectively.

3.4 Model validation

Since there are no experimental or numerical studies with similar geometries, in order to check the simulation validity, the geometry and conditions from the experimental work of Sasaguchi et al. [51, 52] were simulated using the present model. For this purpose, the solidification process around two cylinders with constant temperatures was numerically modeled. The comparison has been shown in Fig. 4. It can be observed that the results from their experimental data and our numerical simulation are in good compatibility and the maximum deviation between their experimental data and this simulation is only 8.13%. The agreement between these results shows the model is highly functional, thus it can be applied to other cases.

Fig. 4
figure 4

The model validation with the experimental work of Sasaguchi et al. [51, 52]

Full size image

3.5 Grid and time step size independence tests

Grid independence tests were done for Case 4 model with 1.3, 1.8 and 2.3 million of grid numbers, the results for these tests have been presented in Fig. 5a. As it can be observed since the deviation between the liquid fraction plots of simulations with two last grid sizes was neglectable, this grid size was selected and it was expanded to other cases.

Fig. 5
figure 5

Grid (a) and (b) time-step size independence tests

Full size image

Figure 6 shows the generated grid for Case 4. It can be seen that the grid sizes near the delicate areas with a high gradient (near the contact region between the tubes and ice storage) was selected to be smaller.

Fig. 6
figure 6

The selected grid for case 4

Full size image

Also, for the time-step size independence tests were done for case 4 with time-step sizes of 2, 1 and 0.5 s. The results for the liquid fraction with these different time-step sizes can be found in Fig. 5b. Since the results for simulations with the time step (TS) size of 1 and 0.5 s were almost the same, the time-step size was selected to be 1 s. It should be noted that for obtaining precise results for each time step 100 number of iterations were done.

4 Results and discussions

4.1 The effects of pitch length (P)

In this section, the effects of pitch length (P) on the liquid fraction, ice formation and the temperature distribution in the storage will be discussed. Figure 7 presents the influence of helical coil pitch length on the rate of reduction in the liquid fraction. Obviously, a lower liquid fraction shows a higher rate in ice formation and would be better for the ice storage charging purpose. Comparing the liquid fractions for the case 1–4, with the pitch lengths of 36, 42, 48 and 54, it can be observed that the case with higher pitch length (P = 54) has a better performance in ice formation. The difference for the liquid fraction for the cases with P = 34 and P = 54 at the end of the simulation (t = 240 min) is 10.70%. Since these values cannot show the quality of ice formation and the distribution of ice in the storage, the liquid fraction contours can be used for this purpose.

Fig. 7
figure 7

The effects of helical coil pitch length on the liquid fraction

Full size image

Figure 8 shows liquid fraction contours for cases 1–4, with helical coil pitch length varying from 36 to 54 mm. Based on these liquid fraction contours it can be observed that with higher pitch length, a better distribution of ice is achieved in the storage. Besides this, the formation of ice block which is an undesirable phenomenon in the melting process can be retarded with using a higher pitch length. Without the formation of the ice block, in the process of external melt, when a flow of water enters the storage to get cooled, the flow can pass through the coils and this way the heat exchange can be done with a higher rate and efficiency. It was also observed that there was no significant difference between ice formation in lower and upper regions of the storage which indicates that the natural convection did not have a considerable influence on the solidification process. The reason can be traced back to the dominance of conduction heat transfer mechanism in solidification processes which was stated in previous researches [53, 54].

Fig. 8
figure 8

The effects of helical coil pitch length on contours of the liquid fraction (λ) at X = 0

Full size image

The contours of temperature for different times and different coil pitch lengths have been shown in Fig. 9. It can be observed that for all of the models, with time passing, more areas have reached to solidification temperature and the phase change has occurred in a wider region. It is obvious that the temperatures in some areas got even lower than the freezing temperatures. It can be demonstrated that with increasing the pitch length of inner and outer coils, a larger region around coils was exposed to the coolant temperature, thus the low temperature was observed in a wider area. In other words, increasing the inner and outer coil pitch lengths, created the freezing temperature in a wider region, so the ice formation and the stored energy could be increased. This result confirms previous results for liquid fraction plot and contours from Figs. 7 and 8.

Fig. 9
figure 9

The effects of coil pitch length on contours of temperature at X = 0

Full size image

It worth noting that in the initial stages of the solidification process the natural convection effect forces the colder fluid to go upward in the storage. The direction of this movement is affected by the density inversion of water in the temperature range of 273–277 K which was applied to the problem by the negative thermal expansion coefficient of water for this range as presented in Table 2. These movements caused by buoyancy forces are due to local temperature difference in the fluid phase and they tend to fade away as the time passes and the fluid phase reaches to a uniform temperature of 273.15 K (solidus temperature). After reaching this temperature natural convection does not affect the phase change process since with a uniform temperature in fluid phase buoyancy forces do not exist. So, the freezing process continues relying on the conduction heat transfer mechanism.

4.2 The effects of inner and outer coil distance (H)

In this section, the effects of inner and outer coil distance (H) on the liquid fraction, ice formation and the temperature distribution in the storage unit will be discussed. Figure 10 shows the influence of the distance between the inner and outer coils on the liquid fraction reduction rate. Comparing the liquid fractions for the case 5–8, with the coil distances of 45, 50, 55 and 60 mm, it can be observed that the case with the higher coil distance (H = 60) has a better performance in ice formation. The fact that the difference for the liquid fraction for the cases with H = 45 and H = 60 at the end of the simulation (t = 240 min) is only 5.11% shows that although the coil distance is effective on the liquid fraction, but its effects on the liquid fraction is not as high as it was for the pitch length.

Fig. 10
figure 10

The effects of inner and outer coil distance on the liquid fraction

Full size image

Figure 11 shows liquid fraction contours for cases 5–8, with inner and coil distances varying from 45 to 60 mm. Based on these liquid fraction contours it can be observed that with an increase in the coil distance, a better distribution of ice is obtained in the storage. Also, as discussed before, the formation of an ice block can be retarded by using a larger inner and outer coil distance.

Fig. 11
figure 11

The effects of inner and outer coil distance on the contours of the liquid fraction (λ) at X = 0

Full size image

The temperature contours for the areas around the coil for different inner and outer coil distances and times have been presented in Fig. 12. According to these contours, for all of the cases, with time passing, a wider region around the coils reaches to the solidification temperature. The temperature reduction even goes beyond freezing temperature. It can be demonstrated that with increasing the distance between inner and outer coils, a larger region gets exposed to the cold temperature from the coil, and the low temperature can be observed in a larger area. In fact, increasing the coil distance affects a wider region with the freezing temperature, thus it increases the stored energy in the form of ice.

Fig. 12
figure 12

The effects of inner and outer coil distance on contours of temperature at X = 0

Full size image

5 Conclusion

In this study, a 3D numerical simulation was done in order to investigate the effects of two geometrical parameters of the double-helical coil heat exchanger, including helical coil pitch length, and inner and outer coils distance, on the phase change process in an ice storage system. The simulation was done with different thermal and physical properties for water and ice, and the simulation time was 4 h. The results indicate that with higher pitch length and inner and outer coil distances, a better distribution of ice can be created in the storage. Furthermore, the rate of ice formation with highest values of pitch length and inner and outer coil distances compared to cases with smallest values for these parameters were increased by 22.81% and 13.99% respectively, and the possibility of ice block formation which is an undesirable phenomenon in the external discharge process was reduced.