Fixed Grid Numerical Models for Solidification and Melting of Phase Change Materials (PCMs)

Phase change materials (PCMs) are classified according to their phase change process, temperature, and composition. The utilization of PCMs lies mainly in the field of solar energy and building applications as well as in industrial processes. The main advantage of such materials is the use of latent heat, which allows the storage of a large amount of thermal energy with small temperature variation, improving the energy efficiency of the system. The study of PCMs using computational fluid dynamics (CFD) is widespread and has been documented in several papers, following the tendency that CFD nowadays tends to become increasingly widespread. Numerical studies of solidification and melting processes use a combination of formulations to describe the physical phenomena related to such processes, these being mainly the latent heat and the velocity transition between the liquid and the solid phases. The methods used to describe the latent heat are divided into three main groups: source term methods (E-STM), enthalpy methods (E-EM), and temperature-transforming models (E-TTM). The description of the velocity transition is, in turn, divided into three main groups: switch-off methods (SOM), source term methods (STM), and variable viscosity methods (VVM). Since a full numerical model uses a combination of at least one of the methods for each phenomenon, several combinations are possible. The main objective of the present paper was to review the numerical approaches used to describe solidification and melting processes in fixed grid models. In the first part of the present review, we focus on the PCM classification and applications, as well as analyze the main features of solidification and melting processes in different container shapes and boundary conditions. Regarding numerical models adopted in phase-change processes, the review is focused on the fixed grid methods used to describe both latent heat and velocity transition between the phases. Additionally, we discuss the most common simplifications and boundary conditions used when studying solidification and melting processes, as well as the impact of such simplifications on computational cost. Afterwards, we compare the combinations of formulations used in numerical studies of solidification and melting processes, concluding that “enthalpy–porosity” is the most widespread numerical model used in PCM studies. Moreover, several combinations of formulations are barely explored. Regarding the simulation performance, we also show a new basic method that can be employed to evaluate the computing performance in transient numerical simulations.


Introduction
Thermal energy storage systems can employ sensible or latent heat. Sensible heat storage technology is more widespread than latent heat storage; however, the use of latent heat storage presents the main advantage of a higher energy storage density [1]. Latent heat thermal energy storage with phase change materials (PCMs) has been extensively studied due to its potential in peak absorption both in energy demand and in energy supply [2]. The study of PCMs was popular during the energy crisis in the 1970s and 1980s, with the primary aim to use PCMs in solar thermal energy storage [3]. As a transient process, the time elapsed by the complete phase change of the PCM (charge or discharge) is essential to defining the viability of the latent heat thermal energy storage system (LHTESS). This period covered by charge or discharge can vary from low peak attenuation to daily base or long-term yearly base [4]. Therefore, the use of PCMs includes a wide range of applications, such as building passive climate control [5][6][7], concentrating solar power stations [8][9][10], electronic device thermal protection [11], textile products [12], and waste heat recovery [13,14].
Several authors have reviewed PCMs with regard to their applications, materials, numerical models, and container shape. Thus, some significant parameters have been established [15] for the application and study of PCMs, such as desirable characteristics of materials, requirements of the LHTESS, initial classification of PCMs and latent heat measurement techniques. Additional contributions have been performed by Zalba et al. [16], Farid et al. [1], Sharma et al. [2], and Agyenim et al. [3]. More specifically, works [17,18] are focused on PCMs for high-temperature applications, while Ge et al. [19] have reviewed the group of metallic PCMs. Salunkhe and Shembekar [20] have outlined the heat transfer characteristics for different size scales of encapsulation, as well as the effects of the PCM container material in the phase change process. Dhaidan and Khodadadi [21] focused their study on the melting process of PCMs and the effects of container shape on heat transfer and the melting process. Fan and Khodadadi [22] and Khan et al. [23] have analyzed the techniques adopted for thermal conductivity enhancement of thermal energy storage using PCMs. Verma et al. [24] have recounted thermal energy storage systems with PCMs: the models were separated according to the first and second laws of thermodynamics. Regarding the numerical approach, Zalba et al. [16] and Sharma et al. [2] have outlined some characteristics of the formulations used to describe phase-change processes. Dutil et al. [25] and Al-Abidi et al. [26] have reviewed several numerical studies of PCMs and have presented the characteristics of computational models according to the geometry of the PCM container. Liu et al. [27] have analyzed the heat transfer mechanisms of phase change as well as analytical and numerical models for phase change. They reviewed PCMs consistently with the division into two categories: considering only heat conduction or considering conduction and convection simultaneously. Ziskind [28] has pointed out physical phenomena associated with PCMs and has reviewed numerical methods used for modeling such phenomena. Some authors have presented specific reviews for numerical models, such as Voller et al. [29] and Samarskii et al. [30], in which phase change formulations and numerical procedures were pointed out. Hu and Argyropoulos [31] focused their study mainly on latent heat formulation. Ma and Zhang [32] have reviewed numerical models for the velocity transition between solid and liquid and performed tests and comparisons between different models.
Despite the large number of review papers regarding PCMs, very few are specifically about numerical modeling features. Among these reviews, the most common approach is to focus on latent heat modeling or velocity transition modeling. In the present paper, we review the formulations used to describe both physical phenomena and relate to the features that occur in solidification and melting processes. Thus, the present review has been divided into four main parts: (1) review of PCM classification and application; (2) analysis of the main features of phase change processes regarding cavity format and thermal conditions; (3) review of fixed grid numerical models used for phase-change processes; (4) review of simplifications and other conditions commonly used in numerical studies of PCMs, performance evolution in computers used for numerical simulations, and compilation of the combinations of phase change numerical models most common in the literature.

PCMs Classification
PCMs are classified according to composition, phase change process, and phase change temperature. More precisely, the latter turns out to be the main parameter for the choice of a PCM, which needs to be suitable for the temperature range that occurs in the considered application. In general, the classification according to the phase change process is divided into four subcategories: solid-solid, solid-liquid, liquid-gas, and solid-gas [19]. As for the classification concerning phase change temperature, three groups are universally accepted among the authors (low, medium, and high temperature), even if there is no clear definition of the temperature ranges for each group, as shown in Table 1. Table 1. Phase change materials classification, according to phase change temperature.

