Effects of Pressure-Induced Density Changes in the Thermal Energy Absorbed by a Micro-Encapsulated Phase-Change Material

Density changes produced by pressure increments during melting of a spherically confined phase-change material have an impact on the thermal energy absorbed by the heat storage unit. Several authors have assumed incompressible phases to estimate the volume change of the phase-change material and the thermal balance at the liquid–solid interface. This assumption simplifies the problem but neglects the contribution of density changes to the thermal energy absorbed. In this work, a thermal balance at the interface that depends on the rate of change of the densities and on the shape of the container is found by imposing total mass conservation. The rigidity of the container is tuned through the coupling constant of an array of springs surrounding the phase-change material. This way, the behavior of the system can be probed from the isobaric to the isochoric regimes. The sensible and latent heat absorbed during the melting process are obtained by solving the proposed model through numerical and semi-analytical methods. Comparing the predictions obtained through our model, it is found that even for moderate pressures, the absorbed thermal energy predicted by other authors can be significantly overestimated.


Introduction
Phase-change materials (PCMs) have provided an extensive line of research due to their appealing applications in areas related to renewable energy systems for the reduction of fossil fuel consumption. These materials are used to provide thermal comfort in homes and buildings by exploiting the isothermal nature of first-order phase transitions [1][2][3][4]. Thermal isolation is provided by using the latent heat of PCMs to absorb thermal energy during the day and releasing this energy during the night hours [5,6]. PCMs may be encapsulated in micro-spheres to enhance thermal transport [7], and may be mixed with the concrete wall in volumetric fractions that do not compromise the mechanical resistance of constructions [8,9]. High-temperature phase-change materials (HTPCMs) constitute another type of application, where PCMs with high fusion temperatures are used to store thermal energy in concentrating solar power (CSP) plants for thermoelectric energy generation [10]. The micro-encapsulated materials are used as a part of the thermal energy storage (TES) unit, which feeds the thermoelectric plant through a heat-transfer fluid (HTF) during those hours without sunlight [11]. Enhancement of thermal conductivities to reduce the charging and discharging times on HTPCMs has also been addressed by studying composite materials [12][13][14]. Finally, studies on other enhancement techniques to reduce the duration of the charging-discharging cycles for domestic heat recovery applications have been addressed as well [15].
When thermal performance of buildings is achieved by mixing the micro-encapsulated material within the concrete matrix, the PCM cannot expand freely during the melting process. PCMs are expected to expand while melting, since most of these materials have lower densities in their liquid state than in the solid state [15,16]. Volumetric expansion of HTPCMs is also constrained when encapsulated in metallic shells, since the elastic modulus of the shell is much higher than the bulk modulus of the PCM in its liquid and solid states [17]. The melting process of the PCM in all these types of applications takes place near the isochoric regime (constant volume). Thermo-mechanical models have been developed for encapsulated HTPCMs in spherical configurations [18][19][20][21]. These proposals constitute a first approach for explaining the behavior of PCMs when encapsulated by different materials. Other authors [22,23], study the one-dimensional liquid-solid phase transition, by coupling the PCM to a linear spring and assuming incompressible phases. This assumption neglects the effects of density variations on the behavior of the phase transition. Other authors avoid the effect of density changes during the phase transition by introducing an air void [21]. More recently, a proposal for incorporating density changes during the melting process has been developed for a one-dimensional plane problem [24]. Large differences between the time evolution of the thermodynamic variables in [18][19][20]24] are predicted at high pressures.
In this work, we study the effects of pressure-induced density changes, on the thermal energy absorbed by a micro-encapsulated HTPCM. As mentioned earlier, these type of PCMs are used as part of the heat storage units for thermoelectric generation during periods without sunlight. Therefore, the thermal storage capability of the PCM has a direct impact on the efficiency of CSP plants.
The thermo-mechanical model presented in this paper has been elaborated through the key idea of total and local mass conservation in a spherically confined PCM. Our proposal leads to an estimation of the sensible heat, latent heat, and total heat absorbed by the PCM that presents significant differences when compared to the predictions of other authors [18][19][20]22,23]. Through mass conservation, an extra term that depends on the shape of the container and the rate of density changes, will be shown to appear on the thermal balance at the interface. To the authors knowledge, this term has been overlooked or neglected by other authors [18][19][20]22,23]. We will show that depending on the mechanical properties of the container, this term represents a small perturbation, or it may have a significant contribution to the thermal energy absorbed by the PCM.
First, to validate the proposed solutions, an isochoric limit will be presented. In this limit, the container is completely rigid and through mass conservation, exact analytical expressions for the dynamic and thermodynamic variables will be found. An implicit finite difference method (FDM) and a refined heat balance integral method (RHBIM) will be developed to solve the proposed model. The numerical and semi-analytical solutions will be validated through the isochoric limit, by probing the phase transition from the isobaric (constant pressure) to the isochoric regimes. The thermo-mechanical models proposed by other authors [18][19][20]22,23] are well behaved close to the isobaric regime; however, it will be shown that these models are unable to reproduce the isochoric limit. Finally, the contributions from density changes to the thermal energy absorbed, will be determined by comparing the proposed solutions with the predictions from other authors. It will be shown that sensible heat must be conceived in four stages, where mass conservation plays a key role. For the chosen PCM, examples will be presented, where the absorbed thermal energy is overestimated by other authors. The difference in the absorbed thermal energy according to the numerical and semi-analytical solutions to both models is observed to grow significantly with the rigidity of the container.

