A Numerical Simulation of an Experimental Melting Process of a Phase-Change Material without Convective Flows

Featured Application: Different simpliﬁed and low-computational-cost models for phase-change CFD simulation were compared with experimental data for geometries where convective ﬂows could be negligible. Abstract: The melting process of lauric acid in a square container heated from the top surface was numerically studied from an experimental case. Knowledge of this process is of special interest for computationally efﬁcient modeling systems, such as PCM-enhanced photovoltaic panels in horizontal positions or energy storage using PCM embedded on ﬂat surfaces. In these systems, the geometric arrangement of the PCM hinders the ﬂuid-phase movements through natural convection, which slows the melting process and can cause overheating in the ﬂuid phase. Using Ansys Fluent Software, three different approaches and two simulation methods, enthalpy-porosity and effective heat capacity, were developed for the numerical study. The results were compared with experimental measurements in a successful evaluation of the accuracy of computational ﬂuid dynamics simulations. It could be observed that the effective heat capacity method presented signiﬁcant advantages over the enthalpy-porosity method, since similar accuracy results were obtained, and a lower computational cost was required.


Introduction
Latent heat thermal energy storage (LHTES) has been proposed as a solution to improve the efficiency of thermal energy storage (TES) systems [1]. LHTES uses a phasechange material (PCM) for energy storage during its melting process (charging), which is retrieved in the solidification process (discharging). The energy storage density of an LHTES system is higher than systems that only store sensitive heat. Phase-change processes in pure substances are isothermal, thereby increasing the thermal inertia of LHTES [2]. The use of PCMs has been extended to numerous applications, such as thermal regulation of photovoltaic systems (PV-PCM) [3], heating and domestic hot water (DHW) production [4,5], and building air cooling and heating circuits [6,7], among other applications [8,9].
However, one of the main problems with the use of PCMs is their low thermal conductivity [10], making storage system design a highly important process. In practice, fins [11] and nanoparticles are employed to improve the heat transfer process of a PCM thereby enhancing its properties [12].
The orientation of the heat flow is a very relevant issue in the process of heat transmission to the PCM. Shokouhmand et al. [13] experimentally investigated the phase-change process of lauric acid heated sideways, and the influence of horizontal fins varying the inclination of the container [14]. In all cases, there were large convective flows in the system that favored the phase-change process.
A particular instance occurs when the heat flow goes down vertically from an upper surface, a condition that can take place in cases using horizontal PV-PCM panels and PCM The paper is structured as follows: in Section 1 there is an introduction. Following this, the experimental facility, the experimental procedure, and the experimental results are described in Section 2. The CFD models of the system are introduced in Section 3, and in Section 4 a summary of the EP and EHC simulation methods that were used are described. Results of the CFD simulations are discussed in Section 5, while the accuracy of the CFD simulation compared to the experimental data is included in Section 6. Finally, Section 7 summarizes the conclusions.

Experimental Section
Experimental testing of the phase-change process of lauric acid yielded results that were used to define the boundary conditions and validate the results of an equivalent CFD model of the system. All experimental tests were developed in an isolated tank in which the volume of PCM was placed and heated until the phase change occurred. Figure 1 shows the experimental facility that was designed for the experiments.
Appl. Sci. 2022, 12, x FOR PEER REVIEW were developed in an isolated tank in which the volume of P placed and heated until the phase change occurred. Figure 1 show perimental facility that was designed for the experiments. A square base volume with sides 220 mm × 220 mm and heigh of PCM were tested, with thicknesses that are considered sufficien panel applications [32]. The entire container was manufactured in ter resin (VER) with a thickness of 10 mm, making heat tran through the walls negligible. An electrical resistance adhered to an ium plate laid on the top of the container was used as a heater proportional-integral-derivative controller (PID) was used to reg temperature of the isothermal process (see Figure 1). The conta thermally insulated using extruded polystyrene (XPS) and K-FL 19 mm [33] to minimize atmospheric heat loss. Additionally, all te performed in an adiabatic chamber ( Figure 2). A square base volume with sides 220 mm × 220 mm and height 67 mm of PCM were tested, with thicknesses that are considered sufficient in solar panel applications [32]. The entire container was manufactured in vinyl ester resin (VER) with a thickness of 10 mm, making heat transmission through the walls negligible. An electrical resistance adhered to an aluminium plate laid on the top of the container was used as a heater, while a proportional-integral-derivative controller (PID) was used to regulate the temperature of the isothermal process (see Figure 1). The container was thermally insulated using extruded polystyrene (XPS) and K-FLEX ST of 19 mm [33] to minimize atmospheric heat loss. Additionally, all tests were performed in an adiabatic chamber ( Figure 2). A small pipe connecting the PCM to the outside of the tank was installed in one of the corners of the top enclosure. The pipe had a dual purpose: to purge the air at filling and to allow the PCM to expand due to the decreased density during the melting process.
Lauric acid (CAS 143-07-7), with a purity greater than 98% from the manufacturer Sigma-Aldrich (Saint Louis, United States) [34], was selected as the PCM. Its thermophysical properties are shown in Table 1. A total of 30 T-type thermocouples (uncertainty temperature, u(T) ± 0.5 °C) were located inside the PCM on the same vertical plane, arranged in an array, separated from the sidewalls, and protected to prevent corrosion [35]. As shown in Figure 1, the thermocouple array was not regular since the tank had a plane of symmetry passing through column A. The thermocouple distribution was selected to obtain a regular and equivalent matrix A small pipe connecting the PCM to the outside of the tank was installed in one of the corners of the top enclosure. The pipe had a dual purpose: to purge the air at filling and to allow the PCM to expand due to the decreased density during the melting process.

Properties
Solid Liquid A total of 30 T-type thermocouples (uncertainty temperature, u(T) ± 0.5 • C) were located inside the PCM on the same vertical plane, arranged in an array, separated from the sidewalls, and protected to prevent corrosion [35]. As shown in Figure 1, the thermocouple array was not regular since the tank had a plane of symmetry passing through column A. The thermocouple distribution was selected to obtain a regular and equivalent matrix by projecting columns B, D, and F onto the plane of symmetry. Together with the mentioned symmetry, this new equivalent matrix was used to simplify the simulated geometries. In Figure 1, red dots represent these thermocouples and Figure 2c shows the thermocouple array. Additionally, other thermocouples were placed for reference temperatures on the heated surface, bottom, container sides, and the outside. All data were Appl. Sci. 2022, 12, 3640 5 of 18 recorded every 15 s by a Campbell-Scientific CR1000 datalogger and a Campbell-Scientific AM16/32B multiplexer.

Experimental Procedure
The filling of the container took place with the PCM in a liquid state, preventing air accumulation in the upper area. Additionally, the container had to be correctly levelled to ensure that all heat transmission from the top surface was only in the vertical direction.
The test consisted of a heating process at a constant temperature of 67 • C. The heating was kept on until all thermocouples inside the PCM indicated a temperature higher than the melting point, ensuring the complete melting of the PCM around the thermocouple array. In this way, the total duration of the experiment exceeded 28 h.

Experimental Results
The experimental temperatures recorded at each depth are represented in Figure 3. All thermocouples placed at the same depth recorded similar temperatures, which meant that the presence of convective flows that increased heat transmission was ruled out. Therefore, the PCM temperature only varied with depth. The experimental value considered for each depth was the average of the recorded temperatures measured by the six thermocouples set at the same height. In addition, the absence of a horizontal temperature gradient made it possible to affirm that the thermal effects of the walls did not influence the measurement area.  The phase-change process is clearly shown in Figure 3 as a change in the trend of the temperature curve. The temperature remained constant until the PCM volume was liquid. The further the PCM was from the hot overhead surface, the longer the PCM phase-change process was. This was due to the low thermal conductivity of the PCM and the absence of convective flows in the liquid phase.

Numerical Model
Three CFD simulation alternatives for the phase-change process in a thermal storage system using PCM, in which heat is applied from the over- The phase-change process is clearly shown in Figure 3 as a change in the trend of the temperature curve. The temperature remained constant until the PCM volume was liquid. The further the PCM was from the hot overhead surface, the longer the PCM phase-change process was. This was due to the low thermal conductivity of the PCM and the absence of convective flows in the liquid phase.

Numerical Model
Three CFD simulation alternatives for the phase-change process in a thermal storage system using PCM, in which heat is applied from the overhead surface, were studied. The Appl. Sci. 2022, 12, 3640 6 of 18 models ranged from a complete and realistic model (approach A) to a simple one based on a solid PCM system with variable properties (approach C). Figure 4 shows the diagrams corresponding to each approach. Geometries, main variables, and thermal fluxes presented in each model are shown.

Approach A-Complete
The first model, approach A, described the geometry considering the elements of the PCM tank. Since the tank was symmetrical, only the l half was simulated. The values obtained in the simulations for columns D', and F' (Figure 4) were equivalent to the results measured experime tally for columns B, D, and F ( Figure 1). The temperature interface betwe the PCM and the hot surface had to be set at 67.15 °C to simulate the exp imental heating. The EP method was used for modeling the phase-chan process. Solid and liquid PCM are part of a single volume of the same m terial (fixed grid). The thermophysical properties of lauric acid for the so and liquid phase were assumed constant, according to the values in Tab 1. However, for the phase-change temperature range (361.65 K-321.35 K each property had a linear variation as a function of temperature, from t solid to the liquid value. The thermophysical properties of the other ma rials in the system were fixed, according to the values attached in Table   Table 2. Thermophysical properties of the materials of the container.

Approach A-Complete
The first model, approach A, described the geometry considering all the elements of the PCM tank. Since the tank was symmetrical, only the left half was simulated. The values obtained in the simulations for columns B', D', and F' (Figure 4) were equivalent to the results measured experimentally for columns B, D, and F ( Figure 1). The temperature interface between the PCM and the hot surface had to be set at 67.15 • C to simulate the experimental heating. The EP method was used for modeling the phase-change process. Solid and liquid PCM are part of a single volume of the same material (fixed grid). The thermophysical properties of lauric acid for the solid and liquid phase were assumed constant, according to the values in Table 1. However, for the phase-change temperature range (361.65 K-321.35 K), each property had a linear variation as a function of temperature, from the solid to the liquid value. The thermophysical properties of the other materials in the system were fixed, according to the values attached in Table 2.  [38,39] Heat losses due to the external convection of the container were modeled as a boundary condition in order to calculate the convection coefficient (h c ), obtained from the Nusselt number (Nu) from Equation (1): where l C is the characteristic length of the surface and k is the thermal conductivity of the material. Several experimental correlations were found to estimate natural convection within different geometries and systems. The implemented correlations, attached in Table 3, were defined in the Fluent software using a user-defined function (UDF). Table 3. Empirical correlations of the Nusselt number for natural convection on different surfaces. For the characteristic length, l C , l is the vertical length of the plate, As is the area of the plate, and p is the perimeter of the plate.

Arrangement
Characteristic Length, l C Ra Limits Nu Ref.
Vertical plate l 10 4 -10 9 10 9 -10 13 Nu = 0.59 Ra 1/4 Nu = 0.1 Ra 1/3 [40] Horizontal plate; top surface of a hot plate As/p 10 4 -10 7 10 7 -10 11 Horizontal plate; bottom surface of a hot plate As/p 10 5 -10 11 Nu = 0.27 Ra 1/4 [42] These expressions are a function of the dimensionless Rayleigh number (Ra), given by Equation (2): where g is gravity, Γ a is the coefficient of thermal expansion of the air, T s is the surface temperature, T ∞ is the temperature of the flow away from the wall, Pr is the Prandtl number, and ν a is the kinematic viscosity of the air. Γ a , v a , and Pr were calculated using the average temperatures of T s and T ∞ . The experimentally recorded temperature from each thermocouple was compared with the simulated temperature of a given point for the evaluation of the simulation results.

Approach B-Simple
In the second model, approach B, the influence of the lateral wall on the melting process was not considered, so the geometry could be simplified by employing a symmetry boundary condition on each vertical side ( Figure 4). The double symmetry implied that there was only unidirectional vertical heat flux. In addition, the temperature gradient was only upwards. Therefore, the array of temperature measurement points was reduced to a single column. Although the phase-change model used in this approach was the EP, the same as in approach A, the simulation was simpler because only the PCM volume was considered. In this approach, the heat flux was not set directly, but it was calculated from the wall temperature defined in the boundary conditions. Experimental temperatures of Figure 3 were used to set the top and bottom wall temperatures, which were defined using time-dependent UDFs.
The assumption for the thermophysical properties of lauric acid was the same as in approach A.
The results of this simulation were evaluated by comparing the average experimental temperature of each row of thermocouples with the simulation temperature calculated at a point at the same height as the row of thermocouples.

Approach C-Solid
The last model proposed in this work, approach C, evaluated the melting process of PCM by the EHC method. This model considered the PCM as a solid in which the phase change only manifested itself as a change in the thermophysical properties, and it could be used when there are unidirectional heat conduction and no convective flows. The specific heat of the PCM was defined by a UDF, where it was modified to consider the latent heat of phase change, while the other properties were defined the same as they were in the previous two approaches. Boundary conditions were kept unvaried form approach B. Vertical walls were set as symmetrical while the top and bottom walls were defined by a time-dependent temperature from the experimentally measured temperatures.
This model reduced the computational cost, which improved the calculation process and permitted higher time steps. The geometry used in this case was the same as approach B's geometry (Figure 4), but it used the EHC method as a phase-change simulation model. The results were evaluated in the same way as in approach B.

Computational Approach
In this section, the thermodynamic methods implemented in this work are described, as well as the studies on the simulation mesh and the fitting of the parameters used for the CFD simulation.

Enthalpy-Porosity Method (EP)
The EP phase-change model, used in approaches A and B, was implemented in Fluent by default for simulations of phase changes (solidification and melting model) [43]. This method solved the phase change by defining the PCM domain as a liquid, for which the solid material was equivalent to a liquid substance of very high viscosity [44]. A finite range of temperatures was established between which the phase change occurred, known as the mushy zone, in which the fluid was treated as a pseudo porous medium. T s and T l were the lower and the higher range temperatures of the mushy zone, according to Table 1. For the mushy zone, β is defined as the liquid fraction, a parameter that represents the fraction of the cell volume that is in the liquid state. When β = 0, the material of the cell was in the solid phase, and if β = 1, all volume was in the liquid phase For this method, the energy equation is defined as [23]: where the enthalpy of the material, H, was calculated as the sum of the sensible enthalpy, h, and the latent heat of the melting process, ∆H. The sensible enthalpy was defined as the sum of the reference enthalpy of the material, h re f , and the integral of the specific heat, C p , between the temperature, T, and the temperature of the reference point of the substance, T re f . On the other side, the enthalpy was equal to the liquid fraction times the latent heat of fusion, L.
In the end, the method used an iterative process between Equation (4) and the enthalpy equations to solve the temperatures.
The EP method added a new term to the Navier -Stokes momentum equation, which defines the motion of a fluid. S is a momentum sink, defined from the Carman-Koseny equation, which extinguishes the velocity of solids due to their high porosity. where ∈ is a smaller constant to prevent division by zero, A mush is the mushy zone constant, and → v p is the solid velocity due to the pulling of solidified material out of the domain. As in this case, the solid material was not being pulled from the domain, → v p = 0. The mushy zone constant regulated the velocity of the phase change: the higher this value, the faster the transition between phases.

Effective Heat Capacity Method (EHC)
The EHC method was used for the evaluation of the phase change through the change in the specific heat of a substance. This method was proposed by Budak et al. [45] and Hashemiand Sliepcevich [46]. To this end, the function of the specific heat concerning the temperature was modified to consider the latent fusion heat. This method has been used in other studies [22,[47][48][49], which had even better results than the EP model when the convective flows were not considered [23]. For the EHC method, the energy equation is defined as: The value of the specific heat, C p , is defined by a piecewise function according to Equation (10): where C p,s and C p,l are constant values of specific heat for the solid and the liquid phase, and T l and T s are the temperatures that delimit the range of the phase change according to Table 1. The value for the zone of the phase change, C p,e f f , was calculated corresponding to Equation (11), as has been proposed in previous work [50,51].
where L is the latent heat of fusion, T l and T s are the temperatures of the phase-change range, C p,s and C p,l are the values of specific heat for the solid and the liquid phase, and T is the temperature of the PCM.

Computational Procedure
All simulations were 2D models carried out in a transient regime using a pressurebased solver and SIMPLE algorithm. The solution spatial discretization was set by Fluent 18.2 by default. However, for the transient formulation, a second-order implicit discretization was selected.
The convergence criteria were set at 10 −10 for the energy equation and 10 −5 for the rest of the equations, with a maximum of 250 iterations per time step. In cases where convective flows were low or imperceptible, the time step could be greatly lengthened, reaching values of several minutes when implemented using the EHC method [22]. This study was executed with fixed time-step sizes: 3 s for approaches A and B and 15 s for approach C.
The viscous model was established as laminar because the Rayleigh number (Equation (2)) was always less than 10 8 . For the EP method, the parameter A mush , which determined the velocity of the melting process, must be defined. The recommended value in studies with lauric acid is 10 5 [52,53]; however, it is not critical when conduction dominates heat transmission [54].

Mesh Studies
A study of the influence on mesh size was performed on each previously described model, considering computation time and convergence criteria. In all cases, the geometry of the study permitted the creation of structured mesh composed of regular quadrilaterals.
The sizes for analysis ranged from 0.2 mm to 2 mm. Table 4 shows, for each approach, the size and number of elements of the different mesh analyzed, and the selected option. The mesh selected for approaches A and B was 0.5 mm since the use of a smaller mesh size did not improve the results. On the other hand, a thicker mesh needs a higher quantity of iterations to converge, lengthening the time for its solution. The mesh for approach C was insensitive to resizing. Being a solid system, convergence was stable and quickly achieved in all cases. Additionally, the time step was varied to values above 15 s, maintaining the stability and accuracy of the calculation. Therefore, the mesh that was selected and applied had the finest mesh to provide the best resolution in the results.

CFD Results
In all the cases under analysis, the main variables were the temperature and the liquid fraction of PCM, β, defining the phase-change process. However, other parameters were computed by FLUENT, such as the total energy provided for the phase change in the PCM, whose estimate is represented in Figure 5. There were no considerable differences between the different approaches, A, B, and C, showing that the simplifications applied did not increase the error. As expected in an isothermal heating process, as the PCM was heated, the absorbed heat flow decreased since the flow was directly related to the temperature gradient. A key validatory aspect of the proposed simplifications was the absence of convective flows in the liquid phase of the PCM. Fluid velocity was therefore analyzed throughout the simulation, registering very low values, smaller than 6 × 10 −6 m/s, in approach A. No convective flows were A key validatory aspect of the proposed simplifications was the absence of convective flows in the liquid phase of the PCM. Fluid velocity was therefore analyzed throughout the simulation, registering very low values, smaller than 6 × 10 −6 m/s, in approach A. No convective flows were observed in approach B. Therefore, the transmission of heat by convection in the system under study was practically non-existent.
Finally, the value of the liquid fraction, β, was represented in Figure 6. This value is useful in order to know the operation of the storage systems and is usually complex to measure in a real process. In the models where the EP method was implemented, it was observed that there were scarcely any differences. More than 70% of the volume of lauric acid was molten at the end of the simulated period.
was therefore analyzed throughout the simulation, registering very low values, smaller than 6 × 10 −6 m/s, in approach A. No convective flows were observed in approach B. Therefore, the transmission of heat by convection in the system under study was practically non-existent.
Finally, the value of the liquid fraction, β, was represented in Figure 6. This value is useful in order to know the operation of the storage systems and is usually complex to measure in a real process. In the models where the EP method was implemented, it was observed that there were scarcely any differences. More than 70% of the volume of lauric acid was molten at the end of the simulated period.

Discussion
The validation of the proposed CFD models was performed by comparing the simulated results with the experimental data in the

Discussion
The validation of the proposed CFD models was performed by comparing the simulated results with the experimental data in the thermocouple array, the hot surface, and the bottom surface of the container. Figures 7-9 show the differences between the experimental and simulated temperature values for each approach (A, B, and C) analyzed throughout the test and for each of the thermocouple probes in the experimental equipment. The root mean square error (RMSD) and the coefficient of determination (R 2 ) parametrical statistics were used for the validation of the simulation of the phase-change process, defined in Equations (12) and (13), and summarized in Table 5: whereŷ n is an experimental temperature value, y n a simulated temperature value, y the average of the experimental temperature values, and N is the total number of experimental values measured (approach A: 957 values; approach B and C: 210 values). The experimental values of approach B and C were lower than those of approach A as the thermocouple matrix was reduced to a single column.     When the models were globally compared, no large deviations were observed, although, at the beginning of the tests, when temperatures were low, all models generated lower temperature values than the experimental ones. Another critical aspect, the melting point temperature, was where   When the models were globally compared, no large deviations were observed, although, at the beginning of the tests, when temperatures were low, all models generated lower temperature values than the experimental ones. Another critical aspect, the melting point temperature, was where  When the models were globally compared, no large deviations were observed, although, at the beginning of the tests, when temperatures were low, all models generated lower temperature values than the experimental ones. Another critical aspect, the melting point temperature, was where greater deviations were also observed. The EP and the EHC models considered a range of temperatures during which the phase change occurred, as shown in Table 6, moving slightly away from the actual behavior of the phase change in a pure substance. All three approaches showed good performance in their simulations of the meltingprocess temperature profiles, obtaining R 2 results higher than 98%. The results of approach A showed the highest deviations in the results, with RMSD in temperatures up to 1 K. The error in the model was greater at the beginning of the test, in the solid phase of the PCM, and the differences between the experimental and simulated temperatures increased in positions furthest from the center of the container; a fact that revealed the existence of horizontal heat flows within the model that did not appear in the experimental test. Approach A was the most complex model under consideration and involved a high number of boundary parameters, which were difficult to evaluate experimentally and were a possible source of error in the simulation. Nevertheless, the results clearly represented the process that occurred in the system.
On the other hand, in approaches B and C, which had geometries simplified by vertical symmetries, the results were more accurate, with RMSE in temperatures lower than 1 K. Only the average temperature for each row is included in Figures 8 and 9. The phase change occurred at a different time for each row, depending on the distance from the hot surface. As in approach A, the largest errors were observed in the solid phase and decreased as the PCM melted. In approaches B and C, the differences between experimental and simulated temperature values in row 5, the closest to the hot surface, increased with time, showing the emergence of convective flows in the liquid phase, which had not been considered in these studies and noted in approach A.
Three simulation approaches showed the highest deviation from the experimental data at the lowest temperatures, around the first hour of the simulation. As seen in Figure 10, the three approaches calculated temperatures lower than the experimental ones. The error was reduced throughout the simulation, resulting in temperatures that were calculated slightly higher than experimentally observed values. The error might be due to close deviations in the properties of lauric acid at temperatures above the melting point. The phase-change time estimated by each approach was close to the experimental one, showing the validity of the models in the study of the melting process.
Three simulation approaches showed the highest deviation from the experimental data at the lowest temperatures, around the first hour of the simulation. As seen in Figure 10, the three approaches calculated temperatures lower than the experimental ones. The error was reduced throughout the simulation, resulting in temperatures that were calculated slightly higher than experimentally observed values. The error might be due to close deviations in the properties of lauric acid at temperatures above the melting point. The phase-change time estimated by each approach was close to the experimental one, showing the validity of the models in the study of the melting process.

Conclusions
In this study, the phase change in lauric acid that was used as PCM in an environment without convective flows was experimentally and numerically analyzed. The experimental melting process of a PCM in a rectangular container only heated from the upper surface was studied. During this process, an array of thirty thermocouple probes distributed in the PCM bulk were used to characterize the temperature of the lauric acid, for defining and validating the CFD simulations carried out in ANSYS Fluent. The shortage of convective flows made it possible to simplify the simulated geometry so that the computational cost was significantly reduced. The study showed that simulations in which the EP model was employed yielded good results with low RMSE temperature values (around 1 K), although the simpler EHC model had clear advantages in this case.

Conclusions
In this study, the phase change in lauric acid that was used as PCM in an environment without convective flows was experimentally and numerically analyzed. The experimental melting process of a PCM in a rectangular container only heated from the upper surface was studied. During this process, an array of thirty thermocouple probes distributed in the PCM bulk were used to characterize the temperature of the lauric acid, for defining and validating the CFD simulations carried out in ANSYS Fluent. The shortage of convective flows made it possible to simplify the simulated geometry so that the computational cost was significantly reduced. The study showed that simulations in which the EP model was employed yielded good results with low RMSE temperature values (around 1 K), although the simpler EHC model had clear advantages in this case.
The absence of convective flows implied that the liquid PCM had no movement, which meant that heat transmission was only by conduction, as for solids. The EHC method simplified the PCM as a solid, which entailed more efficient use of time and computational resources than the EP model, due to it not resolving the velocity equations. For the same geometry, both models showed a similar accuracy from comparisons between the numerical results and the experimental temperatures. Nevertheless, while the simulation of the PCM with the EP method took approximately 10 h to complete, the EHC method barely needed an hour. Another advantage of this method was that it could use high time steps and larger mesh sizes, which allowed the simulation of systems with longer temporal variations and more complex geometries.

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

A mush
Mushy zone constant (kg/s·m 3 ) As Area of the horizontal plate (m 2 ) C p Specific heat (J/kg·K) g Gravity (