PCMs: Basic Characterization
Latent heat is present in every form of phase change: solid-solid, solid-liquid (solidification and melting), liquid-gas (boiling and condensation) or solid-gas (sublimation). In comparison with sensible heat, a large amount of energy is stored in phase changes, with small temperature variation. Phase change processes that involve a gas phase, such as boiling, condensation, and sublimation processes, are characterized by a considerable variation in the material volume [2,19]. Thus, these processes need pressure vessels, which make the thermal energy storage system more complex and less feasible [15]. The transition between two solid phases or between solid and liquid implies less volume variation [1] so that the increase in pressure is also moderate; therefore, solid-solid and solid-liquid phase changes of PCMs make the energy storage system less complex. Furthermore, PCMs with solid-solid transitions present as their main advantage the possibility of direct contact of the heat transfer fluid (HTF) with the PCM, as the shape of the solid remains during the phase change [33,34]. However, it is worth mentioning that fewer materials present a solid-solid transition, having commonly lower latent heat than that corresponding to the solid-liquid transition [13,18].

PCMs Classification
PCMs are classified according to composition, phase change process, and phase change temperature. More precisely, the latter turns out to be the main parameter for the choice of a PCM, which needs to be suitable for the temperature range that occurs in the considered application. In general, the classification according to the phase change process is divided into four subcategories: solid-solid, solid-liquid, liquid-gas, and solid-gas [19]. As for the classification concerning phase change temperature, three groups are universally accepted among the authors (low, medium, and high temperature), even if there is no clear definition of the temperature ranges for each group, as shown in Table 1. Table 1. Phase change materials classification, according to phase change temperature.

Reference
Low Temperature Medium Temperature High Temperature Abhat [15] <120 °C --Zalba et al. [16] <100 °C --Agyenim et al. [3] <60 °C Range 80-120°C >150 °C Gil et al. [8] -->120 °C Nomura et al. [18] -->100 °C Ge et al. [19] <30 °C Range 40-200°C >200 °C Therefore, the classification of PCMs according to composition is still the most common in the literature due to its intrinsic simplicity, since groups of materials with a similar composition also have similar characteristics. Several authors have presented such a classification, based on three main groups: organic, inorganic, and mixtures. Each main group has subgroups, as highlighted in Figure 1, according to the classifications shown in works [2,15,16,18,19].  Therefore, the classification of PCMs according to composition is still the most common in the literature due to its intrinsic simplicity, since groups of materials with a similar composition also have Appl. Sci. 2019, 9, 4334 4 of 22 similar characteristics. Several authors have presented such a classification, based on three main groups: organic, inorganic, and mixtures. Each main group has subgroups, as highlighted in Figure 1, according to the classifications shown in works [2,15,16,18,19].

PCM Applications
The use of PCMs covers several applications, with a higher amount of studies focused on solar energy and Heating, Ventilation and Air Conditioning (HVAC) Systems [16]. In addition, many other areas can benefit from using PCM systems, including the construction field, space crafts and satellites, waste heat recovery, photovoltaic-thermal solar collectors, electronic equipment, personal comfort, and medical systems, among others. Thermal energy storage always plays an essential role in concentrating solar power plants, as it increases the availability of the plant to periods without sunlight. In such systems, part of the solar energy collected during the day is redirected to a reservoir with PCM that keeps the electricity generated during the night [8]. An improvement in the performance of the thermal energy storage and power plant is the use of a cascaded storage system [9,10]. Cascaded PCM systems are composed of several PCMs with different phase change temperatures and are configured with the phase change temperature in ascending order. Thus, the direction of the HTF flow is reversed between charge and discharge. In photovoltaic solar collectors, the usage of PCMs results in a reduction in the peak temperature. Thus, the losses due to overheating are reduced in the collector [14]. Another important application of PCMs is in photovoltaic-thermal solar collectors, where PCMs have been observed to improve the overall efficiency of the collector [12,35,36]. The same occurs in the utilization of PCMs in solar distillation processes in which the operational period after sunset is extended. As a result, both efficiency and productivity are increased [37,38]. Another essential building application of PCMs is the thermal storage in solar heating systems for domestic water, as pointed out in works [15]. Commonly, such systems use water as sensible heat storage due to the low cost of installation. The use of PCMs as latent heat storage in solar heating systems entails the main advantages of the reduction in reservoir size and improvement of the overall system efficiency [39]. However, the cost of the PCM may represent an obstacle to large-scale employment in domestic hot water systems [16]. In air HVAC systems, for example, the utilization of PCMs with phase change temperature lower than ambient temperature can bring about several advantages in order to reduce the peak energy load and slightly increase the efficiency of the equipment [40]. In the construction field, PCMs are used for passive storage daily, in the form of microcapsules contained in thermal insulation boards. The use of these boards provides temperature peak attenuation inside the building, lowering the maximum temperature during the day and increasing the minimum temperature during the night. Therefore, the thermal load is reduced, and consequently, the load of heating or air conditioning is also reduced [2]. Space crafts and satellites orbiting the Earth operate in a cyclical thermal condition that varies between extremes. When the spacecraft is exposed to the sun, the temperature reaches its maximum. Oppositely, when the spacecraft is in the shadow of the Earth, the temperature reaches its minimum. Thus, the use of PCM attenuates the high-and low-temperature peaks and represents great potential to reduce energy consumption [41]. Even waste heat recovery in industries using PCMs represents excellent potential for energy saving [18]. A representative example is the steelmaking industry, in which vast amounts of heat are rejected to cool the alloys produced [42]. According to Yagi and Akiyama [43], after having been stored in the PCM, the waste heat can be directed to a supply center for public use. In vehicles with internal combustion engines and hydraulic machines, waste heat can be directly stored in the PCM and used to preheat the engine, as a cold start represents higher fuel consumption and more significant wear of the moving parts [16]. The use of PCMs is attractive in the automotive industry for battery thermal management, mainly in hybrid and electric vehicles. High and low temperature can influence the battery performance and in more critical cases, overheating of the battery can damage or even cause failure in such a component, causing safety risks [44]. The prevention of battery failures due to overheating is also crucial in electronic equipment. In this context, work [11] has outlined the application of PCMs in such equipment in order to increase the heat dissipation capability. Textile applications of encapsulated PCMs include clothing for thermal comfort in hot and cold environments, space suits, sportswear, medical applications, shoes, and accessories [45]. Heat storage with PCMs also represents an attractive option in solar dryers for agricultural and food products [46]. Additional examples of PCMs applications are in medicine, including hot-cold therapies, food conservation, and transportation [16].

Phase Change Processes
Phase change processes can occur under a single temperature (isothermal phase change) or within a temperature interval, also known as continuous properties phase change [47]. In the latter case, the properties vary continuously within the interval, known as the mushy zone. Figure 2 illustrates the enthalpy variation with temperature for isothermal phase change (dashed line) and continuous properties phase change (solid line).