Description of the Physical System
The system under study consists of an encapsulated PCM in a spherical shell of radius R(t) at any time t, where liquid and solid phases coexist, so that the liquid-solid interface at any time t is located at r = r(t). The boundary of the spherical configuration located at r = R(t) is attached to a set of springs as shown in Figure 1. The system is subjected to the following boundary conditions where λ s is the thermal conductivity of the solid phase and T (r, t) T s (r, t) is the temperature distribution in the liquid(solid) phase at any instant t. The temperature T H at the surface of the PCM is always above the fusion temperature T f (t). Then, the chosen boundary conditions will produce melting of the PCM so that r(t) moves inwards as shown in Figure 1. Within the temperature range that will be considered, thermal expansion effects can be neglected. This assumption implies that pressure will be distributed uniformly along the PCM. The stiffness of the spring array shown in Figure 1 will be tuned through the spring constant k s to study the behavior of the system.

Energy-Mass Balance at the Interface and Heat-Transfer Mechanism
In this part of the present section, the energy-mass balance equation (EMB) at the interface and the proposed model for heat transfer, will be presented for a melting process when the PCM is under mechanical stress. In [18][19][20], the solid phase is assumed to be incompressible. Then, the effects of density changes in the solid, caused by the pressure increment within the PCM can be neglected. In this work, the assumption of incompressible phases [22,23] will be relaxed by coupling the mass conservation of the liquid-solid system with the thermal balance at the interface. First, the melted mass of solid ∆M s over a small time interval ∆t, can be obtained by realizing that ∆M s is transformed into liquid. Additionally, if the total mass of the system is conserved: where M (t + ∆t) is the mass of liquid at t + ∆t, after the melted solid ∆M s (t) is added to the initial mass of liquid M (t). Therefore, the thermal energy needed to melt the mass of solid ∆M s within a small time interval ∆t is given by Here M (t) = 4 π ρ (t) R 3 (t) − r 3 (t) /3 has been used. The right-hand side of the last equation is equal to the net heat flux at r = r(t). Then, the following EMB equation at the interface is obtained: where λ (λ s ) is the thermal conductivity of the liquid(solid) phase. Imposing mass conservation to the entire PCM, another equation for the densities and dynamical variables of motion can be obtained as follows Using this expression, Equation (4) can be written in terms of the density of the solid as This equation is equivalent to the EMB Equation (4) as long as mass conservation is imposed. The pressure increment ∆P within the liquid is equal to the increment in pressure within the solid phase. Using the bulk modulus of each phase [24], a relation between the densities in both media is given by: where B (B s ) is the bulk modulus of the liquid(solid) phase. The pressure increment at any time t can be obtained from the liquid or solid deformation given by the last equation, as follows The deformation that each medium experiences during the phase transition is coupled to the elastic properties of the spring array surrounding the spherical configuration, which is assumed to obey Hooke's law. Then, the deformation of the liquid phase coupled to the elastic constant of the springs is given by: wherek s = k s N s R 0 which will be expressed in multiples of B , and R 0 is the initial radius of the PCM. N s is the concentration of springs over the surface of the system. Changing this parameter, will allow the probing of the behavior of the phase transition for different pressure regimes. Equations (4), (5), (7) and (9) or alternatively, Equations (5), (7) and (9) constitute a set of nonlinear equations for the dynamical variables r(t) and R(t), and the densities of each phase. The liquid-solid saturation line is obtained through a second order approximation in ∆P and ∆T f of the free energy change as described in [18]. In this approximation the fusion temperature is given by: where ∆T f = T f (t) − T f (0). K (K s ) is the compressibility of the liquid(solid) phase, which can be obtained as the inverse of the bulk modulus; α (α s ) is the liquid(solid) thermal expansion coefficient and C (C s ) is the specific heat capacity of the liquid(solid) phase. Then, the latent heat of fusion can be obtained as follows [18] The heat equation for any phase i = liquid(solid) that is consistent with local mass conservation [24,25] is given in spherical coordinates as Finally, the heat equation presented in [18][19][20], and given by: will be used to describe the heat transfer within each phase according to the theory presented in [18].

