Effects of Total Thermal Balance on the Thermal Energy Absorbed or Released by a High-Temperature Phase Change Material

Front tracking and enthalpy methods used to study phase change processes are based on a local thermal energy balance at the liquid–solid interface where mass accommodation methods are also used to account for the density change during the phase transition. Recently, it has been shown that a local thermal balance at the interface does not reproduce the thermodynamic equilibrium in adiabatic systems. Total thermal balance through the entire liquid–solid system can predict the correct thermodynamic equilibrium values of melted (solidified) mass, system size, and interface position. In this work, total thermal balance is applied to systems with isothermal–adiabatic boundary conditions to estimate the sensible and latent heat stored (released) by KNO3 and KNO3/NaNO3 salts which are used as high-temperature phase change materials. Relative percent differences between the solutions obtained with a local thermal balance at the interface and a total thermal balance for the thermal energy absorbed or released by high-temperature phase change materials are obtained. According to the total thermal balance proposed, a correction to the liquid–solid interface dynamics is introduced, which accounts for an extra amount of energy absorbed or released during the phase transition. It is shown that melting or solidification rates are modified by using a total thermal balance through the entire system. Finally, the numerical and semi-analytical methods illustrate that volume changes and the fraction of melted (solidified) solid (liquid) estimated through a local thermal balance at the interface are not invariant in adiabatic systems. The invariance of numerical and semi-analytical solutions in adiabatic systems is significantly improved through the proposed model.


