Thermal Analysis of Pure and Nanoparticle-Enhanced PCM—Application in Concentric Tube Heat Exchanger

In this paper, organic phase change materials are modeled and studied both numerically and computationally. The results are compared with those of other research available in the international literature. Regardless of their heat storing capacity, the low heat conductivity of PCMs constitutes a significant obstacle to the further advancement of technologies pertaining to latent heat energy storage, as it hampers the immediate response of systems both at heat storage and at heat recovery. In this work, two types of nanoparticles are used as enhancing media, and the calculations are repeated in terms of melting front and stored energy enhancement. The results show that a 6.5% value for the melt percentage enhancement and a 5.5% value for the heat stored enhancement are exhibited when copper nanoparticles are utilized.


Introduction
The constantly increasing emissions of greenhouse gases as well as numerous other gases toxic for human organisms necessitate urgent measures in the direction of energy saving and reducing fossil fuel consumption. Significant contributions to achieving both objectives can be expected by improving the energy performance of buildings to secure lower energy consumption for heating and cooling, and by efficiently utilizing the power generated by means of the thermal storage of untapped heat. In addition, the more efficient utilization of renewable energy sources (RES) shall result in an increase in the share of RES in gross final energy consumption and, consequently, in reduced CO 2 emissions.
Among diverse technologies, phase change materials (PCMs) offer considerable advantages, rendering them ideal for thermal energy storage as well as for preserving constant temperatures for longer time spans. PCMs capitalize on the latent heat of fusion and, therefore, store more thermal energy than sensible heat materials. Moreover, during the phase change, PCMs maintain an almost stable temperature, thus contributing to the reduction of the required heat and cooling loads of the buildings [1]. Despite the aforementioned advantages, PCMs require longer time spans to change from the solid to the liquid state and vice versa, due to their low thermal conductivity coefficient. The complexity of the equations describing the phenomenon of the heat flow through the PCM (as a result of the phase change process and the random movement of the fusion-solidification front) does not lend itself to the formation of analytic relations for the calculation of the temperature distribution in the PCM and the heat flow. Consequently, for the purposes of the mathematical analysis and the description of the phenomenon, both the numerical solution and mathematical simulation of the equations are required as well as the experimental calculation of the relevant quantities. Several

The Examined Domain
The examined domain is depicted in Figure 1. This work regards the analysis of two different geometries. Therefore, Figure 1a shows a rectangular geometry of length 50 mm, which consists of three adiabatic sides and a fourth one subject to a heat source. Figure 1b illustrates the domain used in the concentric tube heat exchanger application in which, apart from the pure PCM analysis, nanoparticles are added in order to increase the low conductivity of the paraffins. The dimensions of the inner and outer tubes are 22 mm and 125 mm, respectively.

Pure PCM Mathematical Modeling
The problem of the two-dimensional study, as opposed to the one-dimensional one, requires the solution of all Navier-Stokes equations due to the movement of the liquid PCM caused by the temperature variation in the material. The movement of the PCM facilitates heat flow by natural convection, resulting in the faster melting and distortion of the phase change front. Presented below is the system of equations solving the phase change phenomenon in the material in Cartesian coordinates (Equations (1)-(4)) and cylindrical coordinates (Equations (5)- (8)) [19,20].