Absorbed Sensible and Latent Heat
Sensible heat is absorbed through four different stages that will be described briefly. Mass conservation also plays a key role when considering the thermal energy absorbed by the PCM. This process has been studied for isobaric phase transitions [25]; however, the mechanism by which sensible heat is absorbed in this case is completely different, since the liquid-solid saturation curve must be considered in the process. The heat that enters the PCM between t and t + ∆t is absorbed as sensible heat in the following stages: (a) The mass of solid that will not experience a phase transition between t and t + ∆t absorbs part of the heat by changing its temperature from T s (r, t) to T s (r, t + ∆t), (b) The mass of solid ∆M s absorbs heat before changing to its liquid form by raising its temperature from T s (r, t) to the fusion temperature T f (t), (c) Once transformed to its liquid form, ∆M s absorbs sensible heat by changing its temperature from T f (t + ∆t) to T (r, t + ∆t). At this point, the fusion temperature has changed according to Equation (10), since the inner pressure increases after the phase transition, (d) The original mass of liquid at time t absorbs heat by increasing its temperature from T (r, t) to T (r, t + ∆t).
For each of the processes described above, the absorbed sensible heat will be obtained through the internal energy change. During the first stage, the mass of solid that does not experience a phase transition within the time interval ∆t is equal to M s (t + ∆t) = M s (t) − ∆M s . The internal energy change experienced by this mass is given by: where r s is the radius at time t of the solid sphere with mass M s (t + ∆t) and the interface position r(t + ∆t) is the radius of this same mass of solid at time t + ∆t. Then, applying mass conservation to M s (t + ∆t) between t and t + ∆t, the value of r s can be obtained as During the second stage, the mass of solid ∆M s is distributed along a spherical shell that occupies a volume equal to ∆V s = 4 π r 3 (t) − r 3 s /3, before the phase transition. Therefore, the internal energy change experienced by this mass of solid is given by After the phase transition, the pressure has changed and the fusion temperature has increased to T f (t + ∆t). Therefore, the internal energy change experienced by the melted mass ∆M s is given by: where r − r(t + ∆t) is the thickness of a spherical shell with mass ∆M s , in its liquid form. The melted mass ∆M s = M s (t) − M s (t + ∆t) can be found in terms of ρ s (t) and ρ s (t + ∆t). Additionally, using the volume of the liquid shell associated with this mass, the value of r can be obtained as follows Finally, the mass of liquid present in the system at time t, experiences the following change of internal energy Adding the contributions from Equations (14), (16), (17) and (19), the absorbed sensible heat is obtained. The energy absorbed as latent heat within the time interval ∆t is given by The total heat absorbed by the PCM is obtained by adding the results from Equations (14), (16), (17), (19) and (20) as follows