Phase Change Processes
Phase change processes can occur under a single temperature (isothermal phase change) or within a temperature interval, also known as continuous properties phase change [47]. In the latter case, the properties vary continuously within the interval, known as the mushy zone. Figure 2 illustrates the enthalpy variation with temperature for isothermal phase change (dashed line) and continuous properties phase change (solid line). The phase change temperature is Tm for the isothermal phase change, while the mushy zone in the continuous properties phase change is limited by solidus temperature (Ts) in the lower bound and by liquid temperature (Tl) in the upper bound. Thus, solid and liquid phases coexist in the mushy zone, whose thickness is commonly proportional to the temperature interval, being therefore greater than zero. The temperature interval in continuous properties phase change is commonly approximated to an equivalent Tm to simplify and include parameters in the analysis, as well as apply boundary conditions, perform calculations, and obtain results. The most common and simplest approximation used in works [48,49] is given by: Examples of materials with isothermal phase change are pure substances and eutectic mixtures, as highlighted in work [50]. Continuous properties phase change exists in metallic alloys, mixtures, and amorphous substances, as illustrated in works [47,51].

Solidification Process
In the solidification process, heat transfer is predominantly conductive [52]. However, the effects of convection can be observed in the case of rectangular cavities with heat transfer on one side [53]. In such domains, convection heat transfer affects the shape of the solid along the process. The main characteristic detected in the solidification process is the formation of a solid layer, whose thickness increases gradually from the heat transfer surface [49,54]. The solidification rate is reduced throughout the process, with its effect observed through the smaller variation of the liquid fraction in the final stages. The reduction in the solidification rate is shown in detail on the experimental results obtained by Viskanta and Gau [55] for the solidification of water inside a horizontal cylinder.
The solid fraction is depicted in a dark color and the liquid fraction in a light color, as the darker portions in the center of the images are the thermocouple and overflow tube. The time interval elapsed between images (b) and (c) is higher than between images (a) and (b), even if the corresponding amount of solidified material has an inverse proportion. Also, it is possible to observe in the images a symmetric solidification pattern from the colder surface towards the center. It is worth mentioning that the main effects, as a reduction in solidification rate and symmetrical pattern The phase change temperature is T m for the isothermal phase change, while the mushy zone in the continuous properties phase change is limited by solidus temperature (T s ) in the lower bound and by liquid temperature (T l ) in the upper bound. Thus, solid and liquid phases coexist in the mushy zone, whose thickness is commonly proportional to the temperature interval, being therefore greater than zero. The temperature interval in continuous properties phase change is commonly approximated to an equivalent T m to simplify and include parameters in the analysis, as well as apply boundary conditions, perform calculations, and obtain results. The most common and simplest approximation used in works [48,49] is given by: Examples of materials with isothermal phase change are pure substances and eutectic mixtures, as highlighted in work [50]. Continuous properties phase change exists in metallic alloys, mixtures, and amorphous substances, as illustrated in works [47,51].

Solidification Process
In the solidification process, heat transfer is predominantly conductive [52]. However, the effects of convection can be observed in the case of rectangular cavities with heat transfer on one side [53]. In such domains, convection heat transfer affects the shape of the solid along the process. The main characteristic detected in the solidification process is the formation of a solid layer, whose thickness increases gradually from the heat transfer surface [49,54]. The solidification rate is reduced throughout the process, with its effect observed through the smaller variation of the liquid fraction in the final stages. The reduction in the solidification rate is shown in detail on the experimental results obtained by Viskanta and Gau [55] for the solidification of water inside a horizontal cylinder.
The solid fraction is depicted in a dark color and the liquid fraction in a light color, as the darker portions in the center of the images are the thermocouple and overflow tube. The time interval elapsed between images (b) and (c) is higher than between images (a) and (b), even if the corresponding amount of solidified material has an inverse proportion. Also, it is possible to observe in the images a symmetric solidification pattern from the colder surface towards the center. It is worth mentioning that the main effects, as a reduction in solidification rate and symmetrical pattern from the cold surface, are linked to the different configurations. Thus, regarding geometry, these effects are observed in vertical cylinders [56], horizontal concentric cylinders with different fin configurations [57,58] and spheres [52,59].

Melting Process
In the melting process, heat conduction is dominant only in the initial stages [58]. In any case, convective heat transfer influences most of the melting process [60], as can be detected in the morphology of the solid phase. Regarding rectangular cavities with heat transfer in one side, the liquid phase close to the hot wall flows upwards due to buoyancy forces induced by density reduction in the function of temperature increase. The convective motions heated by the hot wall cause a higher melting rate in the upper region, as the solid assumes a tapered shape that becomes more evident over time, as shown in Figure 3 about experimental results outlined in work [61].

Melting Process
In the melting process, heat conduction is dominant only in the initial stages [58]. In any case, convective heat transfer influences most of the melting process [60], as can be detected in the morphology of the solid phase. Regarding rectangular cavities with heat transfer in one side, the liquid phase close to the hot wall flows upwards due to buoyancy forces induced by density reduction in the function of temperature increase. The convective motions heated by the hot wall cause a higher melting rate in the upper region, as the solid assumes a tapered shape that becomes more evident over time, as shown in Figure 3 about experimental results outlined in work [61]. Regarding the shape of the solid phase, vertical cylinders with an adiabatic base and heat transfer on the cylindrical wall present results similar to rectangular cavities with heat transfer on one side. Vertical cylinders have symmetry around the vertical axis, which induces the solid phase to assume a conical shape, as shown in Figure 4, concerning experimental results obtained by Katsman et al. [62]. Besides the shape of a solid PCM, the adiabatic bottom surface is an additional standard feature regarding PCMs melting in rectangular cavities and vertical cylinders with a one side heated surface. The absence of heat transfer through the bottom surface maintains the position of the remaining solid phase throughout the process. However, most PCM containers present heat transfer through the bottom surface, which is more evident in horizontal cylinders and spheres but can occur as well in rectangular cavities and vertical cylinders. As different densities characterize solid and liquid PCMs, the portions with lower density tend to flow upwards. The density of the PCM is higher for the solid than for the liquid. Thus, in melting processes with heat transfer through the lower surface, the liquid portion flows upwards and the solid is located close to the heated lower surface. For materials with higher density in the liquid phase, such as water, the opposite occurs. Regarding the shape of the solid phase, vertical cylinders with an adiabatic base and heat transfer on the cylindrical wall present results similar to rectangular cavities with heat transfer on one side. Vertical cylinders have symmetry around the vertical axis, which induces the solid phase to assume a conical shape, as shown in Figure 4, concerning experimental results obtained by Katsman et al. [62].