Pure PCM Mathematical Modeling
The problem of the two-dimensional study, as opposed to the one-dimensional one, requires the solution of all Navier-Stokes equations due to the movement of the liquid PCM caused by the temperature variation in the material. The movement of the PCM facilitates heat flow by natural convection, resulting in the faster melting and distortion of the phase change front. Presented below is the system of equations solving the phase change phenomenon in the material in Cartesian coordinates (Equations (1)-(4)) and cylindrical coordinates (Equations (5)-(8)) [19,20].
In order to zeroize the velocity at the points where the PCM is entirely solid [21], the term S is inserted in the equations for the conservation of momentum, which is proportional to the melted material ratio and given by the following expression: Additionally, for the calculation of the pressure field within the liquid PCM, another equation is used [21] to correct the initial hypothetical pressure at any point in time in Cartesian and cylindrical coordinates.
As illustrated by the equations above, density is calculated using the Boussinesq approximation. This method offers satisfactory results for incompressible flows whose other quantities are stable or temperature-dependent variables. In this particular method, the material density is considered to be stable and a thermal expansion constant is used, β (1/K) = 0.001. This method is used to model natural convection provided that the temperature difference is adequately low, i.e.,: The approximation formula is: The operating density (ρ o ) is separate for each material, and the operating temperature (T o ) is set at 50 • C. In cases when there is great temperature difference or the fluid is compressible, a different correlation between the quantities is required, such as a function correlating temperature with density as a result of experimental data [22]. The study of enhanced phase change materials requires the appropriate modification of thermal properties depending on the PCM content in nanoparticles. The density, thermal conductivity, viscosity and heat capacity of the nano-enhanced organic PCM studied above are calculated based on the following expressions [23,24]. The nanoparticles used for analysis are Cu and Al 2 O 3 (Table 1), and the organic PCM is RT50 (Table 2). It should be stressed that for the purposes of the present study, constant nanoparticle properties have been considered while the effect of their size is negligible. Therefore, in the case of nanoparticle-enhanced PCM, the governing equations are the ones described above as well as the following.
where index p refers to nanoparticles and ϕ is the volume concentration of the enhanced PCM. The method of coupling between the momentum and continuity is of vital importance for both the solution process stability and the total computational time required. As to the problem in question, there was difficulty observed in reaching the continuity equation convergence required, which resulted in either multiple iterations per time step or even divergence among iterations. In general, numerical methods are divided into segregated and coupled methods. In the first type, equations are solved successively, while in the second type, equations are solved simultaneously. The successive solution of equations utilizes some default values or the solution products of the previous cycle and solves each equation in sequence, feeding the results of the solution to the next equation. At the end, each equation is controlled for convergence. In the simultaneous solution of equations, all equations are solved in parallel. The first method has lower memory requirements, since one equation is stored at a time, while for the second method, all equations of the problem need to be registered. However, simultaneous Energies 2020, 13, 3841 6 of 19 solution tends to be faster. For the solution of the equation system above, a finite-difference time-domain scheme was used, which was solved by means of code developed in Matlab. One basic assumption for the solution of clipped equations is a stable density value for all terms with the exception of the force term in the momentum equations. The solution procedure is illustrated in the following logical diagram ( Figure 2).
Energies 2020, 13, x FOR PEER REVIEW 6 of 21 need to be registered. However, simultaneous solution tends to be faster. For the solution of the equation system above, a finite-difference time-domain scheme was used, which was solved by means of code developed in Matlab. One basic assumption for the solution of clipped equations is a stable density value for all terms with the exception of the force term in the momentum equations. The solution procedure is illustrated in the following logical diagram ( Figure 2). For the purposes of modeling the phase change phenomenon, a modification of the thermophysical properties of the PCM is required [25,26]. The equations calculating the heat capacity, the density and the thermal conductivity coefficient of the material are given by the equations presented below: where ΔΤ is the temperature range of the phase transition and , are the extreme temperatures of the phase transition range. Finally, Table 3 cites the boundary conditions used in the solution in Cartesian and cylindrical coordinates.  For the purposes of modeling the phase change phenomenon, a modification of the thermophysical properties of the PCM is required [25,26]. The equations calculating the heat capacity, the density and the thermal conductivity coefficient of the material are given by the equations presented below: where ∆T is the temperature range of the phase transition and T s , T l are the extreme temperatures of the phase transition range. Finally, Table 3 cites the boundary conditions used in the solution in Cartesian and cylindrical coordinates. Table 3. Initial and boundary conditions.

Cartesian Coordinates
Cylindrical Coordinates