Initial Conditions
The results obtained from the numerical and semi-analytical methods applied to the proposed model, will be presented in two parts. First, we will obtain the pressure, density, and liquid-solid saturation line for the nitrate salt considered [18], in a wide range of pressures where the approximation established by Equations (10) and (11) is valid. Initially, the spherical configuration will be almost in its solid state with a very thin liquid layer surrounding the solid, and the melting process will be studied until practically all the solid has melted. Next, we will present the results for the sensible heat, latent heat, and total energy absorbed by the salt, for different values of the spring constant. The absorbed thermal energy obtained with the proposed model, will be compared with the results predicted by the thermo-mechanical model presented in [18].
All the results that will be presented in this section correspond to a PCM confined in a spherical configuration with an initial radius of R(0) = 1.0 mm. The solid phase is centered at the origin and has an initial radius of r(0) = 0.99 mm. The solid is surrounded by a thin liquid layer of R(0) − r(0) = 0.001 mm of thickness. Initially, the inner pressure of the system is P(0) = 1 atm. The temperature at the surface of the PCM is fixed at T H = 550 K. The net heat flux vanishes at the origin of the configuration according to Equation (1), and initially the temperature at r = 0 is equal to T C = 420 K. Large values of the bulk modulus of the solid phase B s will be assumed, to compare the predicted results from the proposed model, with the results from the model introduced in [18]. The thermodynamic variables at 1 atm of pressure are given in Table 1. Pressure-induced changes in the thermal conductivities, specific heat capacities, bulk modulus, and thermal expansion coefficients are neglected. Table 1. Thermodynamic data for the nitrite salt KNO 3 /NaNO 3 at 1 atm [18].

State Property Salt
Liquid