Melting Process
In the melting process, heat conduction is dominant only in the initial stages [58]. In any case, convective heat transfer influences most of the melting process [60], as can be detected in the morphology of the solid phase. Regarding rectangular cavities with heat transfer in one side, the liquid phase close to the hot wall flows upwards due to buoyancy forces induced by density reduction in the function of temperature increase. The convective motions heated by the hot wall cause a higher melting rate in the upper region, as the solid assumes a tapered shape that becomes more evident over time, as shown in Figure 3 about experimental results outlined in work [61]. Regarding the shape of the solid phase, vertical cylinders with an adiabatic base and heat transfer on the cylindrical wall present results similar to rectangular cavities with heat transfer on one side. Vertical cylinders have symmetry around the vertical axis, which induces the solid phase to assume a conical shape, as shown in Figure 4, concerning experimental results obtained by Katsman et al. [62]. Besides the shape of a solid PCM, the adiabatic bottom surface is an additional standard feature regarding PCMs melting in rectangular cavities and vertical cylinders with a one side heated surface. The absence of heat transfer through the bottom surface maintains the position of the remaining solid phase throughout the process. However, most PCM containers present heat transfer through the bottom surface, which is more evident in horizontal cylinders and spheres but can occur as well in rectangular cavities and vertical cylinders. As different densities characterize solid and liquid PCMs, the portions with lower density tend to flow upwards. The density of the PCM is higher for the solid than for the liquid. Thus, in melting processes with heat transfer through the lower surface, the liquid portion flows upwards and the solid is located close to the heated lower surface. For materials with higher density in the liquid phase, such as water, the opposite occurs. Melting processes in which the solid material moves towards the heated surface or vice-versa are Besides the shape of a solid PCM, the adiabatic bottom surface is an additional standard feature regarding PCMs melting in rectangular cavities and vertical cylinders with a one side heated surface. The absence of heat transfer through the bottom surface maintains the position of the remaining solid phase throughout the process. However, most PCM containers present heat transfer through the bottom surface, which is more evident in horizontal cylinders and spheres but can occur as well in rectangular cavities and vertical cylinders. As different densities characterize solid and liquid PCMs, the portions with lower density tend to flow upwards. The density of the PCM is higher for the solid than for the liquid. Thus, in melting processes with heat transfer through the lower surface, the liquid portion flows upwards and the solid is located close to the heated lower surface. For materials with higher density in the liquid phase, such as water, the opposite occurs. Melting processes in which the solid material moves towards the heated surface or vice-versa are known as "close contact melting" [63] or "unconstrained melting" [58]. The main characteristic of this process is the formation of a thin liquid layer between the solid and the heated surface. As more material melts, the liquid portion formed in this layer flows continuously towards the side edge of the solid portion, out of the region between the solid and the heated surface [64]. An example of close contact melting in a rectangular cavity is shown in Figure 5, adapted from work [65].
Appl. Sci. 2019, 9, x 7 of 22 known as "close contact melting" [63] or "unconstrained melting" [58]. The main characteristic of this process is the formation of a thin liquid layer between the solid and the heated surface. As more material melts, the liquid portion formed in this layer flows continuously towards the side edge of the solid portion, out of the region between the solid and the heated surface [64]. An example of close contact melting in a rectangular cavity is shown in Figure 5, adapted from work [65]. The shape of the solid (bright color) remains similar to the cavity shape. In comparison with the case of the rectangular cavity highlighted in Figure 3, the influence of the thermal configuration in the shape of the solid portion proves to be evident, even with the different aspect ratios of the cavities. The differences in the phase morphology when comparing close contact melting and fixed solid melting are also observable in horizontal cylinders as well as in spheres. However, other results, such as melting rate, heat flux, or total melting time, are different due to geometric characteristics. When the solid is fixed, the melting process starts symmetrically. However, over time, a wavy profile is formed in the lower region of the solid, due to convective motions of the liquid at the bottom of the cavity [58]. On the other hand, close contact melting presents shapes that are more symmetrical in the solid. Figure 6 shows the shapes of the solid PCM in the melting process in the horizontal cylinder analyzed by Sparrow and Geiger [66]. Regarding the case with the fixed solid, the wavy profile in the lower portion of the solid is visible both longitudinally ( Figure 6a) and transversally (Figure 6b). The case of close contact melting (Figure 6c) does not show the wavy profile, but the solid keeps the shape more symmetrical and similar to the cavity shape. The melting process in spherical containers presents solid shapes like the transversal view of the horizontal cylinder, as depicted in Figure 7, deducted from work [58]. The shape of the solid (bright color) remains similar to the cavity shape. In comparison with the case of the rectangular cavity highlighted in Figure 3, the influence of the thermal configuration in the shape of the solid portion proves to be evident, even with the different aspect ratios of the cavities. The differences in the phase morphology when comparing close contact melting and fixed solid melting are also observable in horizontal cylinders as well as in spheres. However, other results, such as melting rate, heat flux, or total melting time, are different due to geometric characteristics. When the solid is fixed, the melting process starts symmetrically. However, over time, a wavy profile is formed in the lower region of the solid, due to convective motions of the liquid at the bottom of the cavity [58]. On the other hand, close contact melting presents shapes that are more symmetrical in the solid. Figure 6 shows the shapes of the solid PCM in the melting process in the horizontal cylinder analyzed by Sparrow and Geiger [66].
Appl. Sci. 2019, 9, x 7 of 22 known as "close contact melting" [63] or "unconstrained melting" [58]. The main characteristic of this process is the formation of a thin liquid layer between the solid and the heated surface. As more material melts, the liquid portion formed in this layer flows continuously towards the side edge of the solid portion, out of the region between the solid and the heated surface [64]. An example of close contact melting in a rectangular cavity is shown in Figure 5, adapted from work [65]. The shape of the solid (bright color) remains similar to the cavity shape. In comparison with the case of the rectangular cavity highlighted in Figure 3, the influence of the thermal configuration in the shape of the solid portion proves to be evident, even with the different aspect ratios of the cavities. The differences in the phase morphology when comparing close contact melting and fixed solid melting are also observable in horizontal cylinders as well as in spheres. However, other results, such as melting rate, heat flux, or total melting time, are different due to geometric characteristics. When the solid is fixed, the melting process starts symmetrically. However, over time, a wavy profile is formed in the lower region of the solid, due to convective motions of the liquid at the bottom of the cavity [58]. On the other hand, close contact melting presents shapes that are more symmetrical in the solid. Figure 6 shows the shapes of the solid PCM in the melting process in the horizontal cylinder analyzed by Sparrow and Geiger [66]. Regarding the case with the fixed solid, the wavy profile in the lower portion of the solid is visible both longitudinally ( Figure 6a) and transversally (Figure 6b). The case of close contact melting (Figure 6c) does not show the wavy profile, but the solid keeps the shape more symmetrical and similar to the cavity shape. The melting process in spherical containers presents solid shapes like the transversal view of the horizontal cylinder, as depicted in Figure 7, deducted from work [58]. Regarding the case with the fixed solid, the wavy profile in the lower portion of the solid is visible both longitudinally ( Figure 6a) and transversally (Figure 6b). The case of close contact melting (Figure 6c) does not show the wavy profile, but the solid keeps the shape more symmetrical and similar to the cavity shape. The melting process in spherical containers presents solid shapes like the transversal view of the horizontal cylinder, as depicted in Figure 7, deducted from work [58]. Appl. Sci. 2019, 9, x 8 of 22 Based on pure observation, both close contact melting ( Figure 7a) and fixed solid melting (Figure 7b) show similar solid shapes of the transversal view detected in the case of a horizontal cylinder (Figure 6c,b, respectively). Another common characteristic of close contact melting is the higher melting rate for the same geometry and same thermal condition when compared with unconstrained melting. During the close contact melting process, the solid keeps moving towards the heated surface while the liquid layer thickness is stable. This condition provides a nearly constant liquid flow upwards and consequently a nearly constant melting rate, as outlined in work [63]. The liquid layer thickness increases due to the reduction of the body forces from the solid portion, as the solid mass reduces over time [65] at the end of the process. Visualization of the liquid layer below the solid is very difficult to be observed, due to its very small thickness [63,64]. Most of the PCM melted in the entire close contact melting process comes from this liquid layer between the solid and the lower surface. However, some authors have detected a contribution of about 10-15% from natural convection in the upper region for both spheres [64] and horizontal cylinders [67]. Figure 8 shows, schematically, the liquid flow during the close contact melting inside a sphere, where it is possible to observe the liquid streams of natural convection in the upper region and the liquid layer between the solid and the lower heated surface. The liquid amount in the layer below the solid increases from the symmetry axis as the entire surface is heated. This increase causes the variation in the velocity profile of the liquid layer towards the upper region, which is accompanied by an augmentation of the layer thickness [68]. For a given position, the velocity profile does not vary significantly [69] until the upper region of the solid reaches that position. Several studies are based on equations and correlations to describe parameters in the close contact melting such as liquid layer thickness [68,69], melting rate in the liquid layer  (Figure 6c,b, respectively). Another common characteristic of close contact melting is the higher melting rate for the same geometry and same thermal condition when compared with unconstrained melting. During the close contact melting process, the solid keeps moving towards the heated surface while the liquid layer thickness is stable. This condition provides a nearly constant liquid flow upwards and consequently a nearly constant melting rate, as outlined in work [63]. The liquid layer thickness increases due to the reduction of the body forces from the solid portion, as the solid mass reduces over time [65] at the end of the process. Visualization of the liquid layer below the solid is very difficult to be observed, due to its very small thickness [63,64]. Most of the PCM melted in the entire close contact melting process comes from this liquid layer between the solid and the lower surface. However, some authors have detected a contribution of about 10-15% from natural convection in the upper region for both spheres [64] and horizontal cylinders [67]. Figure 8 shows, schematically, the liquid flow during the close contact melting inside a sphere, where it is possible to observe the liquid streams of natural convection in the upper region and the liquid layer between the solid and the lower heated surface. Based on pure observation, both close contact melting ( Figure 7a) and fixed solid melting (Figure 7b) show similar solid shapes of the transversal view detected in the case of a horizontal cylinder (Figure 6c,b, respectively). Another common characteristic of close contact melting is the higher melting rate for the same geometry and same thermal condition when compared with unconstrained melting. During the close contact melting process, the solid keeps moving towards the heated surface while the liquid layer thickness is stable. This condition provides a nearly constant liquid flow upwards and consequently a nearly constant melting rate, as outlined in work [63]. The liquid layer thickness increases due to the reduction of the body forces from the solid portion, as the solid mass reduces over time [65] at the end of the process. Visualization of the liquid layer below the solid is very difficult to be observed, due to its very small thickness [63,64]. Most of the PCM melted in the entire close contact melting process comes from this liquid layer between the solid and the lower surface. However, some authors have detected a contribution of about 10-15% from natural convection in the upper region for both spheres [64] and horizontal cylinders [67]. Figure 8 shows, schematically, the liquid flow during the close contact melting inside a sphere, where it is possible to observe the liquid streams of natural convection in the upper region and the liquid layer between the solid and the lower heated surface. The liquid amount in the layer below the solid increases from the symmetry axis as the entire surface is heated. This increase causes the variation in the velocity profile of the liquid layer towards the upper region, which is accompanied by an augmentation of the layer thickness [68]. For a given position, the velocity profile does not vary significantly [69] until the upper region of the solid reaches that position. Several studies are based on equations and correlations to describe parameters in the close contact melting such as liquid layer thickness [68,69], melting rate in the liquid layer The liquid amount in the layer below the solid increases from the symmetry axis as the entire surface is heated. This increase causes the variation in the velocity profile of the liquid layer towards the upper region, which is accompanied by an augmentation of the layer thickness [68]. For a given position, the velocity profile does not vary significantly [69] until the upper region of the solid reaches that position. Several studies are based on equations and correlations to describe parameters in the close contact melting such as liquid layer thickness [68,69], melting rate in the liquid layer [70,71], and the solid portion sinking [63,72]. However, until now, no direct experimental measurements of the liquid layer thickness or velocity profile have been found in the literature, due to its small thickness, which turns out to be very challenging to observe directly. Regarding numerical investigations on this topic, some authors have studied the liquid layer thickness [73,74] and the downward velocity of the solid [75,76]. Despite the actual numerical results, more effort is necessary to provide more significant results for different geometries and concerning the heat transfer and the flow inside the liquid layer.

