Numerical Assessment on Fin Design Parameters Employed for Augmentation of Natural Convection and Fluid Flow in a Horizontal Latent Heat Thermal Energy Storage Unit ()
1. Introduction
In recent years, PCM-LHTES systems have been widely interested, including their system design, optimization processes, phase transient behavior, and heat-transfer enhancement techniques. Storage systems with PCM can be applied using two different approaches: one latent heat storage (LHS) is to use the latent heat to store a great amount of energy in a very small temperature range, while the other sensible heat storage (SHS) is to use the same property but as a temperature regulator [1]. In both cases, either using the advantage of the PCM to store energy or using it as a temperature regulator, it is very important to consider and study each application carefully and deeply, as each application has different requirements and specifications.
If the principle objective is to store energy, LHTES is the most promising and attractive due to its high energy density and ability to store energy at nearly constant temperature corresponding to the phase transition temperature of the material [2]. It is to be noted that, the phase change temperature of the PCM must be between the charging/discharging temperatures (or between the source temperature and the demand one) in order to be able to complete a full cycle.
However, LHTES has a serious drawback: A slow charging and discharging rate caused by the unacceptably low-thermal conductivity of PCMs, which limits heat transfer rate during both charging and discharging process [3]. To fix this problem heat transfer enhancement techniques are required for most latent heat storage systems. In general, these techniques can be summarized as two categories: one is to improve properties of PCMs itself (density, thermal conductivity, specific heat, viscosity) and another is to optimize the heat transfer conditions, including the increase of heat transfer area by introducing fins (for instance).
Precise evaluation of heat transfer from extended heat transfer area of complex structure to PCM due to convection and conduction is a complicated issue but crucial for the design optimization of thermal storage system. In this context, it needs to be clear that heat transfer in PCM storage is a transient, non-linear phenomenon with a moving solid-liquid interface. Therefore, authentic PCMs properties along with accurate numerical methods are required to predict non-linear behavior (natural convection) during PCM melting process [4].
One of the aspect to understand the melting behavior of PCM is the buoyancy-driven natural convection that dominates the heat transfer processes while melting of PCM body and plays a significant role in the transient heat storage capability of PCMs. The physical state of the PCM molten pool depends on the heat transfer at the pool boundaries. The heat transfer from the pool to the boundaries is mostly dominated by natural convection [5]. Natural convection heat transfer in liquid pools with an internal energy source has been the subject of a considerable number of fundamental and applied research programs. That is why, an intensive knowledge of the natural convection heat transfer is essential for predicting PCMs melt pool thermal hydraulics.
Investigations of natural convection phenomena in a fluid with volumetric heat generation began in the early 1970s with Kulacki and collaborators [6] [7] [8] , who conducted several experiments using Joule heating as a volumetric heat source. In those experiments, heat transfer through a horizontal fluid layer was assessed for different boundary cooling arrangements.
Seeniraj et al. [9] investigated transient behavior of high phase change temperature materials stored in heat pipe oriented LHTES unit. They noticed that by using bare tubes some PCM remains in its solid state near to the outlet (Figure 1).
Agyenim et al. [10] considered three horizontal concentric tubes storage units: One unfinned (control PCM unit) and the other two with circular and longitudinal fins, as shown in Figure 2. The authors compared the thermal performance of the three configurations experimentally and the results showed that the system with the longitudinal fins provided the best thermal response during melting.
Lacroix [11] considered a shell-and-tube storage unit with a PCM on the shell side while HTF at the inner tubes and analyzed numerically the transient behavior of the system. Parametric studies were performed and the effect of various parameters and the presence of fins on the thermal behavior of the storage unit taking into account the effect of natural convection of the melting PCM were investigated. He found that the annular fins resulted in a better performance at moderate flow rates and small inlet temperatures of the heat transfer fluid.
Mahmoud Jourabian et al. [12] conducted a numerical investigation of melting phenomenon with natural convection in a cavity with fin by using enthalpy-based lattice Boltzmann method. The obtained results show that the rate of melting increases when the relative thermal conductivity and the length of the fin become greater. It has also found that the variation of vertical position of the fin from bottom to middle has an insignificant effect on melting while it causes the increase of full melting time when the fin is mounted on the top of the cavity.
Yang et al. [13] reported that fins improve the energy storage and release in a LHTES unit and that smaller fin pitches lead to a shorter melting/solidification period mainly due to the smaller PCM volume and larger heat transfer surface.
Mahmoud Jourabian et al. [14] performed a numerical investigation on melting of phase change material (PCM) enhanced by nanoparticles inside a cylindrical tube using the lattice Boltzmann method. They had found that the melting rate is the same for all regions of the cylinder for a low Rayleigh number, while it intensifies at the top half of the cylinder for a moderate Rayleigh number.
![]()
Figure 1. Enhancement of heat transfer with heat pipe in LHTES unit, Seeniraj et al. [9].
![]()
Figure 2. Enhancement of heat transfer with fin oriented LHTES unit, Agyenim et al. [10].
Bergman et al. [15] reported a numerical study on a high temperature shell-and-tube LHTES system integrated with HPs for solar thermal power plant. Both charging and discharging modes were considered during the simulation and the results showed a significant increase in the charging and discharging rates due to the presence of HPs for two cases of HTF flow pattern. The energy stored within the PCM during melting near a single HP is approximately 13% of the energy stored in PCM next to the HTF tube.
Mahmoud Jourabian et al. [16] carried out a numerical analysis on heat transfer and fluid flow processes of natural convection melting of a phase change material inside a partially heated square cavity, where the momentum and energy equations are solved by using enthalpy-based lattice Boltzmann method combined with multi distribution function model. Finally, the dependence of liquid fraction, temperatures of vertical nodes and average Nusselt number on the positions of heated plates is investigated quantitatively.
Doma fiski and Fellah [17] performed exergy analysis to evaluate the thermal performance of a latent heat storage unit with two PCMs in series for charging and discharging cycles. They recommended that the melting temperature of the downstream unit should be nearly equal to the ambient temperature in order to have the best second-law efficiency.
Mahmoud et al. [18] examined a single cavity, parallel fin arrangement, cross fin arrangement and honeycomb cavity arrangement in a PCM-based heat sink units using PCMs of different melting points. They showed that increasing the number of fins could improve the heat distribution into the PCM, resulting in lower heat sinks peak temperatures.
Even though several experiments have been performed in this regard, a numerical analysis is far more than complementary due to the complexity and limitations in measurement capabilities. Furthermore, much research effort has been directed on the assessment of the effects of fins geometry on the heat transfer. However, only a few studies have focused on fin design parameters like effects of number of fins, fin length and fin thickness on thermal response of the storage units.
In the framework of this research work, close-contact melting process of phase change materials (PCMs) inside a horizontal cylindrical capsule with longitudinal fins of different configurations are numerically investigated, which is based on an enthalpy porosity method. The investigation is conducted by ANSYS FLUENT code, which is an efficient numerical analysis tool for investigating convective heat transfer phenomena during PCM melting process. In this work, several simulation cases are performed where two commercial PCMs are utilized in a 2D cylindrical pipe with their thermo-physical properties (Section 4.0) as input for modelling. The selected modelling approach is then validated against experimental data (Section 9.0) to qualify our model to run the proceeding simulations. Finally, parametric study is carried out to evaluate the effect of length, thickness and number of longitudinal fins on the thermal performance of heat storage and recovery as a latent heat thermal energy storage (LHTES) unit dedicated for performing heating load management.
2. Geometry Model
Cylindrical geometries are considered as the most promising for the devices of commercial heat exchangers, because most engineering systems use cylindrical pipes that have high efficiency in a minimum volume. The preliminary design of the Horizontal Concentric Tube Heat Exchanger (in our case) is a cylindrical tube with inner radius 50 mm (Figure 3). The inner and outer radius of the Inner pipe (made of Steel) is set to 10 mm and 12 mm respectively. The outer domain consists of PCM material and Heat Transfer Fluid (HTF)—boiling water is flowing through the Inner tube, which is assumed to have isothermal temperature. The fins aligned longitudinally with Inner pipe would accelerate the heat transfer from heat transfer fluid (HTF) to PCM depends on the initial temperature conditions.
The two commercialized PCM materials: Rubitherm 50 (RT 50) and Climsel C58 are encapsulated in the concentric heat exchanger.
To investigate the impact of number of fins, fin length and thickness, we have considered a test matrix given in Table 1. However, the case with four (4) fins
![]()
Figure 3. Dimensions of the of heat storage unit.
![]()
Table 1. Test matrix with different fin configuration.
having 15 mm length and 5 mm thickness is considered here as Reference case where, RT 50 paraffin is used as PCM.
The cross section of heat storage unit (reference case: RT 50 with 4 fins-15 mm length-5 mm thickness) is presented here as follows (Figure 4).
3. Computation Domain
Considering the 3D geometry and mesh, CFD method is computationally expensive and not affordable for so many cases in sensitivity analysis of the key parameters. Therefore, a simplified model (2D case) is used to simulate the melting process.
The 2-D model has been developed by state of the art software package ANSYS ICEM and the corresponding computational mesh is generated by ICEM CFD using O-grid method, which guaranteed the unstructured mesh having hexahedral cells. The mesh near the boundaries was refined as shown in Figure 5. The total number of cells and nodes are 76,112 and 76,648 respectively. The maximum cell size of meshes is below 1 mm. The simulations are run on a LENOVO D20 workstation with two Xeon(R) E5620 @ 2.4 GHz CPUs and 20 GB RAM.
4. PCM Properties and Boundary Conditions
Two different PCMs are considered here, identified as paraffin Rubitherm (RT 50) and Climsel (C58). The advantage of paraffin RT 50 as a latent heat energy storage material is that it is chemically stable, non-poisonous and non-corrosive over a large storage period, long life product, with stable performance through the phase change cycles. Compared to RT50, C58 is a commercial PCM product based on Sodium Acetates. It is formulated with nucleating and gelling agents to prevent super cooling and phase segregation effects. In addition, it possesses higher thermal conductivity and latent heat of fusion than most organic commercial products having similar phase change temperatures.
![]()
Figure 4. Cross section of heat storage unit.
![]()
Figure 5. Cross section of the 2-D cylindrical mesh.
The PCM has different thermos-physical properties but similar melting/solidification temperature; therefore, they could be suitable for the same operational temperature in a heating system. The following table shows the PCM thermo-physical properties of the used PCM as follows.
Boundary Conditions Applied
In order to define PCM properties a first assumption has been made considering a single intermediate value for the melting and solidification temperature. For Climsel C-58, a solidification temperature equal to 54˚C (327 K) and a melting temperature of 57˚C (330 K) have been considered. Similarly, for RT the solidification temperature is 48˚C (321 K) while the melting temperature is 53˚C (326 K). Then the variation in density and conductivity of PCMs with the temperature is considered as a linear relationship as shown in Figure 6 & Figure 7.
The thermal boundary and initial conditions are linked to the energy equation. The initial temperature in the whole domain is set to 42˚C (315 K) and an isothermal boundary condition is applied to the inner pipe (373 K) while the outer pipe has been modelled as perfectly insulated (null heat flux) (Figure 8).
5. Numerical Method
5.1. Governing Equations
In order to simulate phase change of melting in a concentric heat exchanger with longitudinal fins, enthalpy-porosity method [19] [20] is used. The flow is considered unsteady, laminar and incompressible. The viscous dissipation term is considered negligible, so that the viscous incompressible flow and the temperature distribution in annulus space are described by the Navier–Stokes and thermal energy equations, respectively. Consequently, the continuity, momentum, and thermal energy equations can be expressed as follows:
Continuity:
Thermal Energy
![]()
Figure 6. Density and conductivity versus temperature (Climsel C58).
![]()
Figure 7. Density versus temperature (Rubitherm-50 Paraffin).
![]()
Figure 8. Temperature initial and boundary conditions.
The enthalpy of the material is computed as the sum of the sensible enthalpy h, and the latent heat, ΔH:
where,
The latent heat content can be written in terms of the latent heat of the material L:
ΔH = λL, where ΔH may vary from zero (solid) to L (liquid). Therefore, the liquid fraction, λ, can be defined as:
if T < TSolidus.
if T < TLiquidus.
if TSolidus < T < TLiquidus.
is the Darcy’s law damping terms (as source term) that are added to the momentum equation due to phase change effect on convection. It is defined as:
The coefficient Amush is a mushy zone constant. This constant is a large number, usually 104 - 107. In present study, Amush is assumed constant and is set to 106.
5.2. Numerical Procedure
The SIMPLE algorithm is utilized for solving the governing equations. A second order upwind scheme has been considered for the momentum and energy equations, as this scheme is largely more stable than a first order scheme. By solving the governing equations at each time step, liquid mass fraction has been updated. The grid size and the time step are chosen after careful examination of the independency of the results. The convergence is checked at each time step, with the convergence criterion of 107 for all variables
6. Analysis of Reference Case
In order to study thermal performance of our concentric heat exchanger, a case study has been performed numerically which is treated as a reference case. As a part of the reference case, RT-50 has been considered as PCM, encapsulated in a cylindrical tube longitudinally connected with 4 fins having 15 mm length and 5 mm thickness from the inner side.
Thus, starting from the inner cylinder through which HTF (boiling water, 373 K) is flowing, the PCM absorbs heat and when the temperature reaches to the phase transition, the melting process began. Actually there are two regions exist during the charging processes: 1) solid region and 2) molten region. At the beginning of the melting process, conduction heat transfer was dominant due to the close contact of solid PCM and the outer surface of the inner tube and fins and because of conduction effect, a thin liquid layer develops all around the contact region and solid PCM receives heat from the melted section by convection.
Then the convection mechanism of the heat transfer drives the recirculation inside the molten region, which is due to the buoyancy forces induced by the density gradients as a result of the temperature differences. Figure 9 shows sharp rise in temperature takes place at location 0˚, located at the upper most section of the storage unit, because of the effect of buoyancy.
This inner recirculation enhances the heat transfer within the melted PCM region, which can be explained by the fact that the points near the upper section
![]()
Figure 9. Density variation at different time steps.
reach higher temperatures than the lower section. As a result, the solid PCM sank to the bottom of the cavity due to its higher density and gravity with respect to the liquid PCM. Then, the liquid layer thickness grows up following cylindrical shape (Figure 9) and the phase transition region has an apparent thickness which is referred to as “mushy region”.
As can be seen from the image (Figure 9), ten minutes after it is detected that the melting front is no more symmetrical on the x-axis; the melting in the upper part of the heat exchanger is facilitated by the natural convection. With the PCM density decreases, the liquid part tends to rise up inside the heat exchanger while the solid part tends downwards. Therefore, the melting process becomes much more visible at the top of the heat exchanger while melting at the bottom of the pipe is found to be inefficient. After 30 minutes, the most of the upper part is melted while the lower part is still solid. After 60 minutes, only the bottom part is still solid; once this happens the melting process slows down. Indeed, the bottom region cannot enjoy the heat transfer improvement brought by the natural convection because the liquid-status of PCM tend to ascend and induce no fluid motion down to the bottom.
It is clear from the series of images that during the first phase, relatively fast melting is achieved while during the second phase the rate of melting is comparatively slower.
Furthermore, as an index of the heat exchange between the HTF and PCM, the Nusselt number at the inner pope surface has been considered. It can be seen that during the first phase a Nu around 50 is achieved with small oscillation due to the effect of natural convection and fluid movement around the pipe. Then in the second phase, the Nu strongly decreases and it becomes constant at around 1 (Figure 10).
![]()
Figure 10. Nusselt number at inner pipe PCM interface vs time (Sec).
7. Mesh Independence Test
Meshing grid size also has a significant impact on simulation accuracy. In practical, the simulation of melting process in LHTES is sensitive to the quality and local features of meshes, and the sharp grids at local region would lead to the divergence of calculations. In general, the transient simulation should satisfy the Courant-Friedrichs-Lewy conditions in case of divergence, that is:
where, u is the magnitude of velocity,
is the time-step,
are the length interval in the respect of x, y direction. The preliminary calculations show that the maximum velocity in our case would not exceed
, and the time steps in the assignment are given as an interval of 0.01 - 0.05 Sec that comes out the optimal length interval would be in the range of greater than 0.35 mm with considering the largest time step and default courant number of 2.0. In contrast, the large length interval would result in the lower accuracy of simulation.
In this section, we considered several meshes with different densities with considering the compromise of accuracy and convergence, and the following table gives the nodes number of meshes, the boundary and initial conditions for the case setups.
The following figures (Figures 11-13) present the variation of density, liquid fraction and total energy over time.
Obviously, the length interval relates to the courant number, which affects to the convergence of calculation. Therefore, with increasing nodes quantity, the courant number increase as well, that might have the possibility of leading to the instability of calculation (Figure 14).
While comparing case 1, 2 and 3, the divergence is observed in case 2 and case 4 due the excess of Courant number. With considering the computational time and limitation of Courant number, 22,936 nodes number can satisfy the requirement of convergence and accuracy for the case of 4 fins of 15 mm for the times step of 0.01 - 0.05s.
![]()
Figure 14. Fluid flow due to Buoyancy driven force.
8. Time-Step Independence Study
In order to capture all of the flow parametric fluctuations during simulations, it is needed to use small time step sizes. Small time steps can make sure the local courant number satisfies the requirement of convergence (C < 2.0 for 2-D model), but consequently it becomes computationally more expensive. Therefore, the time step independence study is necessary to be carried out to choose the optimal value. Here, the setup for the cases with different time steps is shown below.
The variations of density, liquid fraction and total energy with flow time of are shown below in Figures 15-17.
It is observed that the full melting is reached at 70.42, 71.62 and 75.71 min with the time steps 0.01, 0.03 and 0.05 s respectively and the variations of liquid fraction and density with that time steps are very close to each other. The significant variation is found from the first 5 to 15 minutes and then the variation becomes narrow until it fully melt. Comparing among the time steps 0.05, 0.03 and 0.01, it is found that both the time steps (0.01 and 0.03) are acceptable to satisfy the calculation requirement.
The variation of total energy is slightly different with each other and the final energy storage capacities are 230.224, 231.86 and 232.13 kJ/kg for the cases of time step 0.01, 0.03 and 0.05 s, respectively, meaning that the change of time step has a limited influence to the final melting time and the heat storage capacity.
9. Validation
As a part of validation procedure, we have chosen experimental data of total enthalpy change for a cylindrical capsule filled with C58 PCM in the Department of Energy at KTH. The PCM cylinder was vertically placed in a cylindrical storage tank, which is tested by its thermal performance of heat storage and recovery as a latent heat thermal energy storage (LHTES) unit dedicated for performing heating load management.
Figure 18 illustrates an experimental set-up has been designed and built in order to investigate the melting process of a latent heat storage system. The tank contains 50 PCM cylinders in total. It ensures that downward laminar flow of HTF, which is evenly distributed by diffusers before contacting the PCM capsules. HTF was pumped from the top throughout the entire unit (to exchange heat with the PCM) during charging. In the testing condition, HTF temperature
![]()
Figure 18. The illustration of experimental facilities for model validation (a) The schematic of PCM cylinder placement in the storage tank; (b) The geometries of the PCM cylinder; (c) The appearance of the insulated storage tank loaded with PCM cylinders.
is maintained constantly at 65˚C, while the initial temperature of PCMs was kept homogenously at 45˚C. The measuring system consisted of an array of calibrated T-type thermocouples with a mean deviation of 0.05˚C from readings of a precision digital thermometer, a volumetric flow meter with a relative measurement inaccuracy of 0.5% of its reading. A data acquisition system is employed to collect these experimental recordings (Figure 19).
Boundary condition applied:
Top and side of the PCM cylinder: Adiabatic (as per experimental condition).
Outer surface of PCM cylinder: Isothermal boundary condition (338 K).
Initial temperature 318 K (Figure 20).
It shows that there is a good agreement between simulation and experimental result with some minor deviation.
The most probable reason behind the deviation is that constant material properties of the PCM used in the simulation. If the temperature dependent material properties were known, the numerical methods would have shown more precise results for the temperature of the PCM. Thus, the material properties of the PCM should be well known in order to obtain sufficiently accurate results with the numerical methods. Another reason can be numerical errors comes from the numerical code.
![]()
Figure 19. Cross section of the PCM Cylinder mesh.
![]()
Figure 20. Comparison between experiment and simulation result.
10. Results & Discussion
Parametric analysis
After the reference case (number of fin 4 and fin length 15 mm) is done, we now examine the impact of extended surface by increasing the numbers of fins (6 and 8), fin length and fin thickness and comparing heat transfer characteristics in terms of melting time, melt fraction and total energy storage time among these cases.
The results reveals that higher fin number & length, quicker heat transfer rate and less time required to melt and store the same amount of energy in the capsule as shown in Figures 21-23. Although the thermal conductivity of paraffin PCM is lower than that of Climsel C 58, however because of the high influence of latent heat of fusion, RT 50 paraffin (160,000 J/Kg) reaches to the melt fraction 1.0 quicker than Climsel C-58 (212,000 J/Kg). In addition, RT50 exhibits lower viscosity than C58 because C58 is gelled PCM to prevent phase segregation after phase changes. Therefore, less fluid motion is induced in molten C58 by buoyancy effects, which could also lead to the slower melting rates relative to RT50, as depicted in Figure 21.
That is why, RT 50 Paraffin takes almost 50% less time than Climsel C58 to reach melt fraction 1 and to be fully charged as shown in Figure 21. That means the fin number and length definitely influence the thermal behavior of the PCM melting process due to improved natural convection and heat conduction. However, fin thickness shows relatively lesser effects on the melting rate as shown in Figure 24, meaning that the effect of fin thickness is smaller as compared to that of the fin length and number.
Calculation of Total Heat Stored in the Capsule (at Melt Fraction 1)
*Initial energy supplied to paraffin RT 50 PCM is 33 KJ/Kg.
![]()
Figure 21. Variation of liquid fraction with time.
Final energy stored in paraffin PCM: 253 KJ/Kg (6 fins) & 260 KJ/Kg (8 fins).
So, total energy stored in the Paraffin PCM is 220 KJ/Kg (6 Fins) and 227 KJ/Kg (8 Fins).
*Initial energy supplied to Climsel C58 PCM is 61 KJ/Kg.
Final energy stored in Climsel C58 PCM: 410 KJ/Kg (6 fins) & 425 KJ/Kg (8 fins).
So, the total energy stored in Climsel C58 PCM: 349 KJ/Kg (6 fins) & 364 KJ/Kg (8 fins).
Based on the calculation, it is also proven that the number of fins plays a vital role in case of charging process.
11. Conclusions & Future Recommendations
The present study is intended to provide insights of thermal performance of a horizontal concentric heat exchanger, which is numerically investigated to evaluate the heat transfer enhancement process by adding fins with different configurations. As a part of this investigation, the melting process is simulated from the onset of phase change to the offset involving physics of natural convection in PCM fluid pool. The conclusions based on the numerical simulations are summarized as follows:
a) The transient melting process of the PCM depends on the thermos-physical properties and design parameters of the LHTES unit.
b) Slower melting rate is found in case of C58 compare to RT 50 because C58 is gelled PCM to prevent phase segregation after phase changes; meaning C58 exhibits higher viscosity than RT50 and that is why less fluid motion is induced in molten C58 by buoyancy effects, which could also lead to the slower melting rate.
c) Melting and solidification rates are largely influenced by natural convection and conduction.
d) Increasing the number of fins and length increases the overall thermal performance of the heat sink by faster energy extraction from fins to the PCM whereas; increasing fin thickness provided relatively lesser improvement.
e) The sharp rise in temperatures takes place at the uppermost section of the store unit, because of the buoyancy effects.
Future CFD analysis can be conducted with other CFD code like CFX, OPENFOAM, STAR CCM+ or CFD++ and compared with the current code ANSYS Fluent to justify their applicability in the area of study.
Apart from this, the use of extended surfaces (by fins, heat pipe) with multiple layer of PCM having higher conductivity materials by impregnation of porous materials into the base PCM is recommended for future research.
Nomenclature
Cp: Heat capacity [J∙kg−1∙K−1]
H: Height of the melt pool
Pr: Prandtl number (υ/α)
Ra: Rayleigh number (g∙βH5Q/k∙α∙υ)
Greek Symbols
Β: Thermal expansion coefficient [K−1]
ρ: Density [kg∙m−3]
α: Thermal diffusivity [m2∙S−1]
υ: Kinetic viscosity [m2∙S−1]
g: Acceleration due to gravity force [m∙s−2]
k: Heat conductivity [Wm−1∙K−1]
Subscripts
Fig: Figure
Ref: Reference