Numerical and Semi-Analytical Results: Isobaric and Isochoric Regimes
The authors in [18] assume an incompressible solid. In this approximation, the proposed thermal balance at the interface is given by [18]: which neglects the form factor given by Equation (6). The model presented in [18], estimates the volume change of the PCM in terms of the volumetric fraction of melted solid as follows: where is the volumetric fraction of melted solid. In this equation, the changes in the liquid and solid densities are neglected; therefore, it will only be valid when the transition takes place fork s B , which is close to the isobaric regime. Then, the thermal balance at the interface given by Equation (22), total volume expansion described by Equation (23) and Equations (9)-(11) and (13) represent the thermo-mechanical model proposed in [18][19][20].
The solutions to the model introduced in [18] and the proposed model in this work, will be compared for several values ofk s . Although the model proposed in this work does not need to assume an incompressible solid, B s B and B s k s will be assumed, to compare our results with the model proposed in [18]. For each value of the spring constant, the melting process will be studied until the mass fraction of melted solid is The density of the liquid, the inner pressure, the fusion temperature and latent heat of fusion will be obtained from both models at full melting fractions and for several values ofk s . An isochoric limit will be derived to establish the validity of the solutions to each model. In this limitk s B , and the volume of the spherical configuration does not change during the phase transition. At full melting, it is possible to obtain the density of the liquid for an isochoric transition by using mass conservation. When the total mass is conserved, the initial mass of the PCM M PCM (0) is equal to the total mass at full melting M PCM (t max ). Additionally, when the solid melts at constant volume, R(0) = R(t max ). Then, the maximum possible value for the density of the liquid ρ * = ρ (t max ) in this limit is given by The inner pressure in this limit is obtained by substituting Equation (24) in Equation (8); therefore, the isochoric limit for the inner pressure is given by The isochoric values at full melting fractions for the fusion temperature and latent heat of fusion are obtained by substituting ∆P * = P * − P(0) in Equations (10) and (11).  The results shown in Figure 2 confirm that the model proposed in this work is consistent with the isochoric limit and is well behaved near this regime. Since Equations (22) and (23) are valid for weak couplings or small values ofk s , where density changes can be neglected, the approximation presented in [18] is observed to be in good agreement with the model proposed in this work. In this regime, the form factor in Equation (4) and density variations are just a small perturbation to the solution of the proposed model. However, as illustrated in Figure 2, when the spring constantk s is increased and the phase transition deviates from the isobaric regime, the solutions to the model proposed in [18] start to diverge from the isochoric limit given by Equations (24) and (25).
In Figure 3, the results for L f and T f at melted solid fractions of f s = 0.999 are shown. The temperature of fusion and latent heat of fusion are very important thermodynamic variables when estimating the total heat absorbed by the PCM during the melting process. Therefore, the solutions obtained from the proposed model for these thermodynamic variables are also validated through the isochoric limit found by substitution of ∆P * = P * − P(0) in Equations (10) and (11). As illustrated in Figure 3, the numerical and semi-analytical solutions to the model proposed in this work, also show a good agreement with the liquid-solid saturation point in this limit. Also as expected, the solutions to the thermo-mechanical model presented in [18] are well behaved in the isobaric regime but tend to diverge from the isochoric limit as the spring constantk s is increased. Equation (23) is a statement of mass conservation for isobaric transitions, and it can be obtained from Equation (5), when density variations are neglected. The interpretation of the large differences between the two models for P ≥ 50 MPa as observed in Figure 2 can be inferred from the principle of mass conservation. Equation (5) can be expressed in terms of the total volume change experienced by the PCM during the melting process. Assuming that ρ s does not change with pressure increments and keeping only first-order terms in ∆ρ = ρ (t) − ρ (0), the volume change according to Equation (5) is given by The first term on this equation is exactly the volume expansion used in [18], as expressed through Equation (23). The second term represents the contribution to ∆V from the change of density experienced by the liquid phase. On the one hand, it is expected that for negligible density changes, both models show a good agreement fork s B , as illustrated in Figure 2. On the other hand, for moderate and high values ofk s , the contribution to ∆V from the change of density ∆ρ is not negligible. Within this regime, it is expected that the volume change according to Equation (23), leads to solutions with significant deviations from mass conservation. The mass of the PCM was obtained for f s = 0.999 and every value ofk s considered in Figures 2 and 3. According to the model proposed in [18], additional mass is created at moderate and large values ofk s as observed in Figure 4. The total mass values shown in Figure 4 were calculated through the RHBIM by solving the model proposed in [18] and shown in red squares, and the proposed model in this work, whose solutions are shown in blue squares. Each symbol corresponds to the total mass of the PCM at melting fractions of f s = 0.999 and for every value of the elastic constantk s considered in Figures 2 and 3. Finally, the isochoric limit established through Equation (24) has been used to obtain the maximum pressure increment ∆P max where the approximation to the liquid-solid saturation line is valid. According to Equation (24), if the initial radius of the solid sphere r(0) is equal to R(0) (the PCM is initially in its solid phase), ρ * = ρ s (0) in the isochoric limit. This establishes an upper bound for ρ at full melting. On the one hand, at moderate values ofk s = 3.5 B , the semi-analytical solution to the model introduced in [18] for ρ at melted solid fractions of f s = 0.999 is ρ = 2205.77 kg/m 3 , where the corresponding pressure is P max = 267.67 MPa. This value is above the established boundary of ρ * = 2192.00 kg/m 3 . On the other hand, the solution fork s = 3.0 B at f s = 0.999 is ρ = 2189.75 kg/m 3 , which lies below ρ * = ρ s (0). The corresponding pressure at this value ofk s is P max = 230.27 MPa. This criterion was used to establish a maximum value ofk s , above which the thermo-mechanical model introduced in [18] has no applicability.
From the above discussion it follows that for the type of PCM studied in this work, we cannot consider pressure increments above ∆P max = 230.1687 MPa. Additionally, at these pressure increments, the first-order terms in ∆P that appear in the approximation for the liquid-solid saturation temperature given by Equation (10) still dominate over the second order terms. Then, we will assume that the approximation given by Equation (10) is still valid for ∆P max = 230.1687 MPa. According to the numerical and semi-analytical results shown in Figure 2, the proposed model in this work is well behaved within the pressure regime where the approximation to the liquid-solid saturation line is valid. Even at high values ofk s = 1000 B , the numerical and semi-analytical solutions to our model, approach asymptotically to the isochoric limit given by Equations (24) and (25).