Introduction
High-temperature phase change materials (HTPCMs) represent an appealing application as latent heat and sensible heat storage devices used as backup systems for thermoelectric generation [1,2]. Concentrating solar power (CSP) plants use solar power to generate superheated steam during the hours of the day with the highest solar irradiance. Latent heat and sensible heat storage devices are widely used to increase the renewable energy efficiency of CSP plants by providing thermal energy during the hours of the day with unusable solar power. On the one hand, extensive research has been performed to improve mathematical models that can predict the energy density and charging/discharging times in several configurations of thermal energy storage (TES) units [3][4][5]. On the other hand, experimental work has been performed to study the effect of device geometries [6] and composite based heat storage units [7][8][9] in order to enhance the heat transfer rates between the PCM and the heat transfer fluid (HTF). The latent heat and melting temperature of salts used as HTPCMs have been determined through molecular dynamics and ab initio simulations [10]. Mass accommodation methods have been developed to include the effects of density variations produced by pressure increments during melting of encapsulated PCMs [4,[11][12][13][14]. The pressure increments on spherically microencapsulated HTPCMs have been predicted through models that assume incompressible solid phases [4]. The latent heat energy density is diminished by pressure increments, and different models have been proposed to estimate the latent heat during melting of encapsulated PCMs [4,5,14]. The research has been focused on improving the desired thermodynamic properties of PCMs to enhance the heat energy transfer rates between the PCM and the HTF, and increase the thermal energy density of PCMs.
The phase change process related with PCMs lies within the moving boundary type of problems. Two phase one-dimensional problems have been extensively studied in the literature, and there are a few cases for which it is possible to find an exact analytical solution [15,16]. For this reason, different numerical and semi-analytical methods have been used to find approximate solutions to these type of problems [17,18]. The Finite Element Method (FEM) is one of the most popular numerical methods used to solve systems of differential equations. There is an extensive amount of research devoted to develop more accurate and efficient solutions to heat transfer related phenomena through the FEM, including the dynamics of first-order phase transitions [17]. Explicit FEMs have been proposed to find solutions for the temperature distribution in non-linear heat transfer problems, where the temperature dependence of the thermodynamic variables of specific heat and thermal conductivity were taken into consideration [19]. Modifications to the FEM have also been used to find more accurate solutions to the liquid-solid phase transition [17,20,21].
The non-linearity introduced by the thermal balance at the interface during a firstorder phase transition has also been studied through other numerical and semi-analytical methods. Explicit and implicit finite difference methods (FDMs) have also been used to find solutions for the one-dimensional liquid-solid phase transition at constant pressure [18,22,23]. The advantages of FEMs over FDMs rely on the way in which adiabatic boundary conditions are implemented and the interpolation methods needed to calculate the temperature in front tracking problems. Adiabatic systems have been used to find corrections to mathematical models of phase change processes, where total energy conservation plays a key role in these type of systems [23]. On the one hand, adiabatic boundary conditions do not need to be approximated when using FEMs, allowing well-behaved solutions at thermodynamic equilibrium. On the other hand, high-order approximations to the spatial derivative involved in this type of boundary conditions are needed when using FDMs, in order to reach the expected state of the system at thermodynamic equilibrium [23]. Semi-analytical approaches have also been used to find solutions for the liquid-solid interface dynamics through the heat balance integral method (HBIM) [23,24] and the refined heat balance integral method (RHBIM) [13,14,25]. Finally, the accuracy of the proposed solutions to the moving boundary problem when the system is subjected to different types of boundary conditions has been tested with the available similarity solutions, other numerical solutions, or experimental results [26][27][28].
This work is concerned with some of the fundamental aspects during liquid-solid isobaric phase transitions in one-dimensional configurations. The phase change process is held at constant pressure (isobaric phase transition), where volume changes produced by the density difference between liquid and solid phases are taken into consideration by imposing total mass conservation of the liquid-solid sample. The goal is to take a further step into the estimation of the energy capacity of HTPCMs and charging/discharging times. Recently, higher-order corrections to the liquid-solid interface dynamics have been introduced by proposing a total thermal balance through the entire system instead of the more classical local thermal balance at the interface [29]. These corrections are not intuitive and could be interpreted through an equivalent latent heat of fusion using the Leibniz integral rule. In this work, the effects of total thermal balance in the energy capacity of HTPCMs and charging/discharging times are estimated. The precision of the numerical and semi-analytical methods is relevant for estimating the higher-order contributions from total thermal balance in the thermal energy absorbed or released by a PCM sample. The RHBIM and FEM are used to find solutions to the proposed and classical models. The analysis of the thermal energy absorbed/released by the PCM indicates that the equivalent latent incorporates an extra or missing thermal energy contribution during the phase transition when compared with the predictions obtained from the classical model. The equivalent latent heat can be pictured as an apparent latent heat, since it introduces corrections to the thermal energy absorption/release rates, but the latent heat absorbed is shown to be related only with the bulk latent heat of fusion. The consistency of total thermal balance is determined through the invariance of the numerical and semianalytical solutions in adiabatic systems, whether the volume changes of the sample are incorporated by assuming a right or left moving boundary. Additionally, numerical and semi-analytical solutions to invariant quantities such as volume changes and the fraction of melted (solidified) solid (liquid) are found to reach the same thermodynamic equilibrium values in adiabatic (thermally isolated) systems. Finally, the contributions from total thermal balance are determined through the numerical and semi-analytical solutions to the total energy absorbed (released) during a charging (discharging) process in KNO 3 and KNO 3 /NaNO 3 salts.

Description of the Physical System
Consider a sample of size L(t) and cross-section A with a liquid phase in contact with a solid phase, separated by an interface with position ξ(t) at some melting temperature T m . The sample has a left (right) boundary at x (t) x s (t) , which obeys an equation of motion that imposes total mass conservation. The net flux of thermal energy through the interface causes its displacement in time t. The volume changes during the phase transition are incorporated through the dynamic variable x (t) x s (t) . The heat transfer through the system is homogeneous about the perpendicular plane to the heat flux. The temperature of the sample at some position x along the longitudinal axis and at some time t is considered to be uniform along the y − z plane. The thermal flux only takes place along the longitudinal direction x since there are no temperature gradients along the y − z plane. Under this consideration, the problem becomes one-dimensional, and the interface dynamics will only take place along the longitudinal x-direction. The thermodynamic variables of thermal conductivity and specific heat capacity are assumed to be independent of the temperature. Thermal expansion through temperature-dependent densities is also neglected; then, natural convection driven by buoyancy is not considered. The shrinkage or expansion of the sample, due to solidification or melting, is held at constant pressure. Natural convection induced by pressure gradients within the liquid phase can be neglected. The liquid ( ) and solid (s) phases have a temperature profile T i (x, t), where i = , s. Here, the temperature at any point within the liquid (solid) is above (below) the melting temperature T m , which is constant in time and does not depend on the spatial variable x.

Adiabatic Systems
The system is thermally isolated; therefore, it will be subjected to the following adiabatic boundary conditions: and an isothermal boundary condition at the interface where ξ(t) represents the position of the liquid-solid interface at some time t. Equation (1) means that the system is thermally isolated from the surroundings, which implies that thermal energy is only transferred between the liquid and solid phases. The boundary condition at the liquid-solid interface ξ(t), given by Equation (2), is used to describe the isothermal nature of the phase transition. The melting or solidification process will take place at a constant temperature equal to the melting temperature T m of the PCM. The initial temperature profiles will be assumed to be polynomial functions of x: where f i (x) is a polynomial function that represents the initial temperature distribution at phase i. The initial temperature profile can be obtained in order to satisfy the boundary conditions given by Equations (1) and (2). Initially, there is a certain amount of liquid and solid in the sample so that ξ(0) > 0.
Since the thermodynamic variables are independent of the temperature, the heat equation in each medium can be written as follows: where α i = k i /(ρ i C i ) is the heat diffusion coefficient of phase i, and is defined in terms of the thermal conductivity k i , specific heat C i and density ρ i . The thermal energy balance at the liquid-solid interface ξ(t) that is consistent with mass conservation is given by the following equations: , or (5a) . (5b) where L f = (C − C s )T m is the bulk latent heat of fusion, and C (C s ) is the specific heat capacity of the liquid (solid) close to the melting temperature T m . On the one hand, Equation (5a) is the equation of motion for the interface when x s (t) is chosen as a moving boundary and x (t) is constant in time. On the other hand, Equation (5b) is the equation of motion for ξ(t), when x (t) is chosen as a moving boundary and x s (t) is constant in time, as shown in Ref. [23]. Equations (5a) and (5b) result from a local energy balance at the the liquid-solid interface. The volume change experienced by the system due to the melting (solidification) of a small solid (liquid) slab is conceived by imposing total mass conservation of the system. Assuming that heat transfer is homogeneous in the direction perpendicular to the thermal flux, mass conservation can be imposed through the following equations of motion [23]: Equations (5a) and (6a) are used to estimate the amount of melted (solidified) solid (liquid) and the total volume change of the sample when x s (t) is chosen as the dynamical variable. Equations (5b) and (6b) are solved when x (t) is the dynamical variable [23]. The equation of motion for the left (right) boundary, coupled to the equation of motion for the liquid-solid interface position, has been able to reproduce the correct thermodynamic equilibrium values in some cases. Recently, it has been shown that the local thermal balance at the interface cannot predict thermodynamic equilibrium values for ξ eq , system size, and melted (solidified) mass in high-temperature materials [29]. Total thermal balance was proposed by forcing energy conservation in adiabatic systems. The energy and mass balance proposed resulted in a higher-order contribution to the equation of motion for the interface as follows: , or (7a) Equations (7a) or (7b) incorporate the corrections obtained from total thermal balance when x s (t) or x (t) is chosen as the dynamical variable. The above equations of motion take into account an extra amount of thermal energy transfer during the phase change process. The higher-order contribution obtained by imposing total energy conservation in adiabatic systems was pictured through an apparent latent heat [29], given by Using these definitions, the previous form of the equation of motion for ξ(t) shown through Equations (5a) or (5b) is recovered. Total thermal balance contributions are incorporated in the apparent latent heat defined through Equations (8a) or (8b). Total thermal balance contributions depend on the boundary conditions and relative density difference between both phases as shown through Equations (8a) and (8b). On the one hand, materials used for thermoelectric generation, such as salts with high melting temperatures, are subjected to wide temperature ranges, where the contributions from the apparent latent heat are significant. On the other hand, PCMs used in thermal shielding applications are subjected to operating temperatures close to the melting point of the material, and the effects of total thermal balance are practically negligible. Thermodynamic equilibrium values for the volume change of the system and fraction of melted (solidified) solid (liquid) have been previously obtained [23]. These equilibrium values have been used to identify inconsistencies in the predictions of models proposed by other authors. The volume change and fraction of melted (solidified) solid (liquid) can be obtained through an energy balance between the liquid and solid. The energy balance allows to find the amount of energy ∆ U that will be absorbed (released) by a mass of solid (liquid) ∆ M s (∆ M ) to produce the phase transition between an initial state and the stationary state. The energy used to melt (solidify) a given mass of solid (liquid) ∆ M s (∆ M ) is given by where T (x, 0) T s (x, 0) is the initial temperature profile in the liquid (solid) domain. Thermodynamic equilibrium values for the volume change and fraction of melted (solidified) solid (liquid) [23] are given by where f i is the fraction of melted solid f s eq or solidified liquid f eq at thermodynamic equilibrium, and M i (0) is the initial mass of phase i. According to energy conservation in adiabatic systems, the values at thermodynamic equilibrium should be independent of which boundary x (t) or x s (t) is chosen as the dynamical variable. Numerical and semi-analytical solutions will be used to illustrate that ∆ L and f s or f are independent of which boundary is set as the moving variable.

Absorbed/Released Thermal Energy
Numerical and semi-analytical solutions will be found for the dynamics of the phase transition during the charging or discharging process of a HTPCM. Mix boundary conditions can be used to produce full melting or solidification of the system. The solid phase will be completely melted through the following boundary conditions: Here the liquid domain is subjected to isothermal boundary conditions, where the left boundary is kept at the highest operating temperature T H , and the solid domain is thermally isolated. Total thermal balance through the entire system introduces a correction through the apparent latent heat [29], which depends on the type of boundary conditions and the relative density difference between the liquid and solid phase, as follows: On the one hand, Equation (13a) describes the apparent latent heat for a system where x s (t) is the moving boundary. On the other hand, Equation (13b) describes the apparent latent heat when x (t) is the dynamical variable of motion [29]. The discharging process of the PCM is emulated through the following boundary conditions: Here T C represents the lowest operating temperature during the process. The right boundary is in contact with a cold reservoir, which keeps the temperature constant at x s (t).
According to the proposed model in Ref. [29], the apparent latent heat for a solidification process is given by Equation (15a) predicts the apparent latent heat value when x s (t) is the moving boundary, and Equation (15b) is the apparent latent heat for a system where x (t) is chosen as the dynamical variable. Melting and solidification rates are modified by the contributions from total thermal balance according to Equations (13a) and (13b) and Equations (15a) and (15b) [29]. The time evolution of the amount of thermal energy absorbed or released during the charging and discharging processes can be obtained through the enthalpy difference between the initial state of the system and the transient state at some time t. The thermal energy absorbed or released is obtained as follows: The apparent latent heat given by Equations (13a) and (13b) and Equations (15a) and (15b) takes into account the extra energy transfer during the phase transition, which is hidden by applying a local thermal balance at the interface. The corrections introduced through the apparent latent heat may change the predicted amount of melted (solidified) solid (liquid) during the phase change process when compared to the estimated predictions obtained by using the bulk latent heat in Equations (5a) and (5b). However, the latent heat absorbed or released must be determined through the bulk latent heat of fusion L f , as will be shown in the following section.

Finite Element Method
The problem depends on space x and time t; therefore, both directions are approximated. In this section, we present a hybrid method to solve the problem, where the FEM is used to solve the spatial part, and an implicit finite difference method is used to approximate the time derivatives. The space and time dependence of the temperature distributions in the liquid and solid phases will be explicitly used in this part of the section to avoid confusion during the description of the FEM used in this work.

Space Discretization: Finite Element Method
The region Ω = {Ω ∪ Ω s } represents the domain where the problem must be solved, such that Ω and Ω s corresponds to the domain of the liquid and solid phase, respectively. For example, the regions Ω (Ω s ) will be discretized into n (m) cubic elements, in such a way that Here x q (q = 1, 2, · · ·, 3n + 1) and x s p (p = 1, 2, · · ·, 3m + 1) are the coordinates of the nodes at the liquid and solid phases, respectively. The FEM formulation is the same for both phases; therefore, the subscript i that appears in Equation (4) will be omitted.
Suppose that we look for the approximate solution T(x, t) at some time t within certain space of functions H with dimension 4; and Ψ 1 , Ψ 2 , · · ·, Ψ 4 constitutes a basis of H. Then, the temperature field T(x, t) can be written as follows: whereT (e) d (t) is the temperature at node d, associated with element e. Applying the weighted residual method and Galerkin formulation [30,31], Equation (4) becomes where v(x) represents the test function. Integrating (19) by parts, the following is obtained The integrals over Ω can be expressed as the sum of integrals over each element, and Equation (20) becomes where k = n or k = m, depending on the phase. Taking v(x) = Ψ T (x) where the T represents the transpose, and substituting Equation (18) into Equation (21) the following is obtained: where Here T (e) (t) is the time derivative of the temperature at each node. K (e) , M (e) , and F (e) represent the element stiffness matrix, element mass matrix, and the element flux vector, respectively.
After assembling the local matrices the following global representation of the problem is obtained: where K is the global stiffness matrix and M is the global mass matrix. On the other hand, F is the global load vector andT(t) = [T 1 (t),T 2 (t), · · ·,T q (t)] T is the global temperature vector and represents the temperature at each node [x i 1 , x i 2 , · · ·, x i q ] at time t, with i = (s) for the liquid (solid) phase. Dimension q is related to the elements number k = n(m) used to discretize Ω i and matches the total number of nodes at each phase.

Shape Functions
Once the region Ω i is discretized, the shape functions Ψ(x), which will be used to approximate the solution T(x, t), are selected. These functions are piecewise polynomials defined at each element. Lagrange functions are the most common type of shape functions in finite element analysis. In this work, we will use cubic Lagrange functions as shape functions for the FEM implementation. The cubic approach needs four points at every element e: one at each end of the interval and two more points equally spaced. So, for each element e, there are four equidistant nodes when using this approach. Four basic functions (N 1 (x), N 2 (x), N 3 (x), and N 4 (x)) are needed on each element e. So, K (e) and M (e) are 4 × 4 square matrices and F e = [F 1 , F 2 , F 3 , F 4 ] T . Dimension q = 3n + 1 (3m + 1). A higher degree of the shape function implies more exact results but also bigger dimensions for the system (26). Cubic Lagrange shape functions for element e are (23), (24) and (25) we can compute local matrices K (e) , M (e) and local vector F (e) : The integrals over each element are calculated by using Gaussian quadrature.

Time Discretization: Implicit Finite Difference Scheme
Through the FEM, the differential problem (4) has been transformed into a system of ordinary differential Equations (26). These equations must be solved for each time value t. The time derivative of the temperature (T(t)) that appears in Equation (26) will be approximated by using an implicit finite difference scheme. The temperatureT j represents the temperature distribution or global temperature vector at the jth time level. The time derivative of the temperature (T(t)) is approximated by using a backward difference in time as follows: where ∆t is the time step and j represents the jth time level. Substitution of Equation (31) in Equation (26) leads to the following system of equations The system of equations must be solved at each time level j. The temperature dis-tributionT j−1 corresponds to the temperature values at the jth−1 time level, which are determined in a previous step. The temperatures at each node are obtained by using the temperature values from the previous step, and the system of equations is solved in each phase.

Refined Heat Balance Integral Method
The equations of motion for the interface position, total mass conservation, and the heat equation will be solved through a cubic Lagrange implementation of the FEM and a refined heat balance integral method (RHBIM). The energy absorbed or released by the PCM will be determined through the numerical (FEM) and semi-analytical (RHBIM) solutions to the model based on a local energy balance at the interface and the proposed model that considers total thermal balance.
The RHBIM consists of proposing polynomial temperature profiles in each phase with time dependent coefficients that are found by solving the resulting set of ordinary differential equations in time. The method is described in Refs. [13,14], and a brief description of the semi-analytical method will be given in this part of the section. The liquid and solid domains are divided into m and n regions, respectively. Quadratic temperature profiles in the spatial variable x are proposed at each region. Continuity and smoothness of the temperature distributions between adjacent regions are imposed. The temperature profiles within each phase and region are given by Time-dependent coefficients a i (t), b i (t), and T i (t) at each region i within the liquid domain and c j (t), d j (t), and T j (t) at each region j within the solid domain are used to determine the time evolution of the temperature profiles. Continuity and smoothness, along with the isothermal boundary condition at the interface and the boundary conditions given by Equations (12) and (14), are used to find the time-dependent coefficients a i (t), c j (t), T i (t), and T j (t) in terms of b i (t) and d j (t). The interface position is found by substitution of the temperature profiles given through Equation (33) into the equation of motion for the interface, described through Equations (5a) or (5b). Total thermal balance is considered by using the corresponding expression for the apparent latent heat in Equations (5a) or (5b). Additionally, local thermal balance is considered by using the bulk latent heat in Equations (5a) or (5b). The position of the dynamical variable x s (t) or x (t) is found through total mass conservation given by Equations (6a) and (6b). Finally, the heat equation is averaged over the liquid and solid domains by using the temperature profiles given by Equation (33), resulting in a set of m + n ordinary differential equations (ODEs) in time for the coefficients b i (t) and d j (t). A first-order approximation for the time derivatives that appear in the resulting system of ODEs is used to obtain a linear system of algebraic equations for ξ(t), x s (t) x (t) , b i (t), and d j (t) at each region in the next time level t + ∆ t.

Results and Discussion
KNO 3 and the euctectic KNO 3 /NaNO 3 salts are used as HTPCMs to find the thermal energy absorbed (released) during a charging (discharging) process. The thermodynamic properties of the salts used are summarized in Table 1 and are the same properties used in Ref. [29]. The specific heat capacities correspond to their values close to the melting temperature of the PCM, and the latent heat of fusion is obtained as L f = (C − C s ) T m .

Invariance of Solutions in Adiabatic Systems
The first part of this section is devoted to the discussion of the phase change process in adiabatic systems. In this section, the HTPCM considered is the KNO 3 salt. A wide temperature range is used to highlight the breaking of the invariance when x s (t) or x (t) is chosen as the moving boundary, according to the solutions obtained through a local energy balance at the interface. Two examples are shown where the initial energy of the liquidsolid sample produces melting on the one hand and solidification on the other hand. The initial temperature at x (t) x s (t) is T H 0 = 923 K (T C 0 = 535 K) for the melting of KNO 3 with an initial interface position ξ(0) = 0.30 m and initial size of L(0) = 1.0 m. The initial temperature at x (t) x s (t) is T H 0 = 680 K (T C 0 = 300 K) for the solidification example, where ξ(0) = 0.70 m. These conditions are used to determine the initial temperature distribution in the liquid and solid domains as quadratic functions of the spatial variable x. The system is subjected to the adiabatic boundary conditions given by Equation (1) and the isothermal boundary condition at the liquid-solid interface, given by Equation (2). The total energy is a constant of the motion since the system is thermally isolated. The bar is expected to reach thermodynamic equilibrium where the system growth (shrinkage) and the fraction of melted (solidified) solid (liquid) is given by Equations (10) and (11), respectively. Figure 1 shows the time evolution of ∆ L and f s ( f ) for KNO 3 according to the solutions obtained by assuming the classical local thermal balance at the interface. Figure 1a,b shows the system growth and fraction of melted solid f s upon melting of KNO 3 . The numerical solutions to the classical model when x s (t) or x (t) is chosen as the moving boundary do not have the same time-dependent behavior for ∆ L and f s . These solutions are not invariant when the volume of the system is allowed to change from the right or left boundary. Additionally, it is realized that solutions for ∆ L and f s ( f ) reach entirely different thermodynamic equilibrium values when the sample is fixed at the left boundary and when the sample is fixed at the right boundary. The lack of invariance of the numerical solutions to the classical model is in contradiction with the predicted thermodynamic equilibrium values shown through Equations (10) and (11). Figure 1c,d also illustrates this anomalous behavior when the KNO 3 sample shrinks upon solidification of liquid in adiabatic systems. Equations (5a) or (5b) are solved through the FEM by imposing energy conservation through the apparent latent heat L * f given by Equations (8a) or (8b), whether x s (t) or x (t) is the time-dependent boundary. Additionally, total mass conservation is imposed through Equations (6a) or (6b). The numerical solutions to ∆ L and the fraction of melted (solidified) solid (liquid) are independent of which boundary is chosen as the dynamical variable, as illustrated in Figure 2. According to the proposed model, the FEM solutions are observed to reach the predicted thermodynamic equilibrium values given by Equations (10) and (11). The solutions during the melting and solidification processes in KNO 3 are independent of which boundary is chosen as the dynamical variable, as illustrated in Figure 2. Figure 2a  Total energy in adiabatic systems is a constant of the motion. The initial energy of the system E(0) must be conserved throughout the melting or solidification process. Therefore, energy conservation can be used to determine the performance of the numerical solutions. An average energy error (AEE) for the total energy has been defined as follows for this purpose: where E(0) is the initial energy of the system, E(t i ) is the total energy of the system at some time level t i , and n is the number of time partitions. Three versions of the FEM were implemented: the first one labeled as FEM 1 uses linear Lagrange functions, FEM 2 is an implementation of the FEM with quadratic Lagrange functions, and FEM 3 uses the cubic Lagrange shape functions described in Section 2.4. Table 2 shows the AEE for each implementation of the FEM during melting and solidification of KNO 3 . Results are shown for ten elements that were used to discretize the entire spatial domain. The number of elements within the liquid and solid domains was varied during the phase change process. The number of elements at each phase and in a given time level was determined from the volume of the liquid and solid phases in order to have a total number of ten elements in the whole system with the same length. A total of n = 1 × 10 5 time partitions were used on each of the examples illustrated in Figures 1 and 2 and Table 2. Increasing the size of n produces negligible changes to the AEE given by Equation (34) and shown in Table 2. Local matrices obtained for the FEM 1 , FEM 2 , and FEM 3 implementations have dimensions of 2, 3, and 4, respectively. Increasing the degree of the shape functions produces lower values for the AEE, as expected. Table 2 illustrates how total energy conservation is best behaved through a FEM 3 implementation of the FEM.  Table 2. Average energy error (AEE) in Giga-Joules (GJ), upon melting and solidification of KNO 3 in a sample with adiabatic boundary conditions according to the linear, quadratic, and cubic implementations of the FEM. The initial energy of the system during the solidification (melting) example is E(0) = 1.5692 GJ (E(0) = 1.6942 GJ). The AEEs are approximately three, four, and five orders of magnitude smaller than E(0) according to the FEM 1 , FEM 2 , and FEM 3 implementations, respectively. The AEE given by Equation (34) for each implementation of the FEM is illustrated when x s (t) or x (t) is the moving boundary.

Melting
Solidification

Energy Absorbed/Released
The liquid-solid system is subjected to mixed boundary conditions to produce melting on the one hand and solidification of the liquid phase on the other hand. Melting examples are designed through the boundary conditions given by Equation (12), and solidification of liquid phase is achieved by imposing the boundary conditions through Equation (14). The system will absorb (release) thermal energy until the solid (liquid) is almost completely melted (solidified). The total amount of energy absorbed or released by the PCM can be obtained through Equation (16a) or Equation (16b).
The sensible heat absorbed by a PCM can be conceived in four stages, as discussed in Refs. [14,25]. These stages can be used to show that latent heat is absorbed through the bulk latent heat L f of the PCM and not the apparent latent heat previously defined [29]. The first stage considers the thermal energy absorbed between t = 0 and some time t by the initial mass of liquid as where it is assumed that x (t) is constant and equal to x (t) = 0, and x s (t) is the dynamical variable used to impose total mass conservation. During a second stage, the amount of solid mass ∆ m s that will be eventually melted will increase its temperature from its initial value T s (x, 0) to the melting temperature T m , absorbing thermal energy as follows: where x p − ξ(0) is the volume of this amount of solid that will be melted between t = 0 and any later time t and is related with ∆ m s as ∆ m s = ρ s x p − ξ(0) . The value of x p can be obtained through mass conservation as follows: then, solving for x p the following expression is obtained . The third stage considers the energy absorbed by ∆ m s once it has transformed into liquid phase. Through this stage, this mass of liquid absorbs thermal energy from the melting temperature T m , to the temperature of the liquid phase T (x, t) at time t, as follows: where ξ(t) − ξ(0) is the volume of ∆ m s , but now in its liquid state. Finally, the mass of solid that was not melted between the initial state and the state of the system at some later time t, will absorb thermal energy as sensible heat by raising its temperature from its initial value T s (x, 0) to the temperature T s (x, t) at some later time t as follows: Adding the contributions from each stage to the total sensible heat absorbed, the following expression is obtained: The first two terms correspond to the total enthalpy absorbed by a melting PCM as shown by Equation (16a) for x (t) constant and equal to zero, and x s (t) as the moving boundary. The last term is exactly the latent heat absorbed by the PCM during the phase change process, which is proportional to the bulk latent heat L f = (C − C s ) T m and not the apparent latent heat. The above discussion shows that even though solidification or melting rates can be pictured through an apparent latent heat, the latent heat storage capacity of the PCM depends only on the bulk latent heat and not L * f . A similar analysis can be performed for solidification scenarios or when x (t) is chosen as the dynamical variable.
The present work estimates the contributions from the sensible and latent heat stored or released through Equations (16a) and (16b), according to the numerical and semi-analytical solutions to the old and new models. Figure 3 shows the total energy absorbed and released by the KNO 3 salt. Melting (solidification) examples are illustrated until the fraction of melted (solidified) solid (liquid) is f s ≥ 0.999 ( f ≥ 0.999). Numerical and semi-analytical solutions according to a local thermal balance at the interface and the proposed total thermal balance are shown. Melting of KNO 3 is produced through the boundary conditions given by Equation (12). The initial position of the liquid-solid interface is ξ(0) = 0.05 m in Figures 3a,b. The temperature at x (t) is kept at a constant value of T H = 680 K, and the initial temperature at the right boundary is T C 0 = 240 K. Figure 3a,b shows the numerical and semi-analytical estimations of the total energy absorbed by the salt according to both models and for each case, when x s (t) or x (t) is the moving boundary. The relative difference between the estimations of both models is practically negligible upon melting of the solid phase. Alternatively, significant differences are observed when the salt is releasing thermal energy. Figure 3c,d illustrates the numerical and semi-analytical solutions when the liquid experiences solidification and the system is subjected to the boundary conditions given by Equation (14). Initially, a small volume of solid is considered with ξ(0) = 0.95 m. The right boundary is kept at a constant temperature value of T C = 240 K and the initial temperature at x (t) is T H 0 = 680 K. The relative percent difference (RPD) between the old and new models in the total energy absorbed by the KNO 3 salt is shown in Table 3. The small difference can be understood in terms of the apparent latent heat, which depends on the boundary conditions. On the one hand, according to Equation (13a) the corrections to the bulk latent heat are proportional to the difference between the melting temperature and the temperature at the right boundary. The temperature at x s (t) is expected to approach T m as the solid melts; therefore, this term becomes smaller as the system evolves in time. On the other hand, the contributions from the apparent latent heat when x (t) is the moving boundary can be negligible according to Equation (13b) if the temperature at this boundary is close to T m . The RPD between the predictions obtained from the old and new models was determined as follows: where x i with i = , s corresponds to which boundary is considered as the dynamical variable.

FEM RHBIM
Days ∆h Higher RPD values are expected for the solidification of KNO 3 as illustrated in Table 4 since the temperature at the right boundary is much smaller than the melting point of the salt. According to Equation (15a), when the system releases energy and x s (t) is the moving boundary, the corrections introduced through the apparent latent heat are proportional to T m − T C . The example shown in Figure 3c illustrates the solidification of liquid at a very low temperature value of T C = 240 K, well below the melting point of the salt. However, when x (t) is the dynamical variable, the contributions from the apparent latent heat are much smaller since L * f approaches the bulk latent heat of fusion L f as the system evolves in time [29].
Finally, the energy absorbed and released by the KNO 3 /NaNO 3 salt is estimated through the numerical and semi-analytical solutions to both models discussed in this work. The RPD between both models is highly related to the thermodynamic properties of the material and the boundary conditions. The salt can be exposed to higher temperature values, and due to its lower melting temperature as shown in Table 1, the expected difference between both models should be higher when the salt is absorbing thermal energy. Figure 4 shows the numerical and semi-analytical solutions for the thermal energy absorbed and released by the KNO 3 /NaNO 3 salt. The charging process is shown in Figure 4a,b. The initial interface position during the melting process of the salt is 0.05 m. The temperature at x (t) is kept constant at T H = 866 K, well above the melting temperature T m = 496 K of the salt. The left boundary is thermally isolated, and its initial temperature is T C 0 = 300 K. The discharging process starts with a small volume of solid, where the initial interface position is ξ(0) = 0.95 m. Heat is removed from the right boundary which is kept at a constant temperature value of T C = 300 K, and the left boundary is initially set at T H 0 = 300 K. The solutions for the thermal energy released during the discharging process are shown in Figure 4c,d. Table 4. Released energy ∆ h = h(0) − h(t) during solidification of KNO 3 according to the numerical and semi-analytical solutions to the new and old models. The RPDs are obtained from the example shown in Figure 3c,d.

FEM RHBIM
Days ∆h   According to Equation (13b), due to the high difference between the melting temperature of the salt and a maximum charging temperature T H = 866 K, the contributions from the apparent latent heat are significant in the example shown in Figure 4b. The RPD for the energy absorbed by the salt and according to both models is shown in Table 5. The apparent latent heat, according to the new model, predicts lower energy absorption rates than those estimated by assuming local thermal balance at the interface. This behavior is observed in Figure 4b, which according to Equation (13b) the PCM should take longer time values to absorb thermal energy until the sample is completely melted. The asymptotic time behavior of the apparent latent heat when x s (t) is the moving boundary predicts lower RPD values, as illustrated in Figure 4a and Table 5. Table 5. RPD between both models during a charging process of KNO 3 /NaNO 3 salt and according to the numerical and semi-analytical solutions to the absorbed energy ∆ h = h(t) − h(0).

FEM RHBIM
Days ∆h  Lower RPD values are expected upon solidification of the KNO 3 /NaNO 3 salt since thermal energy is drained from the right boundary at temperatures close to the melting point of the salt. According to Equation (15a) the contributions from total thermal balance are significantly lower due to the operating temperature value at the right boundary. The effects of total thermal balance on the released energy are even smaller when x (t) is the moving boundary since the apparent latent heat approaches asymptotically to the bulk latent heat L f , as illustrated in Table 6.
Charging and discharging times are also an important parameter when HTPCMs are used as backup systems in thermoelectric generation applications. The proposed model estimates different energy densities of the salts used in this work and different charging or discharging times. According to the FEM and RHBIM solutions for the thermal energy absorbed or released by the KNO 3 salt, the highest RPD between both models is found during a discharging process. The result is consistent with Equation (15a) that predicts significantly smaller energy releasing rates. Then according to the new model, the PCM should release thermal energy at a faster rate. The behavior can be observed in Figure 3c. The maximum RPD between discharging times according to the numerical and semi-analytical solutions is 9.38 % and 9.17 %, respectively. Thermodynamic properties and operating temperatures used for the example shown in Figure 4 increase the RPD between both models in a charging process. According to the proposed model, the apparent latent heat is increased by total thermal balance, as shown through Equation (13b). Maximum RPD between charging times is expected, as shown in Figure 4b. Equation (13b) predicts lower energy absorption rates, as illustrated in Figure 4b. The estimated RPD between charging times according to the numerical and semi-analytical solutions in this example is 9.76 % and 9.58 %, respectively. Table 6. RPD between both models and according to the numerical and semi-analytical solutions to the released thermal energy during a solidification process of the KNO 3 /NaNO 3 salt.

Conclusions
Volume changes produced by the density difference between liquid and solid phases are taken into account through imposing total mass, as a constant of the motion. Conservation of total mass is incorporated through an equation of motion for the right or left boundary of the system, allowing the sample to expand or shrink during the phase transition. Thermodynamic equilibrium values for the volume changes and fraction of melted (solidified) solid (liquid) in adiabatic systems have been established and shown to be independent of the direction in which volume changes are taking place. The solutions obtained from the local thermal balance at the interface are not invariant in adiabatic systems and predict different thermodynamic equilibrium values, which is in contradiction with the expected behavior. Total energy conservation was imposed as a constant of the motion in adiabatic systems. Energy conservation introduced higher-order corrections to the equation of motion for the interface that were incorporated through an apparent latent heat. The numerical and semi-analytical solutions, considering energy transfer rates through the apparent latent heat, are invariant and consistent with previously established thermodynamic equilibrium values in adiabatic systems.
Thermal balance through the system was considered in samples with other types of boundary conditions. The corrections from total thermal balance were introduced by the apparent latent heat. These corrections are deeply related with the type of boundary conditions, the melting temperature of the material, and the relative density difference between liquid and solid phases. Significant differences between the local thermal balance at the interface and the total thermal balance in the energy absorbed (released) are found in HTPCMs. The magnitude of the correction introduced by total thermal balance to the thermal energy absorbed (released) is increased by the thermodynamic properties of HTPCMs and typical operating temperatures. The apparent latent heat only introduces a change to the energy absorption (release) rates, but it can not be conceived as a correction to the bulk latent heat of the PCM. Finally, charging (discharging) times are also modified by considering the proposed total thermal balance. Relative percent differences of 9.58% between both models in the amount of time needed to melt (solidify) the solid (liquid) phase are found.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflicts of interest.
Sample Availability: Not applicable.

Abbreviations
The following abbreviations and symbols are used in this manuscript: Initial mass of the system T C Temperature at the right boundary T H Temperature at the left boundary ∆ L eq Volume change at thermodynamic equilibrium f s eq Fraction of melted solid at thermodynamic equilibrium f eq Fraction of solidified liquid at thermodynamic equilibrium ∆ U Absorbed/released latent heat at thermodynamic equilibrium