Numerical Models for Phase Change Processes
Isothermal phase change problems are also known as moving boundary problems or Stefan problems, due to the study of Stefan [77] of the ice cap formation in the North Pole region. The model was elaborated considering the interface between solid and liquid as a moving boundary with temperature conditions (T = T m ). The solution consisted of solving the interface position over time. Such a model presented a solution about some specific cases, i.e., considering several simplifications in the properties of the material, boundary, and initial conditions, as pointed out in work [31]. It is worth mentioning that the temperature condition in the solid-liquid interface corresponds to an isothermal phase change process, with a sharp interface. However, when referring to phase changes with continuous properties, both solid and liquid coexist in the mushy zone, and the interface is not a well-defined boundary [78]. Therefore, numerical models are most suitable for continuous properties phase change as the mushy zone represents a steep gradient in the domain. On the other hand, for isothermal phase change, the discontinuity represented by phase change is the main difficulty in managing with numerical models [47]. Numerical phase change models found in the literature are commonly divided into two main groups: fixed grid or continuum models, and variable grid or two-phase models [31]. Variable grid models are based on the implementation of separate computational domains for solid and liquid, i.e., separate equation sets, which are coupled using transference terms between the phases [79]. Thus, such models are more suitable to isothermal phase change processes, as they make possible the definition of an interface whose thickness tends to zero [78]. Besides the higher amount of equations, these procedures need interface tracking methods and grid adjustment methods according to the interface position, consequently presenting more difficulties to solving the problem. The grid position adjustment following the interface can be made through the movement of the grid points (dynamic grid) or through the variation of the time steps (interface fitting). Both adjustments are performed to fit the interface position to the grid points [31]. Variable grid models were implemented by Ismail et al. [54,80] with a dynamic grid and by Gupta and Kumar [81] with interface fitting.
Fixed grid models use a single computational domain with a single set of equations applied to the entire domain: continuity, energy conservation, and momentum equations. Therefore, the implementation of fixed grid models is simpler than variable grid models with the results always being accurate [32]. As fixed grid models have no interface tracking, the use of such models is useful to describe continuous properties phase change processes, since in the produced results a gradient between solid and liquid appears [78]. Therefore, fixed grid models are present in a more significant amount in the literature. The regular equations applied in fixed grid models are usually modified to contemplate the physical phenomena associated with phase change, namely the latent heat and the velocity transition between solid and liquid. The velocity transition occurs as liquid behavior follows the laws of fluid mechanics (with shear stress causing deformation rate), while the solid keeps static (with shear stress producing deformation rate). Fixed grid models are classified according to latent heat modeling and velocity transition modeling. Latent heat models are applied in the energy equation and can be based on source term (E-STM), enthalpy (E-EM) or temperature; in this case, they are also known in the literature as temperature-transforming models (E-TTM) [82]. Temperature-based models are divided into three categories: apparent heat capacity (AHC), effective heat capacity (EHC), and heat integration (HI), as reported in work [31]. Likewise, velocity transition models are divided into three main groups: switch-off method (SOM), source term method (STM), and variable viscosity method (VVM). The SOM and STM groups have in turn some smoothed versions, namely, the ramped switch-off method (RSOM), the ramped source term method (RSTM), and the Darcy STM [32]. Figure 9 resumes the fixed grid methods, based principally on the classification of Hu and Argyropoulos [31] and Ma and Zhang [32].
Appl. Sci. 2019, 9, x 10 of 22 smoothed versions, namely, the ramped switch-off method (RSOM), the ramped source term method (RSTM), and the Darcy STM [32]. Figure 9 resumes the fixed grid methods, based principally on the classification of Hu and Argyropoulos [31] and Ma and Zhang [32]. The following paragraphs are aimed at the description of the main characteristics of each of the above-mentioned fixed grid methods. Besides the main groups of the fixed grid and variable grid methods, some recent models are not exactly suited to either of them. Some examples are concerned with intermediate models in which the grid is fixed, but the solid is displayed as a rigid block, and the interface is continuously tracked to perform a force balance between the solid block and the surrounding liquid [74,83]. Moreover, some authors have arranged approaches based on particle physics, such as the Lattice Boltzmann method [84,85] or the Monte Carlo method [86,87].

