The rapid development of the industry has led to a remarkable increase in electricity demand, which is especially evident in summer. Thermal energy storage (TES) that can balance energy demand and supply could play an important role in achieving carbon peak and carbon neutrality [1,2]. Among common TES technologies, Latent Thermal Energy Storage (LTES) that uses Phase Change Material (PCM) is gathering momentum because it is regarded as a compromise between a sensible and thermochemical type [3,4]. A relatively high thermal storage density and small volume variation make LTES a good candidate for future application [5,6]. As a common type of LTES, Ice Thermal Storage (ITS) systems could be adopted in a variety of industries, which include food processing, air conditioning, drug delivery, and building energy conservation [7,8,9].
Heat transfer enhancement is also a key research direction for all types of LTES technologies, such as ITS [9]. One method is to add materials with high thermal conductivity to form an advanced nanofluid [10,11]. Yang et al., [12] prepared a nanofluid with graphene oxide to improve the thermal properties of water. Their results indicated that there was a large increase in the thermal conductivity of the nanofluid and the maximum value could reach 48.1%. Xing et al., [13] experimentally investigated thermal conductivities of nanofluids with different Carbon Nanotubes (CNTs). It is noted that CNTs–nanofluids could present the maximum improvements in thermal conductivity, potentially reaching 16.2% for long singlewalled nanotubes. Parker et al., [14] developed a nanofluid by using graphene nanoparticles with high thermal conductivity. Their results demonstrated that nanofluid that incorporated oxidized graphene had outstanding properties as heat transfer media. To further improve the performance of nanofluid, Du et al., [15] investigated nanofluid with iron oxide (Fe3O4) plus multiwalled carbon nanotube. Thermal conductivities were improved by 32.76% and 33.23% at 50°C. Moreover, structure variation is another way to improve overall heat transfer performance. Abhishek et al. [16] investigated and analysed the charging and discharging working process of an ITS system. It was indicated that tube diameter, orientation, and outlet location had a great influence on the performance of thermal storage system. Vyshak et al. [17] investigated the discharging characteristic of PCM with three encapsulated configurations. Their results indicated that the heat transfer medium that flows inside the inner tube had the highest energy storage efficiency. Anica et al., [18] conducted a numerical simulation of an LTES device and analysed the heat transfer performance during the melting processes. The enthalpy method was used to couple transient convection and solid–liquid phase change heat transfer. Soltan et al. [19] analysed the water freezing time around a circular pipe. Different technologies were used to investigate transient mass and heat transfer around the pipe. Kousha et al., [20] investigated heat storage performance of cylindrical ShellandTube Latent Heat Storage (STLHS) units at different inclination angles. They demonstrated that the liquid fraction of horizontal units during the melting process was improved when compared with that of vertical units. To ensure the influence of extended surface on the heat transfer performance, Yuan et al., [21] analysed the melting process of an annular Latent Heat Storage (LHS) unit with fins. They showed that heat transfer efficiency with less melting time could be obtained by using fins, while natural convection in heat storage units was weakened.
Except the two main methods of heat transfer enhancement mentioned above, the effect of natural convection cannot be ignored [22]. A solidification process without considering natural convection was investigated by Chiu et al., [23]. It showed that their experimental results were in good agreement with simulated results. Darzi et al., [24] numerically investigated different structures of LHS units in terms of the melting and solidification processes. They observed that natural convection could play a leading role in the heat release process with Neicosane for PCM. The bottom of the annular cross section had a better melting rate. Since the melting process was affected by buoyancy, the area dominated by natural convection increased as the tube moved downward. Cao et al., [25] experimentally and numerically investigated the heatreleasing process in the eccentric structured STLHS units that were filled with lauric acid. It indicated that the liquid fraction rate could be greatly increased with the increase in eccentricity and the area of the region determined by natural convection. They also found that the inlet temperature has a major effect on the heatreleasing performance. However, ITS is a unique process, and it is related to phase changing and natural convection between ice and water. The special property of density inversion of water occurs at around 4°C [26,27]. It is acknowledged that the warm temperature of liquid, e.g., PCM or water with low salinity floats owing to its lighter density, while the denser solid or colder PCM would sink. Density reversal of the water occurs when the peak density is near 4°C, i.e., density is positively correlated with a temperature below 4°C, while it shows a negative trend at a temperature above 4°C [28]. It results in a significant change in the flow direction of the melting process, which is determined by natural convection. It is key to studying the melting process in the water LHS unit.
In our previous research, a passive heat transfer improvement of PCM in an STLHS unit was investigated to make use of natural convection for ITS [29]. However, the performance, by considering different positions and configurations of heat sources, is not covered and also has not been analysed in the literature. Table 1 shows a comparison between this work and some representative literature. To comprehensively investigate cold storage processes, the melting process of ice is numerically investigated by the constant heat source input in this work. The effects of natural convection on PCM in an STLHS unit, solid–liquid interface, and temperature distribution are analysed and compared in terms of different heat source positions and configurations, i.e., central position, eccentric position, and flattube configuration. The framework of this paper is illustrated as follows: physical and numerical models are indicated in Section 2 that are also validated by the experimental data. Results and discussions are shown in Section 3, followed by the conclusions in Section 4.
Physical and Numerical Models of ITS
Physical Model
The experimental rig (Figure 1a) in this study consisted of an STLHS unit, instrumentation devices, thermocouples, and a control loop. The STLHS unit, which included two concentric cylinders, is shown in Figure 1b. The geometrical parameters are as follows: one, with a diameter (Ri) of the inner tube, was 10 mm and was made of copper with a thickness of 2 mm. The other, with a diameter (Ro) shell of 40 mm, was made of Plexiglas. The length of the STLHS unit (L) was 200 mm, which was wrapped with insulation cotton of a thickness within 15 mm. Wires were used as the endothermic source and were arranged inside the copper tube and regulated by a Proportional–Integral–Derivative (PID) controller as well as a DC power supply. The temperature of the heat storage medium was tested by thermocouples, which were placed at four different locations inside the STLHS unit, with an accuracy of ±0.15 K. The thermocouples were arranged as shown in Figure 1b. Water was used as a lowtemperature PCM medium, which was used to fill the circular space between the two cylinders. The initial temperature of the ice was 267.15 K. The testing result was just used for model validation. Figure 2 shows the different heat source positions and configurations, i.e., central position, eccentric position, and flattube configuration, which were used for the simulation in the rest of the paper.
Figure 1. Schematic and photograph of the experimental rig: (a) schematic diagram of the experimental rig; (b) photo and structure of an STLHS unit [29].
Figure 2. Positions and configurations of heat source for PCM in an STLHS unit.
Governing Equations
Pertinent assumptions are defined as follows:
 The endothermic source is considered as a constant temperature.
 The effect of the thickness of tube and shell on the numerical model was ignored, and the shell was regarded as an adiabatic wall.