Numerical and Semi-Analytical Results: Absorbed Thermal Energy
In this section, the results obtained for the sensible and latent heat absorbed by the PCM will be presented. The thermal energy absorbed during the melting process ∆Q is obtained through the stages previously discussed. For large values ofk s , the pressure increments obtained through the model presented in [18] lie outside the pressure domain, where the approximation given by Equation (10) is valid. Therefore, to compare the results of this work, with the model introduced in [18], only weak and moderate values of the spring constant will be considered.
In Figure 5, the sensible heat, latent heat, and total thermal energy absorbed by the PCM is shown for small values ofk s where the transition is within the isobaric regime. The solutions obtained through the numerical and semi-analytical methods applied to both models are used to estimate the total energy according to Equations (14), (16), (17) and (19)  As illustrated in Figure 2, the model presented in [18] estimates larger pressure increments for moderate values ofk s , when compared to the proposed model. According to this result, it is expected that the latent heat of fusion drops faster at these values ofk s than the thermo-mechanical model proposed in this work. This behavior is observed in Figure 3. Therefore, the solution to the model introduced in [18] is expected to underestimate the latent heat absorbed by the nitrate salt for moderate values ofk s , when compared to the proposed model. This effect is more evident in Figure 6, and fork s = 3 B , where the FDM solution to ∆Q o f according to the model proposed in [18] (plotted as a continuous red line) lies slightly below the FDM solution to ∆Q p f obtained from the proposed model (plotted as a black dotted line). Even though, the relative percent difference (RPD) between the inner pressures predicted by each model for f s = 0.999 andk s = 3 B is practically 100%, the corresponding RPD in the latent heat absorbed is 3.36%. Then, for this type of PCM, the difference between in ∆Q f between both models may not be significant.
Although there is a small difference in the predictions for the latent heat absorbed by the PCM, large differences are observed in the sensible heat absorbed, as illustrated in Figure 6. Additionally, large differences in the maximum pressures are observed, since the phase transition for the values of k s considered in this figure is well outside the isobaric regime. The difference in ∆U between both models is related to mass conservation. In Figure 4, it is shown that the solution to the model proposed in [18] tends to create mass. This effect is more evident for higher values ofk s . Then, compared to the model proposed in this work, more energy is stored as sensible heat according to the model introduced in [18]. In this case, to produce a given increment in the temperature of the PCM requires more energy, because there is more mass to store sensible heat. The expected behavior is shown in Figure 6, where the absorbed sensible heat according to the model proposed in [18] (plotted as a continuous blue line) overestimates the absorbed sensible heat obtained with the solutions of the model proposed in this work (plotted with a dotted brown line). In Tables 2 and 3, the RPD for ∆Q is shown. The RPD was obtained from the numerical and semi-analytical solutions for the total energy absorbed, and it is defined as: where ∆Q (o) corresponds to the thermal energy absorbed according to other authors [18][19][20], and ∆Q (p) is the corresponding solution for the thermal energy according to the model proposed in this work. The RPD and the values for the absorbed energy at different melted solid fractions are shown in Tables 2 and 3. For small values of the spring constantk s , the solid melts within the isobaric regime as illustrated in Figure 5. Within this domain of pressures, we expect to observe small values of the RPD, as shown in Table 2.   Table 3 shows the RPD in the energy absorbed, which is illustrated in Figure 6. The RPD between the two models according to the RHBIM and the FDM is calculated at several melted solid fractions f s . According to the results illustrated in Figures 2 and 3, significant differences between ∆Q (o) and ∆Q (p) are expected in this pressure domain. Table 3. Calculated values of ∆Q fork s = B andk s = 3 B .