Latent Heat Modeling
Latent heat is modeled through modifications in the energy equation. The amendment performed in the equation also defines the classification: inclusion of a source term (E-STM), energy equation based on temperature and specific heat (E-TTM), and energy equation based on enthalpy (E-EM). On a general basis, the equation of energy conservation assumes the form of Equation (2), i.e., where ρ is the density, Cp is the specific heat capacity, T the temperature, t the time, V the velocity vector, and k the thermal conductivity.

Source Term Method (E-STM)
The use of an additional source term in the energy equation allows modeling the latent heat in several forms. In such a method, the energy equation (Equation (2)) presents an additional source term that is equivalent to the latent heat. As an example, the source term (F) arranged by Voller and Swaminathan [88] is shown in Equation (3), i.e., where L is the latent heat and β is the liquid fraction, which varies from zero (solid phase) to 1 (liquid phase). The main advantage of source term methods is the possibility of using several equation formats to describe the enthalpy variation during the phase change processes, including cases with The following paragraphs are aimed at the description of the main characteristics of each of the above-mentioned fixed grid methods. Besides the main groups of the fixed grid and variable grid methods, some recent models are not exactly suited to either of them. Some examples are concerned with intermediate models in which the grid is fixed, but the solid is displayed as a rigid block, and the interface is continuously tracked to perform a force balance between the solid block and the surrounding liquid [74,83]. Moreover, some authors have arranged approaches based on particle physics, such as the Lattice Boltzmann method [84,85] or the Monte Carlo method [86,87].

Latent Heat Modeling
Latent heat is modeled through modifications in the energy equation. The amendment performed in the equation also defines the classification: inclusion of a source term (E-STM), energy equation based on temperature and specific heat (E-TTM), and energy equation based on enthalpy (E-EM). On a general basis, the equation of energy conservation assumes the form of Equation (2), i.e., where ρ is the density, C p is the specific heat capacity, T the temperature, t the time, → V the velocity vector, and k the thermal conductivity.

Source Term Method (E-STM)
The use of an additional source term in the energy equation allows modeling the latent heat in several forms. In such a method, the energy equation (Equation (2)) presents an additional source term that is equivalent to the latent heat. As an example, the source term (F) arranged by Voller and Swaminathan [88] is shown in Equation (3), i.e., where L is the latent heat and β is the liquid fraction, which varies from zero (solid phase) to 1 (liquid phase). The main advantage of source term methods is the possibility of using several equation formats to describe the enthalpy variation during the phase change processes, including cases with discontinuities. This approach can be useful in the study of metallic alloys when the influence of the composition and solubility of the components is considered in the phase change [88].

Temperature Transforming Methods (E-TTM)
In the temperature methods, such as the one presented by Morgan [89], the energy equation assumes the basic form shown in Equation (2). Latent heat is modeled as an increase in the specific heat capacity within the temperature interval of the phase change. Such methods are stable in the convergence for continuous properties phase change processes, even if the convergence is difficult in the case of the reduction concerning the phase change interval [32]. As E-TTM needs a temperature interval to model the latent heat, the accuracy of the results is limited when applied to isothermal phase change [27]. The apparent heat capacity method uses the basic formulation for the specific heat capacity shown in Equation (4) [31], i.e., where C ap is the apparent heat capacity and the subscripts "s" and "l" refer to the solid and liquid phases, respectively. The main issue existing in the apparent heat capacity method emerges when significant temperature changes occur (in a time step) in a control volume that is close to the phase change temperature. In such situations, latent heat is not accounted for if the temperature variation is greater than the phase change interval. Therefore, small-time steps are needed, with the implication of higher computational costs [31] to manage this condition. The effective heat capacity method [90] is an improvement on the apparent heat capacity method. The value of the effective heat capacity (C eff ) is obtained through the integration in the control volume, as shown in Equation (5), i.e. , where V is the control volume. It is worth mentioning that the adoption of the effective heat capacity method allows the use of higher time steps and smaller phase change intervals when compared with the apparent heat capacity method. However, the integration necessary for the calculation of effective heat capacity represents an increase in computational cost, as pointed out in work [31]. In the heat integration or post-iterative method, the solution is achieved in two main steps. The first step consists of the resolution of the energy equation without considering the latent heat, and the second step consists of the addition of equivalent latent heat. In the control volumes where the temperature reaches or exceeds the phase change interval, the temperature is set back to the beginning of the phase change interval and the energy is adjusted. The modification consists of the addition of an amount of latent heat equivalent to the energy transferred from the beginning of the phase change interval [90]. As the latent heat is added, the temperature in the remaining control volumes is corrected from the cells in the phase change interval. When the amount of energy added in the control volume reaches the total latent heat, the solution is performed, usually considering only sensible heat [31].

Enthalpy-Based Methods
In enthalpy-based methods [51,91], the energy conservation equation assumes the form shown in the following equation, i.e., ∂(ρh) ∂t where h is the total enthalpy, which is formed by sensible and latent parts. Enthalpy can be modeled in different ways; for example, Swaminathan and Voller [91] have proposed Equation (7), i.e., where T ref is the reference temperature. Enthalpy-based methods have, as their main advantage, the ability to simulate isothermal or continuous properties among phase change properties. However, such methods can be affected by temperature oscillations along the time at a certain point of the domain, as highlighted in work [32].

Velocity Transition Modeling
The objective of modeling the velocity transition is to describe the change between the flow of the liquid phase and the static behavior of the solid. The velocity transition is described through adjustments in the momentum equation. Such modifications are divided into three main groups: SOM, STM, and VVM. The momentum equation in its basic form, without modifications, is shown in Equation (8), i.e., where p is the pressure, µ is the dynamic viscosity, and g is the gravity acceleration.

Switch-Off Methods
Switch-off methods are based on a simple way to model the velocity transition between solid and liquid, as the velocity is set to zero in the solid phase [92]. The velocity modification function can be done directly in the numerical code in order to overwrite the velocity value to zero in the control volumes, which have a temperature below the phase change value [89]. In such methods, the momentum equation has no modifications, as the velocity transition is done as a function of temperature. However, SOM can produce discontinuities in the velocity field and cause inconsistencies in the results. Ramped switch-off methods have been developed to reduce such discontinuities through functions to smooth the transition. The functions usually use numerical coefficients with high value to suppress the solid motion [25]. The adoption of such functions results in an interface region with thickness different from zero, where the velocity transition occurs [32]. Given the condition of zero velocity implied to the solid, SOM and RSOM can provide accurate results for fixed solid melting. However, such methods are not adequately suitable for close contact melting, where the solid velocity is different from zero.

Source Term Methods
In the source term methods, a source term (S) is added to the conventional momentum equation (Equation (8)), assuming a high value when the material is solid [32]. Thus, to solve the momentum equation, the velocity is suppressed in the solid region to reach equilibrium with the high value of the source term. It is worth mentioning that discontinuities can be produced in the velocity field, similarly to SOM methods. Ramped source term methods have been developed in order to reduce problems with discontinuities, similarly to RSOM [25]. Thus, ramped methods are more suitable to be applied with E-TTM than conventional STM and SOM due to the smoother transition provided within the phase change interval [27]. Darcy STM proposed by Voller and Prakash in work [51] represents a variant of RSTM that is generally considered as a separated method, since it is largely used in studies of phase change processes. In such a method, the source term is based on the Darcy law for the flow in porous media, according to the equation where ε is a small value (ε = 0.001) used to avoid dividing by zero and C is the mushy zone constant, which is adjusted according to the case analyzed. The value of C must be large enough to allow significant flow in the mushy region and to suppress the velocity in the solid phase [51].

Variable Viscosity Method
The variable viscosity method, proposed by Gartling [93] uses the momentum equation in the basic form shown in Equation (8): it consists essentially in the increase of the viscosity referred to the solid phase. The higher viscosity of the solid induces the velocity reduction through the augmentation in the diffusive terms of the momentum equation. Fewer studies using VVM are found in the literature based on the comparison with SOM or STM groups. Although accurate results are delivered [75,92], VVM usually requires more calculation time as the solid viscosity value is increased [76], causing some possible discretization difficulties [92]. Different functions are found in the literature to represent viscosity variation using a liquid fraction, temperature, or solid fraction. In this context, Table 2 summarizes the viscosity functions and solid-liquid viscosity ratios (µ s /µ l ) used by several authors [94]: Table 2. Viscosity functions adopted in VVM studies.

Numerical Models Used in Studies Regarding PCMs
As solidification and melting processes have latent heat and velocity transitions between the phases, a proper numerical approach uses a combination of at least one of each group of methods described in the previous sections. In numerical models, the use of simplifying conditions is common to speed up convergence and allow the solution in a reasonable time. Such simplifying conditions can be assumed regarding geometry, material properties, physical phenomena associated with the problem, and boundary conditions. Among the geometrical simplifications, it is worth mentioning the use of a reduced domain (one-or two-dimensional) and the application of symmetry conditions: the former needs fewer equations to be solved and the latter reduces the size of the domain and consequently the number of cells. Therefore, the computational cost is dramatically reduced by adopting the simplifications mentioned above. Regarding material properties, simplifications include the assumption of continuous properties in each phase, independently of temperature, or even continuous properties for both phases, i.e., no property difference between solid and liquid. Another useful simplifying condition is the well-known Boussinesq approach, which considers the density constant for all equations, except for the body forces term in the momentum equation [51]. In such an assumption, convective motions are included in the model with low additional processing. Simplifications regarding physical phenomena contemplate mainly the conductive heat transfer, which represents a massive reduction in computation, as suppresses all the velocity terms, resuming the model to heat transfer equation with only conduction and phase change. However, such simplification limits the accuracy of the results, as primarily shown in the analysis of the phase morphology (see Section 3). Other physical simplifications may include the suppression of surface tension between solid and liquid or to consider the material isotropic. Boundary conditions are normally applied in the heat transfer surfaces of the PCM container. Such conditions include: (i) neglecting the effects of the container wall thickness, (ii) application of constant temperature or heat flux over time, and (iii) ignoring the impact of fluid flow outside of the container, which generates different values of heat flux in different positions. On a general basis, it is evident that fewer simplifications in the model imply more difficulties in the resolution and additional processing time. Computer and processor evolution that has been developed so far has allowed researchers to increase the complexity of the models by including more features as exhaustive properties, boundary conditions, physical conditions, and more detailed or different geometries. Examples of such features may consist of three-dimensional models, geometries with PCM and air taking into account PCM volume variation [52,99], and models considering the difference in heat flux outside the PCM container, according to HTF flow [100].

Formulations
The full phase change models are defined by the combination of one numerical method for latent heat in the energy conservation equation and one numerical method for velocity transition in the momentum equation. Among the abundant possible combinations of numerical methods, the most common are shown in Table 3, in which the latent heat and velocity transition formulations are combined. Several studies are based on the combination of E-EM with Darcy STM, while no literature was found that deals with the combinations of E-STM and STM, E-TTM and STM or E-EM and SOM. Besides the combinations of two single numerical formulations used for phase-change simulations, some authors have tried to adopt "mixed methods," i.e., the combinations of various formulations. Cao and Faghri [82] used E-TTM-AHC, as well as E-STM combined with VVM, for solidification and melting with different mushy zone ranges, reporting an increase in convergence difficulty as the mushy zone range reduces. Erek et al. [132] used E-TTM-AHC and E-STM for solidification outside a horizontal finned tube in which only thermal conduction was contemplated. Regarding models for metallic alloys, Bennon and Incropera [78,133] utilized E-EM with STM coupled with a specific formulation for the alloy composition and solubility.
Several authors have compared different combinations of numerical formulations referred to phase change processes. Most of the tests have been performed on rectangular cavities with adiabatic top and bottom surfaces, and heat transfer in one or both side surfaces. Such domains were tested in the following evaluations. Regarding latent heat formulations, Poirier and Salcudean [90] compared different results of latent heat for solidification: E-TTM-HI, E-TTM-AHC, E-TTM-EHC, E-EM, and a combination of E-EM and E-TTM-AHC. It was demonstrated that the suitability of each method depends on the type of phase change (isothermal or with continuous properties), as well as on the relationship between sensible and latent heat. Lamberg et al. [134] compared E-EM and E-TTM-EHC for solidification and melting, reporting slight nonconformities in the results for different temperature ranges used in E-TTM-EHC. Regarding the velocity transition formulation, Voller et al. [92] compared SOM, VVM, and Darcy STM, all combined with E-EM for solidification in a rectangular cavity. Implementation was considered more accessible with Darcy STM when compared to the other formulations. Ma and Zhang [32] performed a comparison of RSOM, RSTM, and VVM, all using the E-TTM latent heat formulation in a rectangular cavity. Ramped methods proved to improve the convergence if compared with non-ramped (SOM and STM) formulations, and the results of RSOM and RSTM were more accurate than VVM. Finally, Samara et al. [96] compared VVM and a mixed method of VVM and Darcy STM, both with E-TTM-AHC in a rectangular cavity: the mixed method had faster convergence than VVM.

Performance Evaluation of Computers and Numerical Models
The performance evaluation method shown in the present work can be applied to compare numerical models or different computers. Regarding geometrical characteristics, the analysis can be performed in different geometrical configurations, since they have the same number of dimensions (i.e., two-dimensional domains). It is worth mentioning that several characteristics of the computers are not included in the method, such as frequency, number of cores, or memory bandwidth. Such a method consists of two simple normalization steps based on the number of mesh elements and on the number of time steps evaluated. The first step consists of the mean number of time steps solved per time unit, through Equation (10), i.e., dt = n dt t e (10) where dt is the rate of time steps solved, n dt is the number of time steps, and t e is the computational time elapsed to solve n dt . The second normalization step is based on the grid size, aimed to obtain the rate of units solved (η), using Equation (11), i.e., η = n e dt (11) where n e is the number of elements of the grid. As a brief comparison of the evolution of processing capabilities of computers, we choose two studies from the 1980s, namely Bennon and Incropera [133] and Brent et al. [103], and two other studies dated after year 2000, namely Assis et al. [48] and Hosseinizadeh et al. [109]. Despite the greater complexity in the most recent studies, all are characterized by two-dimensional geometries and consider natural convection. The obtained outcomes are shown in Table 4.  Table 4 shows that the performance of modern computers to solve phase change problems is 1-2 orders of magnitude higher than that of supercomputers in the late 1980s and at least three orders of magnitude higher than that of regular computers in the late 1980s. However, it is worth mentioning that the considered parameter is not precisely a processing capacity indicator, as many relevant factors in the studies are different, such as geometries, numerical configurations, solution methods, and the number of iterations per time step.
The implementation of the numerical methods for solving melting and solidification problems in latent heat storage, as well as details about the accuracy, grid independence and validation codes using commercial software, as well as self-developed codes can be found in works [135,136]. It is worth emphasizing that the in-lab-scale or full-scale heat exchangers have also been discussed in the literature (see Zauner et al., 2016 [137]).

Conclusions
The present paper has focused on a literature review concerning numerical models of solidification and melting processes. Starting from the fundamental notions regarding PCM classification and application, the treatment was aimed at the analysis of the physical characteristics of solidification and melting processes. Regarding numerical models for phase change processes, the review describes the existing fixed grid methods, as well as compiles the combinations of various techniques used in PCM studies. Additionally, a new primary method to evaluate computer performance in numerical simulations has been proposed.
Phase change processes occur in a single temperature or over a temperature range. The solidification process is predominantly conductive, while in the melting process, convection heat transfer usually plays an essential role in defining the shape of the solid portion along the process. Close contact melting represents a specific case, as the motion of the solid phase exerts a significant effect in the form of the solid and increases the melting rate. The main characteristic in close contact melting is the formation of a skinny liquid layer between the solid and the heated surface. As reported in the literature, most of the melting rate occurs in the liquid layer. However, due to the small thickness, the liquid layer was hardly studied, given the difficulties to observe directly or even to measure such a region. In the growing scenario of PCM literature [5][6][7][138][139][140], further investigation efforts have to be made to clarify the physical phenomena that occur in the liquid layer, as well as to enforce and validate the existing theoretical equations.
Finally, a comparison of computational methods for phase change phenomena has been pointed out, to help researchers to choose the most suitable formulation for each specific case. First, the combination of E-EM and Darcy STM, also known as "enthalpy-porosity," have proved to be the most widespread in the literature. On the other hand, several combinations of methods were barely explored, such as E-STM or E-TTM associated with STM or RSTM, or even E-EM with SOM or RSOM. Therefore, such combinations, added to mixed methods, represent a large area to be explored in future research.
Author Contributions: This work is the result of collaborative research in which all the co-authors contributed equally to the paper.