Computational Solution Scheme
Further to the numerical solution of the Navier-Stokes equations, which describes the phase change phenomenon, simulations of the phase change phenomenon were also performed by means of the COMSOL and Fluent software. For the purposes of the solution for the phenomenon, appropriate modifications to the thermal properties were made to separate the liquid from the solid phases [26,27].
As is generally required in computational methods, within the surface or domain to be studied, computational nodes are specified by means of various methods, thus determining a finite number of computational cells, for each of which a series of calculations is performed. The formation of the appropriate mesh is particularly significant for the proper simulation of the behavior of the various models and phenomena to be studied in each case. In general, the denser the mesh is, the better the approximation but with higher demands on computational resources (CPU time and memory). Therefore, standard practice provides for the determination of meshes that are denser near the critical regions (complex geometries, orifices, boundary layers and plumes) and thinner in areas expected to be of less interest. In the particular case studied, the natural convection occurring due to the temperature difference in the paraffin melt highlights the importance of the mesh for the results. The mesh was defined in terms of triangle cells and was structured using the patch-dependent method. According to this method, the mesh is determined by considering the number of nodes defined at the edge-curves and the maximum and minimum distances defined between the nodes [28].

Results and Discussion
The PCM selected for the two-dimensional analysis of the phase change was paraffin RT50, which was also used for the one-dimensional analysis. A constant time step at 0.1 s was defined for the solution of the problem. For the completion of each time step and the attainment of convergence, constantly increasing iterations are required. This is due to the fact that as PCM melting proceeds, more and more computational nodes or cells (in the case of COMSOL) are affected by natural convection and acquire the speed and properties of fluid, thus leading to an increase in residuals and, consequently, in the number of iterations required. Therefore, the computational time corresponding to each real second increases and the model solution decelerates. It is worth mentioning characteristically that while, initially, the iterations required are 20-30 per time step, after 900 s and 9000 steps, the iterations required are over 100 per iteration. The total time required for each of the cases so as to reach 30 min in the evolution of the phenomenon is over 5 days. Moreover, in the case of the pure PCM, in which the kinematic connectivity is lower and higher speeds are developed, the convergence of the continuity equation is even harder to attain and the total iterations per time step further increase.

Grid Independence
A solution by a numerical or computational method is controlled in terms of validity by the grid independence in the computational domain. In the case of the solution code in Matlab, the temperature variation at one particular position within the material is used for selecting the density of the grid. Thus, in the case of the numerical solution, an 81 × 81 grid is selected as shown in Figure 3.
As far as the computational solution by means of software is concerned, the effective simulation of models in CFD software (Multiphysics 5.2) is largely dependent on the mesh density and the time step. Further to setting these parameters with a view to calculation convergence, it is also imperative that the grid density and the size of the time step do not affect, within certain boundaries, the results of the model solution. For this particular reason, prior to the final solution of the model, through trials, with increasingly denser grids, the number of nodes is defined, a further increase in which does not significantly affect the results [28]. Figure 4 presents the melt percentage as function of time for the different grids as well as the evolution over time of the percentage variation from the melt percentage with the highest density. The grids studied exhibit densities of 70,000 and 100,000 computational nodes, respectively. It can be observed that out of the two grids, the one with the fewest nodes exhibits a significant difference in relation to the dense grid, which actually increases over time. On the other hand, the 100,000-node grid exhibits less difference in relation to the dense grid, which actually remains relatively stable over time. Based on the above, the 100,000-node grid was selected as being both accurate and computationally cheaper.
Energies 2020, 13, x FOR PEER REVIEW 8 of 21 As far as the computational solution by means of software is concerned, the effective simulation of models in CFD software (Multiphysics 5.2) is largely dependent on the mesh density and the time step. Further to setting these parameters with a view to calculation convergence, it is also imperative that the grid density and the size of the time step do not affect, within certain boundaries, the results of the model solution. For this particular reason, prior to the final solution of the model, through trials, with increasingly denser grids, the number of nodes is defined, a further increase in which does not significantly affect the results [28]. Figure 4 presents the melt percentage as function of time for the different grids as well as the evolution over time of the percentage variation from the melt percentage with the highest density. The grids studied exhibit densities of 70,000 and 100,000 computational nodes, respectively. It can be observed that out of the two grids, the one with the fewest nodes exhibits a significant difference in relation to the dense grid, which actually increases over time. On the other hand, the 100,000-node grid exhibits less difference in relation to the dense grid, which actually remains relatively stable over time. Based on the above, the 100,000-node grid was selected as being both accurate and computationally cheaper.     least temporal step that can be controlled in such problems; thus, it is used as a reference in our approach. Through this step, Figure 5 was created, as it depicts the deviation of the solution from the reference one for other time steps. The discrepancy is minimized for time step 0.1 s, which was finally selected. As for the selection of the grid and the time step, the cost in terms of computational Energies 2020, 13, 3841 9 of 19 Figure 5 presents the evolution of the melt percentage over time for different time steps and their difference from the smallest time step controlled, which was 0.01 s. The time step of 0.01 s is the least temporal step that can be controlled in such problems; thus, it is used as a reference in our approach. Through this step, Figure 5 was created, as it depicts the deviation of the solution from the reference one for other time steps. The discrepancy is minimized for time step 0.1 s, which was finally selected. As for the selection of the grid and the time step, the cost in terms of computational resources is also considered. Melting fraction for δt=0.01s