Variation of PCM properties of ice and water was ignored except for the density. It implied that only the effect of density variation of PCM with the temperature on buoyancy was considered.

Melting of the PCM was isotropic and homogeneous.

The experimental system had four sets of thermocouples distributed uniformly. The geometry and place occupied by these probes were not defined.

The enthalpy–porosity method, as a common method for PCM simulation [31–34], was applied to simulate the icemelting process.
For PCM, based on the above assumptions, the governing equations of the model were continuity, momentum, and energyconservation equation, which are presented below. The continuity equation was expressed as Equation (1) [29].
The momentum equation was expressed as Equations (2) and (3) [29].
where Vr and Vq denote the flow velocity of the melted water along r and q directions (ms^{1}). The first term on the right side of the equation represents natural convection term according to the Boussinesq approximation, i.e., the density variation is only considered in the gravity term. b is the expansion coefficient (K^{1}),denotes dynamic viscosity, S_{r} and Sq are the damping source terms which are used to vary fluid velocity during the melting process (ms^{1}), and r and q are introduced into the momentum equation to describe the effect of the phase transition on convective region of the PCM, and are calculated by Equation (4) [29].
where A_{mush} represents the mushy region constant, which describes the steepness of the velocity gradient in the mixing region (kg.m^{1}s^{1}), ranging between 10^{4} and 10^{10}, f denotes the liquefaction rate in the range from 0 (solid) to 1 (liquid), which represents the proportion of liquid PCM in total PCM. A_{mush} is assumed to be constant and set to 10^{6}. e is a constant that prevents the denominator from being zero, which is normally taken as 0.001.
Energy conservation equation could be expressed as Equation (5) [29].
where l is thermal conductivity (W.m^{1}1.K^{1}), and H represents the sum of sensible and latent heats formulated by the enthalpy method (kJ.kg ^{1}). Instead of identifying a precise solid–solution interface, this method distributes a given liquefaction rate to each calculation cell according to enthalpy balance [31,35], where H‘ (kJ.kg^{1}) is the sum of sensible enthalpy, and DH (kJ.kg^{1}) is latent heat, as in Equation (6).
Water density varied nonlinearly with temperature and reverses at 277.15 K, which can be considered by the following equation [36]. It can be derived for water densities at different temperatures, as shown in Figure 3.
where r is the transient density (kg.m^{3}), rl,max is the maximum water density with temperature (kg.m^{3}), and T_{max} is the temperature of the maximum density of water. In this study, r1,_{max} is 999.97 kg.m^{3}, T_{max} is 277.15 K, and g is 9.3 10^{6}.
As shown in Equations (2) and (3), the effect of density variation on natural convection was transferred to that of the expansion coefficient. Thus, b in the model is described as Equation (11).
Heat flux on the tube surface was set as constant and the shell was adiabatic from the surrounding environment. Therefore, boundary conditions of the PCM in the STLHS unit are expressed as Equations (12) and (13).
To evaluate the general heat transfer performance of the STLHS unit, the overall heat transfer coefficient, which simultaneously considered heat conduction and convection, was calculated for discussion as in Equation (14).
Figure 3. Variation of water density with the increase in temperature.
Physical properties of water in the simulation are listed in Table 2 [29].
Numerical Model and Validation
The coupled equations above were solved by ANSYS FLUENT software, which developed the userdefined functions to match properties of density dynamics with temperature. Boundary conditions of the external wall of the STLHS unit were assumed to be adiabatic for insulation simulation. A PISO (Pressure Implicit with Split Operator) algorithm was adopted to solve the pressure–velocity coupling, the skewness correction and the neighbourhood correction were set to 1, the second order scheme could be adopted to discretize the momentum and energy option, and the convergence criteria used for continuous equation, momentum equation, and energy equation were all less than 1x 0^{6}.
Figure 4a shows the temperature at various measuring positions of the STLHS unit during the melting process. The dashed line represents the numerical simulation results, whereas the continuous line represents experimental measurements. It was indicated that the model of PCM was generally reliable and well indicated that the experimental situation was reliable. Furthermore, Figure 4b,c show the total melting time of ice when the input heat flux was 3000 W.m^{2} with different simulation time steps and cell numbers. It indicates that there was no significant change when reducing the time steps and cell numbers. Considering computational cost and accuracy, 0.2 s and 12,500 cells were used for the numerical study. For more details of validation, please refer to our previous research [29].
Figure 4. Validation results of melting model: (a) comparison between simulation and experiments [29], (b) independence test of time step, (c) independence test of cell number.
Results and Discussion
Heat Source in Central Position
Figure 5 shows thermal characteristics of PCM in an STLHS unit based on a heat source in the central position. In Figure 5a, the melting evolutions of PCM are presented in terms of different heat flux from 500 W.m^{2} to 3000 W.m^{2} and liquefaction rate distribution when the global average of liquefaction rate increases from 0.1 to 0.9. The left part of each circle is the phase and velocity vector distributions, while the right part is the temperature distribution. Initially, heat conduction plays a leading role between the heat source and the ice. Thus, the melted liquid is symmetrically distributed in the parts that are close to the central position. As the melting process proceeds, the thickness of the melted part increases, and the melting of warm water starts to flow in the STLHS unit. It demonstrates that natural convection is gradually dominant in the melting process with circulation vortices and further promotes the heat exchange process. It is worth noting that the direction of the flow above and below the STLHS unit are determined by different levels of heat flux ranging from 500 W.m^{2} to 3000 W.m^{2}, which is mainly due to the phenomena of water density, as shown in Figure 3. When heat flux is below 1000 W.m^{2}, heat convection in the melting process flows to the bottom of the unit while the flow starts to go to the up part with a heat flux higher than 2000 W.m^{2}. Thus, the critical point for heat convection of water flow ranges between 1000 W.m^{2} and 2000 W.m^{2}. The reason for the influence of heat flux on convection direction is described as follows: due to the unique density feature, the hot water moves downwards, driven by natural convection when the temperature is lower than 277.15 K, and moves upwards when the temperature is higher than 277.15 K. For different levels of heat flux, the temperature of liquid water is different. When heat flux is relatively small, e.g., 500 W.m^{2}, the water temperature mainly ranges from 273.1 K to 277.15 K, while, when heat flux is relatively large, e.g., 3000 W.m^{2}, the water temperature is always higher than 277.15 K. Thus, the natural convection directions of these two conditions are different, which means that the direction of heat convection changes with different heat flux. The overall heat transfer coefficient is presented in terms of different heat flux, as shown in Figure 5b.
Figure 5. Thermal characteristics of PCM in an STLTH unit with a heat source in the central position vs. different heat flux: (a) temperature distribution; (b) the overall heat transfer oefficient; (c) the liquefaction rate.
It is noted that the overall heat transfer coefficient of different heat flux from 500 W.m^{2} to 3000 W.m^{2} decreases sharply with the increase in the melting time at the beginning of the melting process. This is mainly because heat conduction initially plays a leading role in the STLHS unit. Then, the heat transfer coefficient tends to be constant at about 250 W.m^{2}.K^{1} based on heat convection. Figure 5c shows the liquefaction rate of the unit with the increase in the melting time. The results show that the liquefaction rate increases with the increase in the melting time. The higher heat flux is, the larger the liquefaction rate becomes, which leads to a reduced melting time of PCM in the STLHS unit. The melting time of PCM with a heat flux of 500 W.m^{2} is more than five times longer than that of 3000W.m^{2}.
…To be continued
This article was first published in the MDPI journal, MDPI, Basel, Switzerland. Authors retain the copyright. © 2022 by the authors.
Chunwei Zhang and Dongdong Chai belong to Beijing Institute of Aerospace Testing Technology, Beijing 100074, China. Yubin Fan, Wenyun Zhang and Long Jiang belong to Institute of Refrigeration and Cryogenics, Zhejiang University, Hangzhou 310027, China. Meng Yu belongs to Special Equipment Safety Supervision Inspection Institute of Jiangsu Province, Nanjing 210036, China. Zhenwu Wang belongs to The School of Architecture and Civil Engineering, Jinggangshan University, Jinggangshan 343009, China.
References
 Wang, L.W.; TamainotTelto, Z.; Metcalf, S.J.; Critoph, R.E.; Wang, R.Z. Anisotropic
thermal conductivity and permeability of compacted expanded natural graphite. Appl.
Therm. Eng. 2010, 30, 1805–1811. [CrossRef]  Li, Q.; Bai, F.; Yang, B.; Wang, Z.; El Hefni, B.; Liu, S.; Kubo, S.; Kiriki, H.; Han, M.
Dynamic simulation and experimental validation of an open air receiver and a thermal
energy storage system for solar thermal power plant. Appl. Energy 2016, 178, 281–293.
[CrossRef]  Jiang, L.;Wang, R.Z.;Wang, L.W.; Roskilly, A.P. Investigation on an innovative
resorption system for seasonal thermal energy storage. Energy Convers. Manag. 2017,
149, 129–139. [CrossRef]  Jiang, L.; Li, S.; Wang, R.Q.; Fan, Y.B.; Zhang, X.J.; Roskilly, A.P. Performance analysis
on a hybrid compressionassisted sorption thermal battery for seasonal heat storage in
severe cold region. Renew. Energy 2021, 180, 398–409. [CrossRef]  Jiang, R.; Yu, X.; Chang, J.; Yu, X.; Wang, B.; Huang, R.; Li, Z. Effects evaluation of fin
layouts on charging performance of shellandtube LTES under fluctuating heat sources.
J. Energy Storage 2021, 44, 103428. [CrossRef]  Hu, C.; Li, M.;Wang, Y.; Li, G.; Ma, X.; Du,W.; Zhou, X.; Zhag, Y. Preliminary
investigation on pilotscale photovoltaicdriven cold storage with ice thermal storage
based on vapor compression refrigeration cycle. Sustain. Energy Technol. Assess. 2021,
45, 101187. [CrossRef]  Khadiran, T.; Hussein, M.Z.; Zainal, Z.; Rusli, R. Encapsulation techniques for organic
phase change materials as thermal energy storage medium: A review. Sol. Energy Mater.
Sol. Cells 2015, 143, 78–98. [CrossRef]  Sarkar, S.; Mestry, S.; Mhaske, S.T. Developments in phase change material (PCM)
doped energy efficient polyurethane (PU) foam for perishable food coldstorage
applications: A review. J. Energy Storage 2022, 50, 104620. [CrossRef]  Wang, X.; Li, W.; Luo, Z.; Wang, K.; Shah, S.P. A critical review on phase change
materials (PCM) for sustainable and energy efficient building: Design, characteristic,
performance and application. Energy Build. 2022, 260, 111923. [CrossRef]  Fan, Y.; Zhang, C.; Jiang, L.; Zhang, X.; Qiu, L. Exploration on twostage latent thermal
energy storage for heat recovery in cryogenic air separation purification system. Energy
2022, 239, 122111. [CrossRef]  Fan, Y.; Yu, M.; Zhang, C.; Jiang, L.; Zhang, X.; Zhao, Y. Melting performance
enhancement of phase change material with magnetic particles under rotating magnetic
field. J. Energy Storage 2021, 38, 102540. [CrossRef]  Yang, L.; Ji, W.; Zhang, Z.; Jin, X. Thermal conductivity enhancement of water by
adding graphene Nanosheets: Consideration of particle loading and temperature effects.
Int. Commun. Heat Mass Transf. 2019, 109, 104353. [CrossRef]  Xing, M.; Yu, J.;Wang, R. Experimental study on the thermal conductivity
enhancement of water based nanofluids using different types of carbon nanotubes. Int.
J. Heat Mass Transf. 2015, 88, 609–616. [CrossRef]  Park, S.S.; Kim, N.J. Influence of the oxidation treatment and the average particle
diameter of graphene for thermal conductivity enhancement. J. Ind. Eng. Chem. 2014,
20, 1911–1915. [CrossRef]  Du, C.; Nguyen, Q.; Malekahmadi, O.; Mardani, A.; Jokar, Z.; Babadi, E.; D’Orazio, A.;
Karimipour, A.; Li, Z.; Bach, Q.V. Thermal conductivity enhancement of nanofluid by
adding multiwalled carbon nanotubes: Characterization and numerical modeling
patterns. Math. Methods Appl. Sci. 2020. [CrossRef]  Abhishek, A.; Kumar, B.; Kim, M.H.; Lee, Y.T.; Chung, J.D.; Kim, S.T.; Kim, T.; Lee,
C.; Lee, K. Comparison of the performance of iceoncoil LTES tanks with horizontal and
vertical tubes. Energy Build. 2019, 183, 45–53. [CrossRef]  Vyshak, N.R.; Jilani, G. Numerical analysis of latent heat thermal energy storage
system. Energy Convers. Manag. 2007, 48, 2161–2168. [CrossRef]  Trp, A.; Lenic, K.; Frankovic, B. Analysis of the influence of operating conditions
and geometric parameters on heat transfer in waterparaffin shellandtube latent
thermal energy storage unit. Appl. Therm. Eng. 2006, 26, 1830–1839. [CrossRef]  Soltan, B.K.; Ardehali, M.M. Numerical simulation of water solidification phenomenon
for iceoncoil thermal energy storage application. Energy Convers. Manag. 2003, 44,
85–92. [CrossRef]  Kousha, N.; Hosseini, M.J.; Aligoodarz, M.R.; Pakrouh, R.; Bahrampoury, R. Effect
of inclination angle on the performance of a shell and tube heat storage unit—An
experimental study. Appl. Therm. Eng. 2017, 112, 1497–1509. [CrossRef]  Yuan, Y.; Cao, X.; Xiang, B.; Du, Y. Effect of installation angle of fins on melting
characteristics of annular unit for latent heat thermal energy storage. Sol. Energy 2016,
136, 365–378. [CrossRef]  Liu, Z.; Quan, Z.; Zhao, Y.; Jing, H.; Liu, X.;Wang, L. Experimental research on the
performance of ice thermal energy storage device based on micro heat pipe arrays.
Appl. Therm. Eng. 2021, 185, 116452. [CrossRef]  Chiu, J.N.W.; Martin, V. Submerged finned heat exchanger latent heat storage design
and its experimental verification. Appl. Energy 2012, 93, 507–516. [CrossRef]  Darzi, A.R.; Farhadi, M.; Sedighi, K. Numerical study of melting inside concentric
and eccentric horizontal annulus. Appl. Math. Model. 2012, 36, 4080–4086.
[CrossRef]  Cao, X.; Yuan, Y.; Xiang, B.; Haghighat, F. Effect of natural convection on melting
performance of eccentric horizontal shell and tube latent heat storage unit. Sustain. Cities
Soc. 2018, 38, 571–581. [CrossRef]  Yu, C.; Peng, Q.; Liu, X.; Cao, P.; Yao, F. Role of metal foam on ice storage
performance for a cold thermal energy storage (CTES) system. J. Energy Storage 2020,
28, 101201. [CrossRef]  Lou, X.;Wang, H. Role of copper foam on solidification performance of icecool
storage sphere system. J. Energy Storage 2022, 47, 103552. [CrossRef]  Zhao, Y.; Li, Z.; Utaka, Y.; Chen, Z.; Ohkubo, H. Adhesion characteristics of ice in urea
aqueous solution for efficient slurry formation in cold storage. Int. J. Refrig. 2019, 100,
335–342. [CrossRef]  Wu, F.; Fan, Y.B.; Zhang, X.J.; Zhang, H.; Wang, Z.L.; Wang, Z.W.; Jiang, L.
Performance prediction on ice melting process for cold energy utilization: Effect of
natural convection. J. Energy Storage 2022, 55, 105638. [CrossRef]