Numerical and Semi-Analytical Methods
The proposed model described by Equations (4), (5), (7) and (9)-(12) is solved through a second order FDM and the semi-analytical RHBIM. Also, the thermal balance at the interface and volume expansion given by Equations (22) and (23), and Equations (9)-(11) and (13) which represent the model proposed in [18][19][20] are solved through the numerical and semi-analytical methods. All these methods were implemented in our own Fortran codes and Maple for mathematical manipulation of long expressions. Figure 1 was created with the Inkskape drawing software and graphs shown in Figures 2-6, were produced by using the OriginLab software.
The initial temperature profile considered is a second order polynomial in r that satisfies the boundary conditions given by Equation (1), and the following initial conditions The initial temperature distribution in the liquid and solid phase is given by the following second order polynomials Here, the coefficients at phase i = (s) liquid(solid), a i and b i can be determined by applying the boundary conditions and initial conditions given by Equations (1) and (28), respectively. The fusion temperature T f (0), corresponds to the value of the liquid-solid saturation temperature at 1 atm of pressure.
After performing the time discretization in both methods, Equations (4), (5), (7) and (9) constitute a set of nonlinear equations for the dynamical variables of motion r(t) and R(t), and the densities of each phase. The time discretized EMB equation at the interface is then, given by: where n r n+1 r , n R n+1 R and n ρ n+1 ρ are the interface position, radius of the sphere, and liquid density at the nth(nth + 1) time level. Also, in Equation (30), the net heat flux per unit area at the interface has been defined aṡ The heat fluxq n r is obtained by using a fourth order approximation in ∆r to the spatial derivative when using the FDM [26]. The semi-analytical approach is to estimateq n r , through continuous temperature distributions [26]. The equation for mass conservation (5) is discretized as follows: dM PCM (t) dt = n r 2 n ρ s n+1 r− n r ∆t where n ρ s n+1 ρ s is the density of the solid phase at the nth(nth + 1) time level. Finally, the deformation of each phase (7), and the coupling between the liquid and spring array deformation (9) are discretized as: and The same discretization is used to solve the model introduced in [18]. In this case, the dynamic variables and the density of the liquid are uncoupled, since density changes are neglected in the thermal balance at the interface and in the estimation of ∆V. Then, the interface position can be solved directly from Equation (22)  From the volume change given by Equation (23), the radius of the PCM is given by: where the interface position at the nth + 1 time level can be substituted from Equation (35). Equation (33) is not necessary, given that in this approximation the solid is assumed to be incompressible. Finally, the density of the liquid phase can be solved from Equation (34).
The solutions for n r, n R, n ρ and n ρ s are used to obtain the pressure increment n ∆P = n P − P atm at the nth time level. Finally, the estimated pressure increment according to each model is used to obtain the fusion temperature n T f and the latent heat of fusion n L f from Equations (10) and (11), respectively.

FDM
The difference between the numerical and semi-analytical methods used in this work, lies in the way in which the temperature field is estimated. The spatial derivatives that appear in Equations (12) and (13) are approximated up to second order terms in ∆r. An implicit finite difference scheme is used [26], with a backward difference definition for the time derivatives that appear in Equations (12) and (13). The liquid(solid) phase domain is discretized by using a mesh with a constant number of nodes M 1 + 1(M 2 + 1). Additionally, the mesh is moved according to the dynamics of the interface and the system size [26]. A central definition for the first and second order spatial derivatives is used, so a system of M 1 (M 2 ) equations for the temperature n+1,m T ( n+1,m T s ) at each node is obtained [27]. The mesh size for the FDM used in this work was chosen with 161(161) nodes in the liquid(solid) phase. For a solid sphere with an initial radius of r(0) = 0.99 mm and a liquid shell with a thickness of 0.001 mm, the initial spatial step ∆r in the solid(liquid) is 0 ∆r s = 6.19 µm ( 0 ∆r = 0.0063 µm). The time step used in the simulations was ∆t = 0.5 µs.