Rectangular Cross Section Analysis
The geometry selected for the accuracy of the computational model is rectangular, which is the most common geometry for PCM study. The dimensions of the cross section are 5 cm × 5 cm, while the field depth is not defined. The validation of the COMSOL model is accomplished by means of comparison with the available literature. The simulation parameters are adopted by the respective work in Souyafane et al. [29], and the COMSOL model in which the modified properties are inserted runs until 5000 s. The complexity of the problem is illustrated by the fact that in order for the model to store as a result only the melt front up to the time of 5000 s, 98 h of continuous simulation was required. Figure 6 represents the melting front as isothermal curves and includes the comparison with the literature results. As observed, there is agreement between the two methods, which renders the modified COMSOL model reliable.

Rectangular Cross Section Analysis
The geometry selected for the accuracy of the computational model is rectangular, which is the most common geometry for PCM study. The dimensions of the cross section are 5 cm × 5 cm, while the field depth is not defined. The validation of the COMSOL model is accomplished by means of comparison with the available literature. The simulation parameters are adopted by the respective work in Souyafane et al. [29], and the COMSOL model in which the modified properties are inserted runs until 5000 s. The complexity of the problem is illustrated by the fact that in order for the model to store as a result only the melt front up to the time of 5000 s, 98 h of continuous simulation was required. Figure 6 represents the melting front as isothermal curves and includes the comparison with the literature results. As observed, there is agreement between the two methods, which renders the modified COMSOL model reliable.
Following the validation of the COMSOL model, the simulation parameters are inserted in the Matlab mathematical model. The calculations cover the first 250 s and have a total duration of 16.8 h, which is indicative of the problem complexity. Figure 7 presents the melt front for material RT50 of the geometry selected. The buoyancy effect is obvious and deforms the melt front, creating an area in the northwestern corner, which moves at higher velocity than the rest of the front. The melting front is in complete agreement with the available literature, which renders the solution of the Matlab developed model acceptable. In the case in which the effect of gravity through the buoyancy flow is not taken into consideration, the front moves in uniform flow and the one-dimensional model would suffice for the analysis of the phase change materials. Buoyancy creates a velocity field, the maximum value of which is 10 −3 m/s. In order to illustrate the evolution of the phenomenon, the mathematical model is again used up to the time of 7200 s with a thinner grid so as to slightly reduce the time required to reach the solution. Figure 7b represents the melt front two hours after the onset of the process. Due to natural convection, the melt PCM flows at the upper side of the cross section, thus accelerating non-uniform melting. The total duration of the simulation for the material studied up to the specific time is 160.6 h.
Following the validation of the COMSOL model, the simulation parameters are inserted in the Matlab mathematical model. The calculations cover the first 250 s and have a total duration of 16.8 h, which is indicative of the problem complexity. Figure 7 presents the melt front for material RT50 of the geometry selected. The buoyancy effect is obvious and deforms the melt front, creating an area in the northwestern corner, which moves at higher velocity than the rest of the front. The melting front is in complete agreement with the available literature, which renders the solution of the Matlab developed model acceptable. In the case in which the effect of gravity through the buoyancy flow is not taken into consideration, the front moves in uniform flow and the one-dimensional model would suffice for the analysis of the phase change materials. Buoyancy creates a velocity field, the maximum value of which is 10 −3 m/s. In order to illustrate the evolution of the phenomenon, the mathematical model is again used up to the time of 7200 s with a thinner grid so as to slightly reduce the time required to reach the solution. Figure 7b represents the melt front two hours after the onset of the process. Due to natural convection, the melt PCM flows at the upper side of the cross section, thus accelerating non-uniform melting. The total duration of the simulation for the material studied up to the specific time is 160.6 h.  (a) (b) Figure 6. COMSOL results (a) and comparison with available data from literature (b) [29].
Following the validation of the COMSOL model, the simulation parameters are inserted in the Matlab mathematical model. The calculations cover the first 250 s and have a total duration of 16.8 h, which is indicative of the problem complexity. Figure 7 presents the melt front for material RT50 of the geometry selected. The buoyancy effect is obvious and deforms the melt front, creating an area in the northwestern corner, which moves at higher velocity than the rest of the front. The melting front is in complete agreement with the available literature, which renders the solution of the Matlab developed model acceptable. In the case in which the effect of gravity through the buoyancy flow is not taken into consideration, the front moves in uniform flow and the one-dimensional model would suffice for the analysis of the phase change materials. Buoyancy creates a velocity field, the maximum value of which is 10 −3 m/s. In order to illustrate the evolution of the phenomenon, the mathematical model is again used up to the time of 7200 s with a thinner grid so as to slightly reduce the time required to reach the solution. Figure 7b represents the melt front two hours after the onset of the process. Due to natural convection, the melt PCM flows at the upper side of the cross section, thus accelerating non-uniform melting. The total duration of the simulation for the material studied up to the specific time is 160.6 h. The validity of the developed numerical model is confirmed through comparison with the respective results calculated by the COMSOL model in Figure 8. The deviation observed in the first seconds and for the low values of the melt percentage is due to the density of the grid selected for the numerical scheme in Matlab. The rapid front expansion for the areas near the heat source can be easily observed. Then, the effect of the heat source diminishes and the front evolves at a low rate. For The validity of the developed numerical model is confirmed through comparison with the respective results calculated by the COMSOL model in Figure 8. The deviation observed in the first seconds and for the low values of the melt percentage is due to the density of the grid selected for the numerical scheme in Matlab. The rapid front expansion for the areas near the heat source can be easily observed. Then, the effect of the heat source diminishes and the front evolves at a low rate. For the same time span, the heat stored in the PCM is also represented. The heat is determined by the PCM surface integral, and the results are represented by field depth. The relatively low values of heat stored as shown in Figure 9 are justified by the fact that they correspond to the heat charging of the material with sensible heat. The discrepancy between the numerical solution and the COMSOL solution is due to the reasons already stated for the case of the melt percentage. the same time span, the heat stored in the PCM is also represented. The heat is determined by the PCM surface integral, and the results are represented by field depth. The relatively low values of heat stored as shown in Figure 9 are justified by the fact that they correspond to the heat charging of the material with sensible heat. The discrepancy between the numerical solution and the COMSOL solution is due to the reasons already stated for the case of the melt percentage.  the same time span, the heat stored in the PCM is also represented. The heat is determined by the PCM surface integral, and the results are represented by field depth. The relatively low values of heat stored as shown in Figure 9 are justified by the fact that they correspond to the heat charging of the material with sensible heat. The discrepancy between the numerical solution and the COMSOL solution is due to the reasons already stated for the case of the melt percentage.

Pure PCM Results
The geometry studied consists of a concentric tube heat exchanger with tube diameters of Φ125 and Φ22, respectively. The phase change material RT50 is located in between the concentric tubes while the high temperature (80 • C) is applied to the internal tube. The analysis was conducted in the Fluent environment with the modification of properties similar to those effectuated in COMSOL. The following figures represent the temporal development of the phenomenon and, in particular, the evolution of the melt over the thermal field and the velocity fields for different points in time. Figure 10 presents the variation in the paraffin temperature up to 800 s, while Figure 11 presents the isothermal curves for the whole domain at 900 s. Up to the first 100 s, conduction is the driving force; therefore, the isothermal curves are uniform. After that, gravitational forces and buoyance are starting to affect the phase change and heat is transferred from the source to the material, not only through conduction but via convection as well. Small vortices are created, which enhance the natural convection phenomenon. The shape of the isothermal curve is also typical and refers to the shape of the melt as a result of the buoyancy and the vortices produced by the velocity field.
Fluent environment with the modification of properties similar to those effectuated in COMSOL. The following figures represent the temporal development of the phenomenon and, in particular, the evolution of the melt over the thermal field and the velocity fields for different points in time. Figure 10 presents the variation in the paraffin temperature up to 800 s, while Figure 11 presents the isothermal curves for the whole domain at 900 s. Up to the first 100 s, conduction is the driving force; therefore, the isothermal curves are uniform. After that, gravitational forces and buoyance are starting to affect the phase change and heat is transferred from the source to the material, not only through conduction but via convection as well. Small vortices are created, which enhance the natural convection phenomenon. The shape of the isothermal curve is also typical and refers to the shape of the melt as a result of the buoyancy and the vortices produced by the velocity field.  Farhadi et al. [30] studied the melting of organic PCM with a two-dimensional model of the same geometry as the one analyzed here. The melting formulation of this work is in agreement with the results of Farhadi.
The following figures represent the temporal evolution of the velocity field as well as magnified areas from the field at the point 900 s. Figure 12 illustrates the mode in which velocities are developed as a result of natural convection in the melt as well as the significant differentiation of the measured velocity depending on the area and given the specific scale. In every case, the measured velocity acquires low values-typical of the buoyancy flow-which vary from 10 −5 m/s in the first −3 Figure 11. Temperature field at 900 s.
Farhadi et al. [30] studied the melting of organic PCM with a two-dimensional model of the same geometry as the one analyzed here. The melting formulation of this work is in agreement with the results of Farhadi. The following figures represent the temporal evolution of the velocity field as well as magnified areas from the field at the point 900 s. Figure 12 illustrates the mode in which velocities are developed as a result of natural convection in the melt as well as the significant differentiation of the measured velocity depending on the area and given the specific scale. In every case, the measured velocity acquires low values-typical of the buoyancy flow-which vary from 10 −5 m/s in the first 100 s up to 10 −3 m/s at 900 s. It is also important to stress the observed symmetry as to the vertical axis y in Figure 13. Farhadi et al. [30] studied the melting of organic PCM with a two-dimensional model of the same geometry as the one analyzed here. The melting formulation of this work is in agreement with the results of Farhadi.
The following figures represent the temporal evolution of the velocity field as well as magnified areas from the field at the point 900 s. Figure 12 illustrates the mode in which velocities are developed as a result of natural convection in the melt as well as the significant differentiation of the measured velocity depending on the area and given the specific scale. In every case, the measured velocity acquires low values-typical of the buoyancy flow-which vary from 10 −5 m/s in the first 100 s up to 10 −3 m/s at 900 s. It is also important to stress the observed  Within the areas at which melt circulation occurs, there are areas with near zero velocity at specific points. This can be observed at Figure 14, which represents a magnified area of the vector velocity field at point t = 900 s; it illustrates the vortex formed due to the natural circulation of the melt, which justifies the zero velocity at the center. As a result of this circulation, in which the hot Within the areas at which melt circulation occurs, there are areas with near zero velocity at specific points. This can be observed at Figure 14, which represents a magnified area of the vector velocity field at point t = 900 s; it illustrates the vortex formed due to the natural circulation of the melt, which justifies the zero velocity at the center. As a result of this circulation, in which the hot and light part of the melt rises until it reaches the solid edge or loses its heat load and cools down, several similar vortices and equal zero velocity areas are developed. Finally, the highest velocities per meter are observed at the confluence of individual vortices. Within the areas at which melt circulation occurs, there are areas with near zero velocity at specific points. This can be observed at Figure 14, which represents a magnified area of the vector velocity field at point t = 900 s; it illustrates the vortex formed due to the natural circulation of the melt, which justifies the zero velocity at the center. As a result of this circulation, in which the hot and light part of the melt rises until it reaches the solid edge or loses its heat load and cools down, several similar vortices and equal zero velocity areas are developed. Finally, the highest velocities per meter are observed at the confluence of individual vortices.
Finally, Figure 15 represents the temporal variation in the energy stored in the mass of the phase change material. This heat results from the completion of the enthalpy, the melt percentage fluctuations at the incremental surface of the material studied, and the ensuing conversion to the total surface of the PCM.  Finally, Figure 15 represents the temporal variation in the energy stored in the mass of the phase change material. This heat results from the completion of the enthalpy, the melt percentage fluctuations at the incremental surface of the material studied, and the ensuing conversion to the total surface of the PCM.

Nanoparticle-Enhanced PCM
The following figures represent the effects of adding nanoparticles on the melt percentage ( Figure 16) as well as the heat stored in the enhanced PCM ( Figure 17). The percentage increase in the melt appears to be unstable at the beginning of the process because the melt quantity in the first seconds is infinitely small, which causes abrupt fluctuation. With the increase in the melt quantity comes the stabilization of the percentage or at least a smoother variation. When melting begins and heat transfer by conduction prevails, the enhanced materials with the highest thermal conductivity, i.e., PCM-Cu (3%) and PCM-Al 2 O 3 (3%), have better performance in terms of both heat storage and material melting. After 300 s of operation and the production of an efficient melt quantity, natural circulation intensifies and natural conduction becomes the primary heat transfer mechanism. The outcome is that copper nanoparticles are mostly advantageous in enhancing the thermal behavior of the PCM, whereas PCM-Al 2 O 3 (3%), from the establishment of natural convection and onwards, falls behind in terms of melt growth even when compared with PCM-Al 2 O 3 (1%). The fact that PCM-Al 2 O 3 is less effective than PCMs-Cu may be attributed mainly to its lower thermal conductivity. One major parameter influencing the speed of melting and heat storage is the kinematic viscosity of the enhanced material, which increases with the addition of nanoparticles. The increase in kinematic viscosity inhibits the growth of melt circulation, thus suspending natural convection and, consequently, heat transfer from the hot wall to the enhanced PCM. Therefore, heat is transferred faster in the lower heat conductivity PCM-Cu (1%) than in PCM-Al 2 O 3 (3%), which is attributed to its low kinematic viscosity, especially after 300 s, when the circulation of the superheated melt has reached satisfactory levels. The importance of kinematic viscosity has even been stressed in the literature but is also prominent in the slight prevalence of PCM-Al 2 O 3 (1%) over PCM-Al 2 O 3 (3%).

Conclusions
Considering the two-dimensional analysis accomplished by means of developed code and available software, and which focuses on both pure organic PCM and the enhanced material with metal and metal oxide nanoparticles, the main conclusions are summarized as follows: • The complexity of the problem, especially in the "mushy zone" area, necessitates particularly long simulation periods and increased computing power for the solution. For example, the evolution of the phenomenon up to 7200 s requires approximately 160.6 h of simulation.

•
Buoyancy induces velocity fields with measured velocity of about 10 −3 m/s. This natural circulation of the melt intensifies the heat transfer and the phase change as it transfers heat from the melt to the solid by a combination of natural convection and conduction.

•
The velocity field takes the form of vortex with areas of zero velocity (vortex center) and areas of maximum velocity (confluence with neighboring vortex).

•
The addition of Cu and Al 2 O 3 nanoparticles to the compositions studied enhances both the melt percentage and the heat stored in the material. However, this effect declines over time for the case of the enhanced PCM bearing Al 2 O 3 nanoparticles. On the other hand, for the case of copper nanoparticles, this enhancement exhibits intense and irregular fluctuations for the first seconds of the process and then maintains stable values of about 6.5% for the melt percentage and 5.5% for the heat stored.

•
The quantity of the nanoparticles added cannot be increased infinitely, and the reason for this is the kinematic viscosity of the enhanced material, which increases with the addition of nanoparticles. The increase in kinematic viscosity inhibits the growth of melt circulation, thus suspending natural convection and, consequently, heat transfer from the hot wall to the enhanced PCM. Acknowledgments: Michael Nitsas would like to thank the State Scholarships Foundation (IKY) for the financial support (MIS-5033021). Both authors acknowledge with appreciation the Laboratory of Applied Thermodynamics, Thermal Engineering Section, School of Mechanical Engineering, National Technical University of Athens, for the support of the work on which this paper is based.

Conflicts of Interest:
The authors declare no conflict of interest.