RHBIM
The RHBIM used in this work has been fully described elsewhere [24,25,28,29] for one-dimensional problems in rectangular coordinates. However, a brief description of the methodology used in this work will be given. The solid sphere and liquid shell are divided in M and N regions, respectively. Continuous and smooth functions for the temperature distributions are proposed in each phase. The accuracy of the method depends on the type of temperature profile proposed, which classically is a quadratic polynomial [28,30]. However, depending on the nature of the boundary conditions and the geometry of the system, other profiles have been used [28,30]. In this work, parabolic profiles have been proposed within each region of the solid and liquid phase, as follows m T s (r, t) = m a s (t) r − m r s (t) + m b s (t) r − m r s (t) 2 + m T ( t), with 1 ≤ m ≤ M and n T (r, t) = n a (t) r − n r (t) + n b (t) r − n r (t) 2 + n T(t), with 1 ≤ n ≤ N.
Here the time dependent coefficients m a s (t) and n a (t) within each region of the solid and liquid phases, and the temperature values m T(t) and n T(t) for each of the solid and liquid regions, can be expressed in terms of the coefficients m b s (t) and n b (t). This can be done through the boundary conditions given by Equation (1) and by imposing continuity and smoothness at each boundary between adjacent regions, m r s (t) and n r (t) [29]. These boundaries are coupled to the dynamical behavior of the interface and radius of the PCM [24,25]; therefore, the time dependence of the boundaries between adjacent regions is given by: for the solid sphere, and m r (t) = r(t) for the liquid shell.
The temperature distributions given by Equation (37) are used to average Equations (12) and (13) over the spatial variable, as described in [25]. After the averaging process, a system of M(N) ordinary differential equations (ODEs) in time is found for the time dependent coefficients, m b s (t)( n b (t)). The resulting system of ODEs is discretized by using a forward difference approximation for the time derivative. Finally, each coefficient can be found for the next time level by solving the resulting system of algebraic equations. In this work, Equations (12) and (13) were averaged over 3 regions in each phase, resulting in 3 ODEs for m b s (t) and 3 ODEs for n b (t). The nonlinear system of equations described at the beginning of this section, and the linear system of equations for the time dependent coefficients m b s (t) and n b (t), were solved by using a time step of ∆t = 0.5 µs until the system reached melted solid fractions of f s = 0.999.

Conclusions
Several results can be outlined from this work. From the theoretical and physical points of view, mass conservation constitutes the central idea for estimating the contribution of density changes to the absorbed energy. Total mass conservation leads to an EMB equation at the interface, which for a spherical configuration is found to have an extra term that depends on the geometry of the system and rate of density changes. This extra term represents a small perturbation, when the melting process takes place near the isobaric regime. Within this pressure domain, the proposed solutions in this work, show a good agreement with the solutions from other authors. At high pressures, the contributions from density changes in the EMB at the interface and mass conservation equation are significant. Large differences between the proposed solutions to the densities and liquid-solid saturation line, and the predictions from other authors are observed in this pressure domain. Isochoric limits for the density, pressure, and liquid-solid saturation point, were obtained by assuming that total mass is also conserved in this regime. The solutions obtained through the proposed model are observed to be well behaved near the isochoric regime, and found to approach asymptotically to the predicted limit.
The most important finding in this work, corresponds to the effects of density changes in the thermal energy absorbed by the phase-change material. The sensible heat absorbed during the melting process was conceived through total and local mass conservation. According to the numerical and semi-analytical solutions to the model proposed in this work, the energy absorbed by the nitrate salt can be significantly overestimated by the solutions from other authors. On the one hand, the sensible and latent heat absorbed near the isobaric regime, show a good agreement between both models. On the other hand, at large pressures, where the approximation used for the liquid-solid saturation line is still valid, significant differences in the absorbed energy were observed. Within this pressure domain, and according to the numerical and semi-analytical solutions, an approximate difference of 16% between the predictions from both models, can be expected. This difference in the estimation of the heat storage capacity of TES units has a direct impact on the design and efficient use of renewable energy in CSP plants and domestic house heating applications.
Finally, experimental validation of these results is still missing. However, by using materials with low thermal expansion coefficients or keeping a small temperature range, the solutions to the proposed model are expected to be well behaved at high pressures when compared to experimental results. We have still to elaborate a model which is capable of taking into account thermal expansion effects on encapsulated PCMs. These effects may also have an impact on the thermal energy absorbed by the system. However, the EMB equation at the interface will be completely different when including these effects. Additionally, the role of mass conservation for obtaining the liquid and solid deformation during the melting process when the densities depend on the temperature profile is still not clear.

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

Abbreviations
The following abbreviations and symbols are used in this manuscript: