Next Article in Journal
Discovery of Cinnamylidene Derivative of Rhodanine with High Anthelmintic Activity against Rhabditis sp.
Next Article in Special Issue
Combination of Passive and Active Solar Heating with Thermal Energy Storage
Previous Article in Journal
Secondary Metabolites from Hericium erinaceus and Their Anti-Inflammatory Activities
Previous Article in Special Issue
Thermophysical Characterization of Paraffin Wax Based on Mass-Accommodation Methods Applied to a Cylindrical Thermal Energy-Storage Unit
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Effects of Volume Changes on the Thermal Performance of PCM Layers Subjected to Oscillations of the Ambient Temperature: Transient and Steady Periodic Regimes

by
Rubén D. Santiago-Acosta
1,
Ernesto M. Hernández-Cooper
1,
Rolando Pérez-Álvarez
2 and
José A. Otero
1,*
1
Tecnologico de Monterrey, Escuela de Ingeniería y Ciencias, Carr. al Lago de Guadalupe Km. 3.5, Atizapán de Zaragoza 52926, Mexico
2
Centro de Investigación en Ciencias, Universidad Autónoma del Estado de Morelos, Cuernavaca 52900, Mexico
*
Author to whom correspondence should be addressed.
Molecules 2022, 27(7), 2158; https://doi.org/10.3390/molecules27072158
Submission received: 24 February 2022 / Revised: 23 March 2022 / Accepted: 24 March 2022 / Published: 27 March 2022
(This article belongs to the Special Issue Phase Change Materials 2.0)

Abstract

:
The Stefan problem regarding the formation of several liquid–solid interfaces produced by the oscillations of the ambient temperature around the melting point of a phase change material has been addressed by several authors. Numerical and semi-analytical methods have been used to find the thermal response of a phase change material under these type of boundary conditions. However, volume changes produced by the moving fronts and their effects on the thermal performance of phase change materials have not been addressed. In this work, volume changes are incorporated through an additional equation of motion for the thickness of the system. The thickness of the phase change material becomes a dynamic variable of motion by imposing total mass conservation. The modified equation of motion for each interface is obtained by coupling mass conservation with a local energy–mass balance at each front. The dynamics of liquid–solid interface configurations is analyzed in the transient and steady periodic regimes. Finite element and heat balance integral methods are used to verify the consistency of the solutions to the proposed model. The heat balance integral method is modified and adapted to find approximate solutions when two fronts collide, and the temperature profiles are not smooth. Volumetric corrections to the sensible and latent heat released (absorbed) are introduced during front formation, annihilation, and in the presence of two fronts. Finally, the thermal energy released by the interior surface is estimated through the proposed model and compared with the solutions obtained through models proposed by other authors.

1. Introduction

Liquid–solid phase transitions have industrial applications, where on the one hand, phase change materials (PCMs) can be used as backup systems for thermoelectric generation [1,2] and on the other hand, as thermal barriers for air conditioning applications in the building industry [3,4]. The energy density of materials is enhanced through the latent heat absorbed in thermal energy storage units for thermoelectric generation during periods of low solar irradiance [5]. In contrast to applications related to thermoelectric generation, PCMs with low thermal conductivities are highly desirable when used as a thermal barrier to provide thermal comfort in building applications. The performance of PCMs used as thermal shields through PCM wallboards [4,6], PCM ceilings [7] and energy storage units for day/night cooling to reduce energy consumption [8], has been previously studied. Numerical simulations have been used to analyze the thermal performance of PCM walls, where the exterior surface is subjected to daily periodic boundary conditions and under real weather data [9]. The effects of natural convection on PCM layers with different orientations have been previously studied [10]. The authors conclude that energy transfer rates are reduced in PCM ceilings when the exterior surface is above the melting temperature T m and the liquid phase is in the conductive regime.
Numerical methods have been used to incorporate volume displacements produced by the density change during the phase transition [11,12,13,14]. The authors of Refs. [11,12] use several ways to accommodate mass during freezing of liquid water in confined systems. The authors perform numerical simulations with models that conserve mass properly, while considering the density difference between ice and liquid water. They found that solidification times are 19% higher and the energy absorbed by the PCM is 9% higher, when compared with models where the total mass of the system is not conserved. Finally, the authors assume incompressible phases, so the methods proposed can only be applied for isobaric transitions or materials with low isothermal compressibility. Thermo-mechanical models that consider the volume changes of a melting salt when confined in an elastic spherical shell have been developed [15]. The effects of the salt volume expansion played a significant role on the melting dynamics, thermal energy stored, and mechanical response of the shell. Experimental estimations of the latent heat, sensible heat, and mechanical energy stored in a micro-encapsulated KNO 3 / NaNO 3 salt were performed [16]. The thermo-mechanical model proposed in Ref. [15] was tested against the experimental estimations of the melting dynamics [16]. The results show a qualitative agreement between the the proposed model and the experimental estimations. The observed difference can be attributed to density changes produced by the thermal stress during the melting process and that are not considered in the original model [15]. Experimental studies of micro-encapsulated phase change materials (MEPCMs) with voids have also been performed [17]. The voids prevent the loss of latent heat storage capacity produced by the thermal stress during the melting process. Density variations produced by the thermal stress during melting of confined salts have also been studied. The model proposes a mass balance that considers the compressibility of each phase [18], extending the applicability of the original models [11,12,15] to the melting of confined PCMs.
Semi-analytical solutions to the Stefan problem on finite systems and with different types of boundary conditions have been reported in Refs. [19,20]. Periodic boundary conditions were preferred, since PCM layers used for thermal shielding applications are subjected to temperature oscillations produced by daily thermal variations [20]. The thermal performance of PCM walls and ceilings for building applications with periodic boundary conditions has also been addressed. Semi-analytical solutions to the Stefan problem in a finite PCM layer subjected to periodic boundary conditions have been previously reported [20]. The problem is simplified by assuming temperature oscillations on the external surface that are always above the melting temperature of the PCM. Additionally, temperature variations on the interior surface and below the melting temperature of the PCM were assumed. The situation produces the formation of one liquid-solid interface in the PCM layer that oscillates with the same frequency as the thermal oscillations on the layer boundaries [20]. The problem was also addressed in Ref. [21], where the authors use an enthalpy method and neglect volume changes produced by the oscillations of the interface. Experimental and numerical evidence in PCM layers subjected to temperature oscillations above the melting temperature on the external surface has been reported [22]. The temperature on the interior surface was fixed to a constant value and below the melting temperature of the PCM through a cooling system. Good agreement between the experimental results and a model where volume changes were not considered, was observed in the liquid phase. Discrepancies between the numerical predictions and the experimental results were observed on the temperature oscillations in the solid phase. According to the authors, the observed discrepancies could be produced by the assumption of one dimensional heat transfer. The weather or experimental conditions may produce the formation of several liquid–solid interfaces or fronts. Temperature oscillations around the melting temperature of the PCM layer on the external surface and homogeneous temperatures below the melting point of the PCM on the internal surface, have also been considered [23]. The authors found that during the heating process, the solid PCM melts on the exterior surface and the interface or front moves towards the interior surface. Additionally, when the ambient temperature is below T m , a thin solid layer in contact with the external surface is formed. The PCM layer is then divided into three regions and in the presence of two fronts during the cooling process. Formation of three liquid-solid interfaces may also occur, dividing the PCM layer into four regions. The configuration of phases and number of existing fronts depends on the amplitude of the temperature oscillations, the thickness of the PCM layer and thermodynamic properties of the PCM [23]. Periodical and non-sinusoidal boundary conditions have also been considered to address the effects of weather variations due to season changes during an entire year [24,25]. The phase configuration in the PCM layer when the exterior surface is subjected to variable periodic boundary conditions may change during each season. The Stefan problem with periodic boundary conditions around or above the melting temperature of the PCM, has been extensively studied. The problem studied by the mentioned authors has a one-dimensional treatment and no volume changes are incorporated. The literature that is focused on the Stefan problem with applications on thermal shielding, frequently neglects volume changes when describing the phase transition dynamics. It is commonly stated that volume changes in liquid-solid transitions represent at most 10% of the total volume of the system and can be neglected in the modeling by assuming phases with equal densities [26]. Recently, the issue of volume changes during isobaric transitions and in several types of organic PCMs was addressed to test this claim [26]. The authors found that, depending on the chemical composition of the PCM, volume changes can represent as much as 24% of the original volume in polymers. According to the authors, the treatment of volume changes experienced by the PCM during the phase transition represents a challenge that needs to be addressed in practical applications.
The effects of volume change on the thermal performance of PCM layers subjected to temperature oscillations about the melting temperature of the PCM has not been addressed. In this work, we consider the effects of volume changes on the energy released (absorbed) and energy transferred by the PCM in the presence two fronts. The energy transferred by the PCM is significantly influenced by the volume changes of the system and mainly in the presence of the two-front configuration, during the cooling stage of the cycle. A model that incorporates volume changes by promoting the layer thickness to a dynamic variable, is proposed. The oscillations of the layer thickness in the presence of several fronts obey an additional equation of motion that results from imposing total mass as a constant of the motion. The equation of motion for each interface is obtained through a local energy–mass balance that must be consistent with conservation of mass. The density used in each equation of motion results from the mass balance proposed in this work. The thermal performance of the PCM layer in the transient and steady periodic regimes, is analyzed through the thermal energy released (absorbed). The volumetric effects on the latent heat and sensible heat released (absorbed) during the solidification (melting) process are discussed in detail. Numerical and semi-analytical methods consisting on the finite element method (FEM) and heat balance integral method (HBIM), are used to solve the proposed equations of motion. Two different methods are used to verify the consistency of the solutions for the dynamical variables, released (absorbed) energy by the PCM layer and the energy released by the interior surface. The collision of two fronts produces a continuous temperature profile that is not smooth at the collision site. The HBIM is adapted by introducing a local energy balance at the collision site and predicting the time evolution of the temperature field just after the annihilation of two interfaces. The FEM is used to verify the consistency of the semi-analytical solutions. The total mass of the PCM layer is registered at every time interval to guarantee that no mass is created or destroyed during the entire phase change process. The time evolution of the released (absorbed) sensible and latent heats, and the energy released by the interior surface into the room are compared with models proposed by other authors. Finally, we find that the efficiency of a PCM layer as a thermal barrier close to the steady periodic regime may be reduced when considering the volume changes during the phase transition.

2. Description of the Physical System and Mathematical Model

The volumetric effects produced by the density difference between liquid and solid phases is considered when the external surface is subjected to periodic boundary conditions and the temperature of the internal surface is constant and below T m . The temperature dependence of liquid and solid densities is negligible within the temperature range and PCM considered in this work [27]. Octadecane is used as the PCM, and within the experimental error, the density of the solid is practically constant in a wide temperature range [ 260 , 301.13 ] K . The largest variation of the liquid density is less than 1.5% within the maximum temperature range in the liquid phase [ 301.13 , 313.15 ] K , that will be used in this work. Thermal expansion effects are negligible and only volume variations produced by the density difference between liquid and solid phases will be considered. Two case scenarios will be discussed:
  • Temperature oscillations on the external surface above the melting temperature of the PCM: one-front dynamics, and
  • Temperature oscillations around the melting temperature of the PCM: two-front dynamics with three phase coexistence, one-front dynamics with two phase coexistence and no phase change presence.
The system under consideration consists of a PCM that constitutes a layer of thickness L and cross section S with a left (outer) boundary in contact with the ambient air and the right (inner) boundary in contact with the air at the interior of a room. The thermal flux is perpendicular to the surface of the PCM layer and the temperature is uniformly distributed throughout the exterior and interior surfaces. Isothermal boundary conditions are employed at each liquid–solid interface with a temperature value equal to the melting temperature T m of the PCM. The right boundary at x = L ( t ) is subjected to isothermal boundary conditions and the left boundary at the outside air–PCM interface ( x = 0 ) is subjected to periodic boundary conditions, as follows:
T 2 0 , t = T 0 + δ sin ω t + ϕ and , T 2 ξ 1 ( t ) , t = T 1 ξ 1 ( t ) , t = T m , T 1 L ( t ) , t = T C ,
where T 0 represents the average daily temperature, δ the amplitude of the temperature oscillations, ω = 2 π / T the angular frequency with a 24 h period and ϕ is the phase angle. The liquid solid front at any time t is located at x = ξ 1 ( t ) . The oscillations of the ambient temperature were obtained through fitting weather data to a periodic function, in a tropical region located at Villahermosa, Tabasco in Mexico [28]. The temperature on the right boundary at x = L ( t ) is represented by T C and is below the melting temperature of the PCM. The temperature is uniform and homogeneously distributed throughout the interior surface. Maximum temperature differences in the liquid are T H T m T m , and in the solid phase are T m T C T m . The net heat flux at the interface is small in these temperature ranges and the phase transition takes place close to thermodynamic equilibrium; then, supercooling effects are not considered. Additionally, octadecane has a low supercooling degree within the operating temperature range considered in this work [29].

2.1. One-Front Dynamics: Transient and Steady Periodic Regime

The thickness of the PCM layer must incorporate volume changes during the phase transition according to total mass conservation and it is promoted to a dynamic variable defined as L ( t ) . Volume changes have been previously considered during phase transitions at constant pressure [14]. We briefly describe the corresponding equations of motion, which can be applied to any type of boundary conditions. In a PCM layer of thickness L ( t ) , when the temperature at x = 0 is above T m , the domain of the liquid phase lies on the interval: 0 x ξ 1 ( t ) , where ξ 1 ( t ) is the position of the liquid-solid interface. Additionally, the domain of the solid phase lies on the interval ξ 1 ( t ) x L ( t ) , where the boundary at x = L ( t ) is always below T m . The local energy-mass balance at x = ξ 1 ( t ) that is consistent with volume displacements on these type of configurations [14], is given by:
ρ L f d ξ 1 ( t ) d t = k T 2 ( ) ( x , t ) x | x = ξ 1 ( t ) + k s T 1 ( s ) ( x , t ) x | x = ξ 1 ( t ) ,
where T 1 ( s ) ( x , t ) ( T 2 ( ) ( x , t ) ) is the temperature distribution in the solid(liquid) phase as illustrated in Figure 1, ρ is the density of the liquid phase and L f is the latent heat of fusion of the PCM. The equation of motion for L ( t ) can be obtained straightforward through the total time derivative of the PCM’s mass per unit area A, and promoting this quantity as a constant of the motion, as follows:
ρ s d ξ 1 ( t ) d t + ρ s d L ( t ) d t d ξ 1 ( t ) d t = 0 .
Equations (2) and (3) are valid for melting and solidification processes, and consider volume changes during the phase transition produced by the density difference between liquid and solid. The density of the liquid phase appears in Equation (2) when the moving boundary is L ( t ) , as shown in Ref. [14]. The local energy–mass balance within the liquid and solid phases is taken into consideration through the following heat equation:
ρ i C i T i ( x , t ) t k i 2 T i ( x , t ) x 2 = 0 ,
where ρ i , C i and k i is the density, specific heat capacity and thermal conductivity of phase i, with i = ( s ) for liquid (solid). The one-front dynamics problem is defined through Equations (2)–(4), with the corresponding boundary conditions, given by Equation (1).
Homogeneous isothermal boundary conditions have been applied to the one-front dynamics problem where, on the one hand, the temperature at x = 0 is constant and higher than T m . On the other hand, the temperature at x = L ( t ) is constant and below T m . Steady state solutions for the thickness L ( t ) and liquid-solid interface position ξ 1 ( t ) have been found when the system is subjected to this type of boundary conditions [14]. Periodic boundary conditions are applied at x = 0 , which emulate the ambient temperature on the exterior surface of the PCM layer. The steady state values found in Ref. [14] can be applied to define upper and lower bounds to the interface position and layer thickness when the system is subjected to the boundary conditions shown in Equation (1). According to the steady state values for the liquid–solid interface position found in Ref. [14], an upper and lower bound for ξ 1 when the system is in the steady periodic regime and subjected to the boundary conditions given by Equation (1) are given by:
ξ s p ( u ) = k T 0 + δ T m ρ s L ( 0 ) ρ s ρ ξ 1 ( 0 ) ρ k T 0 + δ T m + ρ s k s T m T C , and
ξ s p ( l ) = k T 0 δ T m ρ s L ( 0 ) ρ s ρ ξ 1 ( 0 ) ρ k T 0 δ T m + ρ s k s T m T C ,
where ξ s p ( u ) and ξ s p ( l ) , represent the upper and lower bound of the interface position in the steady periodic regime. Highest and lowest temperature values are represented by T 0 + δ and T 0 δ . The initial interface position and thickness of the PCM layer are represented by ξ 1 ( 0 ) and L ( 0 ) , respectively. Additionally, we can apply the analytical expression for the thickness of the PCM layer in the steady state [14], to determine the upper and lower bounds of L in the steady periodic regime as follows:
L s p ( u ) = k s T m T C + k T 0 + δ T m ρ s L ( 0 ) ρ s ρ ξ 1 ( 0 ) ρ k T 0 + δ T m + ρ s k s T m T C , and
L s p ( l ) = k s T m T C + k T 0 δ T m ρ s L ( 0 ) ρ s ρ ξ 1 ( 0 ) ρ k T 0 δ T m + ρ s k s T m T C .
Here, L s p ( u ) ( L s p ( l ) ) represents the upper(lower) bound of L in the steady periodic regime. Mass conservation in Equation (2) was applied by promoting the boundary in contact with the interior at x = L ( t ) to a dynamical variable. Consequently, Equations (5)–(8) are only valid when the exterior surface of the PCM layer at x = 0 is fixed in time and the right boundary in contact with the interior at x = L ( t ) is the moving boundary.

2.2. Two-Front Dynamics

The formation of two or more liquid-solid fronts and coexistence of several adjacent liquid and solid phases, depends on the oscillations of the ambient temperature, the thermodynamic properties of the PCM and layer thickness. In this part of the section, the equations of motion that describe the effects of volume displacements on the dynamics of two liquid–solid interfaces will be introduced. The thermodynamic properties of the PCM, layer thickness and ambient temperature used in this work, will produce two liquid-solid fronts during the periods of the day when the ambient temperature is lower than T m . The two fronts will collide at some instant during the cooling stage, and one liquid-solid front will be formed and evolve in time when the ambient temperature oscillates above T m .
Initially, the configuration of the PCM layer will consist of two thin solid slabs and the ambient temperature is below T m . The ambient temperature lies below the melting temperature of the PCM for the first few hours of the initial part of the cycle, when the two initial solid fronts propagate and collide during the cooling stage. The PCM layer will be in its solid state after the collision when only sensible heat is absorbed, while the ambient temperature is increasing but is still below T m . The heating process starts when the ambient temperature is increasing and reaches T m . A thin liquid slab is formed on the exterior surface and a single front dynamics is observed during the time interval where the ambient temperature is above T m . The single front position is bounded by the value shown in Equation (5); therefore, a one-front configuration will always be present during the heating process of the cycle. The cooling process starts when the ambient temperature is decreasing and reaches T m . Finally, a thin solid slab is formed on the exterior surface while the existing front within the PCM layer is moving towards the recently formed interface. The PCM layer is under the presence of a two-front configuration during the cooling stage, and the process is repeated.
Figure 2, is an schematic illustration of the front and phase configuration during the periods of the day where two liquid-solid fronts are present. The temperature distribution within each phase is labeled according to regions: 1, 2 and 3 from right to left, as shown in Figure 2. Regions 1, 2 and 3 represent a solid-liquid-solid phase configuration with a temperature distribution T 1 ( s ) ( x , t ) , T 2 ( ) ( x , t ) and T 3 ( s ) ( x , t ) , respectively. The temperature on the interior surface (right boundary) at x = L ( t ) is constant and equal to T C < T m . Figure 2 assumes that the temperature on the exterior surface is below the melting temperature of the PCM; therefore, the liquid layer between ξ 2 ( t ) and ξ 1 ( t ) will gradually transform into solid. The mass of liquid Δ m ( 1 ) in close contact with ξ 1 ( t ) that is transformed into solid phase between t and t + Δ t can be obtained through a mass balance of the solid phase in region 1 as follows:
Δ m ( 1 ) = ρ s L ( t + Δ t ) x i 1 ( t + Δ t ) ρ s L ( t ) x i 1 ( t ) .
Here, the first term on the right hand side represents the mass of solid in region 1 at t + Δ t and the second term is the mass of solid at time t. Mass balance implies that Δ m ( 1 ) is equivalent to the mass of solid in region 1 that is formed between t and t + Δ t . The last equation can be simplified as follows:
Δ m ( 1 ) = ρ s Δ L ( t ) Δ ξ 1 ( t ) ,
where Δ L ( t ) = L ( t + Δ t ) L ( t ) and Δ ξ 1 ( t ) = x i 1 ( t + Δ t ) x i 1 ( t ) represent the displacement of the layer thickness and interface motion between t and t + Δ t , respectively. The solidification rate of liquid mass Δ m ( 1 ) in contact with the interface at x = ξ 1 ( t ) is produced when Δ m ( 1 ) releases thermal energy in form of latent heat during a small time interval Δ t . Consequently the thermal flux d Q s ( 1 ) / d t released at the interface position ξ 1 ( t ) is higher than the thermal flux absorbed from the liquid layer d Q ( 2 ) / d t at x = ξ 1 ( t ) , as shown in Figure 2. The energy–mass balance at x = ξ 1 ( t ) is then, given by
ρ s L f d L ( t ) d t d ξ 1 ( t ) d t = d Q s ( 1 ) d t d Q ( 2 ) d t .
The temperature profile in region 1 must be a decreasing function of x and must have a negative concavity at region 2 within the liquid layer, as illustrated in Figure 2. Consequently, the net thermal flux at x = ξ 1 ( t ) can be obtained as follows
d Q s ( 1 ) d t d Q ( 2 ) d t = k s T 1 ( s ) ( x , t ) x | x = ξ 1 ( t ) + k T 2 ( ) ( x , t ) x | x = ξ 1 ( t ) .
The rate of energy released by the liquid mass Δ m ( 1 ) is equal to the net thermal flux at ξ 1 ( t ) , and the energy–mass balance at ξ 1 ( t ) is given by
ρ s L f d L ( t ) d t d ξ 1 ( t ) d t = k s T 1 ( s ) ( x , t ) x | x = ξ 1 ( t ) + k T 2 ( ) ( x , t ) x | x = ξ 1 ( t )
Additionally, some of the liquid mass Δ m ( 2 ) adjacent to the solid phase at region 3, will transform into solid phase between t and t + Δ t . The mass of solid at region 3 will increase during this time interval, since the temperature on the exterior surface is below T m as illustrated in Figure 2. According to mass balance, Δ m ( 2 ) should be equal to the mass of solid that forms during this time interval, and is given by
Δ m ( 2 ) = ρ s ξ 2 ( t + Δ t ) ξ 2 ( t ) .
The energy released by this mass of liquid as latent heat, results from the net thermal flux at x = ξ 2 ( t ) . The mass of liquid Δ m ( 2 ) releases thermal energy to the solid phase at region 3 and absorbs thermal energy from the liquid layer. The rate of energy released d Q s ( 3 ) / d t at ξ 2 ( t ) , exceeds the rate of energy absorbed from the liquid layer d Q ( 2 ) / d t close to the interface ξ 2 ( t ) when solidification takes place. Additionally, the temperature profile at region 3 is an increasing function of x since the temperature on the exterior surface is below the melting temperature of the PCM, as shown in Figure 2. The temperature distribution in the liquid layer close to x = ξ 2 ( t ) shows a positive slope due to its concavity as illustrated in Figure 2. Consequently, the net amount of thermal flux at x = ξ 2 ( t ) is given by:
d Q s ( 3 ) d t d Q ( 2 ) d t = k s T 3 ( s ) ( x , t ) x | x = ξ 2 ( t ) k T 2 ( ) ( x , t ) x | x = ξ 2 ( t ) .
The rate of thermal energy transfer at x = ξ 2 ( t ) is equal to the rate of latent heat energy released by Δ m ( 2 ) . Applying a thermal balance between the last two equations, the local energy–mass balance at ξ 2 ( t ) can be obtained as follows
ρ s L f d ξ 2 ( t ) d t = k s T 3 ( s ) ( x , t ) x | x = ξ 2 ( t ) k T 2 ( ) ( x , t ) x | x = ξ 2 ( t ) .
The density that appears in Equations (13) and (16) must be the density of the phase that incorporates the correct mass balance. Conservation of the PCM’s total mass may now be applied to obtain the equation of motion for the thickness of the PCM layer. In the presence of two fronts, as shown in Figure 2, the total mass of the PCM is given by
m ( t ) = ρ s ξ 2 ( t ) + ρ ξ 1 ( t ) ξ 2 ( t ) + ρ s L ( t ) ξ 1 ( t ) .
The time derivative of m ( t ) is equal to zero when the mass of the PCM layer is conserved, and the equation of motion for L ( t ) is given by
ρ s d ξ 2 ( t ) d t + ρ d ξ 1 ( t ) d t d ξ 2 ( t ) d t + ρ s d L ( t ) d t d ξ 1 ( t ) d t = 0 .
Equations (13), (16), and (18) represent the equations of motion for ξ 1 ( t ) , ξ 2 ( t ) and L ( t ) that incorporate the volume changes induced by the presence of two moving fronts during the cooling stage of the cycle (when the exterior surface is below T m ).
The liquid mass Δ m ( 1 ) used to obtain Equation (13), can also be estimated by subtracting the mass of liquid Δ m ( 2 ) close to ξ 2 ( t ) from the total mass of liquid Δ m that will transform into solid between t and t + Δ t . The mass of liquid that will change into its solid state is given by:
m ( t ) m ( t + Δ t ) = ρ ξ 1 ( t ) ξ 2 ( t ) ρ ξ 1 ( t + Δ t ) ξ 2 ( t + Δ t ) ;
therefore, Δ m ( 1 ) can be obtained from Equation (14) and the last equation, as follows
Δ m ( 1 ) = ρ ξ 1 ( t ) ξ 2 ( t ) ρ ξ 1 ( t + Δ t ) ξ 2 ( t + Δ t ) ρ s ξ 2 ( t + Δ t ) ξ 2 ( t ) .
The rate of transformed mass Δ m ( 1 ) / Δ t when Δ t 0 is then:
d m ( 1 ) d t = ρ d ξ 2 ( t ) d t d ξ 1 ( t ) d t ρ s d ξ 2 ( t ) d t ,
which also results by solving ρ s d L ( t ) / d t d ξ 1 ( t ) / d t from Equation (18). Consequently, through a mass balance at the liquid layer or from total mass conservation, an equivalent equation of motion for ξ 1 ( t ) is obtained as follows
ρ L f 1 ρ s / ρ d ξ 2 ( t ) d t d ξ 1 ( t ) d t = k s T 1 ( s ) ( x , t ) x | x = ξ 1 ( t ) + k T 2 ( ) ( x , t ) x | x = ξ 1 ( t ) .
Equation (22) is completely equivalent to Equation (13) due to the mass balance at the liquid layer and total mass conservation. The two-front dynamics problem is described through Equations (13), (16), and (18) or equivalently through Equations (16), (18), and (22). The net thermal flux that is shown on the right hand side of Equations (13), (16), and (22) can be obtained through the temperature field at each phase, which is found by solving the local energy balance shown in Equation (4).

2.3. Volume Adjustments on Front Formation and Annihilation

The problem of several front formation takes place when the temperature on the exterior surface oscillates about the melting temperature of the PCM. Three possible scenarios will take place due to the thermodynamic properties, thickness of the PCM layer and the ambient temperature oscillations considered in this work. Solid phase will form when the ambient temperature reaches T m , decreasing towards its daily minimum, a liquid phase will form when the ambient temperature reaches T m and increases to its daily maximum, and the two fronts, ξ 1 ( t ) and ξ 2 ( t ) will collide during the time intervals when the ambient temperature is below the melting temperature of the PCM. The thickness of the PCM layer used in this work is such that ξ 1 ( t ) and ξ 2 ( t ) will meet at some instant when the temperature on the exterior surface is still below T m .
In this work, volume adjustments will be introduced during the creation of one-front and during the annihilation of ξ 1 ( t ) and ξ 2 ( t ) . Volume displacements of the PCM layer will be incorporated so that mass is not created or destroyed during the creation of a new phase. In this work, the supercooling of liquid is not considered, and it is assumed that the liquid–solid saturation temperature is equal to its value at thermodynamic equilibrium T m . A new solid phase will be formed and will be thermodynamically stable when the rate of energy released to the exterior d Q s ( 3 ) / d t by the newly formed solid phase of unknown thickness ξ 2 ( t a ) is equal to the energy released by the liquid phase d Q ( 2 ) / d t in contact with the new phase, as shown in Figure 3. The time instant t a when the new solid phase is stable, corresponds to any time value when the temperature on the exterior surface is equal or just below the melting temperature of the PCM. The thermal flux within the new solid phase and close to x = ξ 2 ( t a ) is small enough and equal to the thermal flux in the liquid layer close to ξ 2 ( t a ) , so that the new solid phase is thermodynamically stable at temperature values close to T m . Then, the thickness of the new phase can be found as follows
k s T 3 ( s ) ( x , t a ) x | x = ξ 2 ( t a ) = k T 2 ( ) ( x , t a ) x | x = ξ 2 ( t a ) .
The left hand side of the last equation represents the rate of energy released to the environment. We assume a linear temperature distribution within the new solid phase, and the last equation is reduced to the following:
k s T m T 0 + δ sin ω t a + ϕ ξ 2 ( t a ) = k T 2 ( ) ( x , t a ) x | x = ξ 2 ( t a ) ,
where the temperature on the exterior surface is T 0 + δ sin ω t a + ϕ as shown in Figure 3b. The volume of the PCM layer changes due to the formation of this phase, and the thickness of the PCM layer changes to avoid mass creation during the formation of the solid phase. The temperature at x = 0 , when t = t a Δ t is just above T m , and the total mass of the PCM layer is given by
m ( t a Δ t ) = ρ ξ 1 ( t a Δ t ) + ρ s L ( t a Δ t ) ξ 1 ( t a Δ t ) .
The total mass of the system m ( t a Δ t ) at t = t a Δ t should be equal to m ( t a ) ; therefore, imposing total mass conservation during the creation of the new solid phase, an equation for the thickness of the PCM layer is obtained as follows:
ρ ξ 1 ( t a Δ t ) + ρ s L ( t a Δ t ) ξ 1 ( t a Δ t ) = ρ s ξ 2 ( t a ) + ρ ξ 1 ( t a ) ξ 2 ( t a ) + ρ s L ( t a ) ξ 1 ( t a )
where the interface at x = ξ 1 ( t a ) = ξ ( t a Δ t ) L ( t a Δ t ) L ( t a ) is shifted to the left due to the volume shrinkage of the system during the formation of the solid phase as illustrated in Figure 3b. Substituting x = ξ 1 ( t a ) = ξ ( t a Δ t ) L ( t a Δ t ) L ( t a ) in the last equation, total mass conservation of the system between t a Δ t and t a is reduced to the following
ρ s ξ 2 ( t a ) 1 ρ / ρ s ρ L ( t a Δ t ) L ( t a ) = 0 .
Finally, Equations (24) and (27) can be solved for the thickness of the new phase ξ 2 ( t a ) and PCM layer L ( t a ) .
Additionally, a thin layer of liquid phase will form when the ambient temperature reaches T m , and increases towards the daily maximum temperature T 0 + δ . During the formation of the liquid layer, volume displacements are considered to avoid loss of total mass. The process is illustrated in Figure 4a,b, where at some time t = t a Δ t , the temperature on the exterior surface is just below T m and a thin layer of solid will transform into liquid phase when the temperature shifts to some value just above T m at some time t a . We do not consider overheating of solid phase during the transformation, and assume that a new liquid phase of unknown thickness ξ 1 ( t a ) is at thermodynamic equilibrium with the adjacent solid close to the new interface. The rate of energy absorbed by the newly formed liquid layer d Q 2 / d t is equal to the rate of energy released to the solid phase d Q s ( 1 ) / d t at x = ξ 1 ( t a ) , as follows
k T 2 ( ) ( x , t a ) x | x = ξ 1 ( t a ) = k s T 1 ( s ) ( x , t a ) x | x = ξ 1 ( t a ) .
The temperature field within the new liquid phase is close to the melting temperature of the PCM, and a linear profile is assumed; then, last equation is simplified as follows
k T 0 + δ sin ω t a + ϕ T m ξ 1 ( t a ) = k s T 1 ( s ) ( x , t a ) x | x = ξ 1 ( t a ) .
Approximating the temperature distribution in the solid phase (region 1) to a linear function at t = t a , as illustrated in Figure 4b, the last equation can be reduced as follows:
k T 0 + δ sin ω t a + ϕ T m ξ 1 ( t a ) = k s T m T C L ( t a ) ξ 1 ( t a ) ,
where T C is the temperature at the interior surface x = L ( t a ) . The corresponding expansion of the system produced by the formation of the new liquid phase can be obtained through total mass conservation. The expansion of the PCM layer is then, given by
ρ s L ( t a Δ t ) = ρ ξ 1 ( t a ) + ρ s L ( t a ) ξ 1 ( t a ) .
Solving the last two equations for the thickness of the new liquid phase ξ 1 ( t a ) and the thickness of the PCM layer L ( t a ) , the following approximate expressions are obtained:
ξ 1 ( t a ) = ρ s k T amb ( t a ) T m ρ k T amb ( t a ) T m + ρ s k s T m T C L ( t a Δ t ) and
L ( t a ) = ρ s k T amb ( t a ) T m + ρ s k s T m T C ρ k T amb ( t a ) T m + ρ s k s T m T C L ( t a Δ t ) ,
where T amb ( t a ) = T 0 + δ sin ω t a + ϕ is the temperature on the exterior surface x = 0 as shown in Figure 4b.
Finally, the collision between ξ 1 and ξ 2 will take place at some time t a during the daily lows, when the temperature on the exterior surface is below T m as illustrated in Figure 5. The two fronts will meet at some time t = t a , when a thin layer of liquid with thickness Δ ξ ( t a ) L ( t a ) remains as saturated liquid. The temperature distribution in the liquid layer of thickness Δ ξ ( t a ) = ξ 1 ( t a ) ξ 2 ( t a ) is practically equal to T m during the phase transition. The saturated liquid is assumed to be at thermodynamic equilibrium when the phase change starts. The transformation takes place at some time t a , when the thermal energy released by the liquid layer is exactly equal to the amount of latent heat that must be released to transform the mass of the remaining liquid. The total amount of thermal energy released by the liquid layer is given by
ρ L f Δ ξ ( t a ) = k s T 3 ( s ) ( x , t a ) x | x = ξ 2 ( t a ) k s T 1 ( s ) ( x , t a ) x | x = ξ 1 ( t a ) .
The thickness of liquid layer Δ ξ ( t a ) , when the phase change takes place can be found from the last equation as follows
Δ ξ ( t a ) = k s ρ L f T 3 ( s ) ( x , t a ) x | x = ξ 2 ( t a ) T 1 ( s ) ( x , t a ) x | x = ξ 1 ( t a )
The thickness of the PCM layer must decrease to avoid mass creation during the transformation to its solid state. Mass conservation can be used to find the thickness of the system once the liquid layer is transformed into its solid form as follows:
L ( t a ) = L ( t a ) Δ ξ ( t a ) 1 ρ / ρ s ,
where L ( t a ) represents the thickness of the system, once the liquid layer changes to its solid state.

3. Thermal Energy Released (Absorbed): Transient and Steady Periodic Regimes

The energy released (absorbed) by the PCM layer, when the exterior surface is subjected to temperature oscillations about T m will be discussed. First, we describe the energy released during the formation of a solid front ξ 2 and in the presence of two liquid–solid interfaces. Later, the energy released during the collision of two fronts and the energy released (absorbed) by the PCM layer in its solid state will be described. Finally, the energy absorbed during the formation of a liquid phase and the thermal energy absorbed (released) in the presence of a one interface will be discussed.
In this work, the sensible and latent heat released (absorbed) by the PCM layer are obtained independently. The thermal energy released (absorbed) by the system is estimated as the sum of the sensible and latent heat. The sensible heat released (absorbed) corresponds to internal energy differences in the PCM, produced by temperature changes from an initial state at t = t a to a final state at t = t b . Sensible heat estimations are performed through the integral of the entire temperature profiles in the PCM layer and provides complete information on the temperature distribution within the system. The FEM will be used to solve the proposed model and will be compared with the solutions estimated with the HBIM. The sensible heat released (absorbed) will be used as an indirect comparison between the temperature fields according to each method, and verify the consistency of the numerical and semi-analytical solutions. Additionally, it is observed that most of the contributions to the thermal oscillations come from the latent heat released (absorbed) by the PCM layer. The authors of Ref. [25] estimate the thermal energy released(absorbed) by the PCM layer through the time integral of the thermal flux that enters and exits the system. The net thermal flux through the layer only provides information on the behavior of the spatial derivative of the temperature at the exterior and interior surface.

3.1. Thermal Energy Released: Two-Front Configuration

The three phase configuration with two liquid–solid interfaces is present when the exterior temperature is below T m . A thin solid slab will form when the ambient temperature reaches T m and evolves towards the daily lowest temperature values. The system releases thermal energy as latent heat during the formation of a thin solid layer of thickness ξ 2 ( t a ) . The formation of the solid layer takes place at some time t a when the ambient temperature is just below T m (supercooling of liquid is not considered). The thickness of the thin solid layer ξ 2 ( t a ) can be found by solving Equations (24) and (27). Consequently, the latent heat released during the formation of this solid layer is given by
Δ Q f = ρ s L f ξ 2 ( t a )
Figure 6, shows the two-front configuration of the PCM layer at two different time instants t a and t b . The sensible heat released between t a and t b can be conceived through differences between the internal energy of the solid, the internal energy of the liquid mass that will be transformed into solid phase and the internal energy changes of the liquid that will not transform during this time interval Δ t = t b t a .
In Figure 6, Δ ξ 1 ( s ) = L ( t a ) ξ 1 ( t a ) represents the thickness of the solid phase in region 1 at time t = t a . After the phase transition, a fraction of the mass in the liquid layer between ξ 2 ( t a ) and ξ 1 ( t a ) will transform to its solid phase. The interior surface will shift to the left a distance equal to Δ L = L ( t a ) L ( t b ) as shown in Figure 6. The shift represents the volume displacement produced by the transformation of liquid into solid phase at constant pressure. The original solid in region 1, will be shifted to the left a distance ξ 1 ( t a ) ξ 1 ( t b ) as illustrated in Figure 6b. The shift can be found by applying mass conservation to the mass of solid in region 1, and solve for ξ 1 ( t b ) as follows
ξ 1 ( t b ) = ξ 1 ( t a ) + L ( t b ) L ( t a ) .
The internal energy change of the solid phase Δ U 1 between t a and t b can now be obtained, and is given by:
Δ U 1 = ρ s C s 0 ξ 2 ( t a ) T 3 ( s ) ( x , t b ) T 3 ( s ) ( x , t a ) d x + ρ s C s ξ 1 ( t b ) L ( t b ) T 1 ( s ) ( x , t b ) d x ξ 1 ( t a ) L ( t a ) T 1 ( s ) ( x , t a ) d x ,
where the value of ξ 1 ( t b ) is given by Equation (38).
The internal energy change of the liquid mass that is not transformed into solid between t a and t b is given by
Δ U 2 = ρ C ξ 2 ( t b ) ξ 1 ( t b ) T 2 ( ) ( x , t b ) T 2 ( ) ( x , t a ) d x
The total thickness in the liquid layer that will be transformed into solid during the time interval Δ t = t b t a is shown in Figure 6a. The temperature of the liquid mass that belongs to the regions in direct contact with ξ 1 ( t a ) and ξ 2 ( t a ) shown in Figure 6a will change from an initial value at t = t a to the melting temperature T m , just before the phase transition at some time t = t p t , between t a and t b . The total thickness of this fraction of liquid is: Δ ξ = Δ ξ 2 ( ) + Δ ξ 1 ( ) , where Δ ξ 2 ( ) = ξ 2 ( t p t ) ξ 2 ( t a ) and Δ ξ 1 ( ) = ξ 1 ( t a ) ξ 1 ( t t p ) . Mass conservation can be applied to the fraction of liquid close to ξ 2 ( t a ) as follows:
ρ ξ 2 ( t p t ) ξ 2 ( t a ) = ρ s ξ 2 ( t b ) ξ 2 ( t a ) ,
where the right hand side represents the mass of this fraction of liquid in its solid state. The expression for ξ 2 ( t p t ) can be found in terms of the known variables ξ 2 ( t a ) and ξ 2 ( t b ) , as follows
ξ 2 ( t p t ) = ξ 2 ( t a ) + ρ s ρ ξ 2 ( t b ) ξ 2 ( t a ) .
Additionally, mass conservation can be applied to the fraction of liquid close to ξ 1 ( t a ) in the following manner:
ρ ξ 1 ( t a ) ξ 1 ( t p t ) = ρ s ξ 1 ( t b ) ξ 1 ( t b ) ,
where ξ 1 ( t b ) is given by Equation (38). Consequently, according to the last equation ξ 1 ( t p t ) is given by
ξ 1 ( t p t ) = ξ 1 ( t a ) ρ s ρ ξ 1 ( t a ) ξ 1 ( t b ) + L ( t b ) L ( t a ) .
The internal energy change of the liquid mass that will transform to its solid phase, from an initial state with temperature T 2 ( ) ( x , t a ) to a state at the melting temperature of the PCM (saturated liquid), is given by:
Δ U 3 = ρ C ξ 2 ( t p t ) ξ 2 ( t a ) T m ξ 2 ( t a ) ξ 2 ( t p t ) T 2 ( ) ( x , t a ) d x + ρ C ξ 1 ( t a ) ξ 1 ( t p t ) T m ξ 1 ( t p t ) ξ 1 ( t a ) T 2 ( ) ( x , t a ) d x ,
where ξ 2 ( t p t ) and ξ 1 ( t p t ) are given by Equations (42) and (44), respectively. Finally, when this fraction of liquid is transformed into solid, it will release sensible heat from an initial state at t = t t p as saturated solid to its final state with a temperature distribution T 3 ( s ) ( x , t b ) between x = ξ 2 ( t a ) and x = ξ 2 ( t b ) , and T 1 ( s ) ( x , t b ) between x = ξ 1 ( t b ) and ξ 1 ( t b ) . The internal energy change of this fraction of transformed liquid and now in its solid state, is given by:
Δ U 4 = ρ s C s ξ 1 ( t b ) ξ 1 ( t b ) T 1 ( s ) ( x , t b ) d x ξ 1 ( t b ) ξ 1 ( t b ) T m + ρ s C s ξ 2 ( t a ) ξ 2 ( t b ) T 3 ( s ) ( x , t b ) d x ξ 2 ( t b ) ξ 2 ( t a ) T m ,
where ξ 1 ( t b ) is given by Equation (38).
The latent heat released between t a and t b can be obtained by estimating the fraction of liquid that is transformed into solid. The liquid mass Δ m that experiences the phase transition may be obtained by subtracting the liquid mass at t = t b from the liquid mass at t = t a , as follows:
Δ m = ρ ξ 1 ( t a ) ξ 2 ( t a ) ρ ξ 1 ( t b ) ξ 2 ( t b ) ;
therefore, the latent heat released during this time interval is given by
Δ Q f = ρ L f ξ 1 ( t b ) ξ 1 ( t a ) + ξ 2 ( t a ) ξ 2 ( t b ) .
According to the initial and final positions of ξ 1 and ξ 2 shown in Figure 6a,b, the latent heat Δ Q f given by Equation (48) will be negative.

3.2. Thermal Energy Released (Absorbed): Single Solid Phase

The collision of the two fronts ξ 1 and ξ 2 takes place at some time t = t a , as described in the previous section. The collision, implies the release of thermal energy through latent heat, that results from the transformation of the remaining liquid layer of thickness Δ ξ = ξ 1 ( t a ) ξ 2 ( t a ) , into its solid state. The latent heat released during the transformation of a thin liquid layer of thickness Δ ξ is given by
Δ Q f = ρ L f ξ 2 ( t a ) ξ 1 ( t a ) .
Here, Δ Q f is also negative, since liquid is transformed into solid phase. The thickness of the remaining liquid layer is: Δ ξ ( t a ) L ( t a ) , and we assume that the phase transition to its solid state is almost instantaneous. The thickness of this layer in its solid form can be obtained through mass conservation as follows:
ρ ξ 1 ( t a ) ξ 2 ( t a ) = ρ s ξ 1 ( t a ) ξ 2 ( t a ) ,
where ξ 1 ( t a ) is shown in Figure 5b and represents the position of the shifted interface after solidification of the liquid layer. The last equation can be solved to obtain an expression for ξ 1 ( t a ) , which is given by
ξ 1 ( t a ) = ξ 2 ( t a ) + ρ ρ s ξ 1 ( t a ) ξ 2 ( t a ) .
The solid at region 1 and illustrated in Figure 5a is shifted by an amount Δ L = ξ 1 ( t a ) ξ 1 ( t a ) as shown in Figure 5b, and the thickness L ( t a ) of the PCM layer after the phase change is given by Equation (36). The PCM layer only releases/absorbs thermal energy as sensible heat between t = t a and t = t b during the cooling stage of the cycle. Then, the internal energy change of the solid phase is given by:
Δ U s = ρ s C s 0 ξ 2 ( t a ) T 3 ( s ) ( x , t b ) T 3 ( s ) ( x , t a ) d x + ξ 2 ( t a ) ξ 1 ( t a ) T 2 ( s ) ( x , t b ) T 2 ( s ) ( x , t a ) d x + ρ s C s ξ 1 ( t a ) L ( t a ) T 1 ( s ) ( x , t b ) T 1 ( s ) ( x , t a ) d x ,
where ξ 1 ( t a ) is given by Equation (51) and shown in Figure 5b. The thickness of the PCM layer in its solid state L ( t a ) and the positions x = ξ 1 ( t a ) and x = ξ 2 ( t a ) only represent coordinates or reference positions when the PCM is releasing/absorbing sensible heat during the cooling stage. The values of ξ 1 ( t a ) and ξ 2 ( t a ) are constant in time and represent the location of the interfaces at the time of the collision and after the transition of the liquid layer, as shown in Figure 5.

3.3. Thermal Energy Absorbed (Released): One-Front Configuration

The two phase configuration in the presence of one liquid-solid interface is observed when the exterior temperature is above T m . A thin liquid slab will form when the ambient temperature reaches the melting temperature T m of the PCM, and evolves towards the maximum temperature values of the cycle. The system absorbs latent heat during the formation of a thin liquid layer of thickness ξ 1 ( t a ) , given by Equation (32). The latent heat absorbed during the formation of the liquid layer can be obtained as follows:
Δ Q f = ρ L f ξ 1 ( t a ) ,
where t a is the time value when the liquid layer is already formed. In the presence of one moving front, the PCM will absorb (release) energy through sensible and latent heat. The PCM layer will experience melting (solidification) with approximately the same frequency as the temperature oscillations.

3.3.1. Thermal Energy Absorbed: Melting

The PCM layer will absorb thermal energy from the environment, when the liquid-solid interface ξ 1 ( t ) moves towards the interior surface, as illustrated in Figure 7a,b. The sensible heat absorbed by the initial liquid mass at t = t a , during the time interval Δ t = t b t a is given by:
Δ U 1 ( m ) = ρ C 0 ξ 1 ( t a ) T 2 ( ) ( x , t b ) T 2 ( ) ( x , t a ) d x ,
where the super index ( m ) in Δ U ( m ) is used to specify internal energy changes during the melting stage, when ξ 1 ( t b ) > ξ 1 ( t a ) . Next, we must estimate the internal energy change of the solid layer, that will not be melted during the time interval Δ t = t b t a . In order to obtain Δ U 2 ( m ) , mass conservation is used to determine the thickness of solid Δ ξ 1 ( t a ) = ξ 1 ( t a ) ξ 1 ( t a ) at time t a and shown in Figure 7a, that will change to its liquid state. Applying mass conservation, the thickness of this fraction of solid can be obtained as follows:
ρ s ξ 1 ( t a ) ξ 1 ( t a ) = ρ ξ 1 ( t b ) ξ 1 ( t a ) ,
where ξ 1 ( t a ) is shown in Figure 7a, and solving for ξ 1 ( t a ) from the last equation; then,
ξ 1 ( t a ) = ξ 1 ( t a ) + ρ ρ s ξ 1 ( t b ) ξ 1 ( t a ) .
The internal energy change of the unmelted mass of solid is given by:
Δ U 2 ( m ) = ρ s C s ξ 1 ( t b ) L ( t b ) T 1 ( s ) ( x , t b ) d x ξ 1 ( t a ) L ( t a ) T 1 ( s ) ( x , t a ) d x ,
where ξ 1 ( t a ) is given by Equation (56).
The mass of solid that will melt during this time interval, will absorb sensible heat from an initial state at t a , where the temperature distribution is T 1 ( s ) ( x , t a ) to a final state at some time t p t between t a and t b , when the solid is at the melting temperature (saturated solid). The sensible heat absorbed by the fraction of melted solid between t a and t p t is given by
Δ U 3 ( m ) = ρ s C s ξ 1 ( t a ) ξ 1 ( t a ) T m ρ s C s ξ 1 ( t a ) ξ 1 ( t a ) T 1 ( s ) ( x , t a ) d x .
Finally, this mass of solid will absorb thermal energy from an initial state after the phase transition at some time between t p t as saturated liquid, to a final state at t b with a temperature distribution T 2 ( ) ( x , t b ) . Then, the sensible heat absorbed by this fraction of melted solid in its liquid phase, is given by
Δ U 4 ( m ) = ρ C ξ 1 ( t a ) ξ 1 ( t b ) T 2 ( ) ( x , t b ) d x ρ C ξ 1 ( t b ) ξ 1 ( t a ) T m .
The latent heat absorbed can be easily found by using the mass of melted solid (newly formed liquid), as follows
Δ Q f = ρ L f ξ 1 ( t b ) ξ 1 ( t a ) .
Equations (54) and (57)–(60), can only be applied during the melting stage, and in the presence of one interface.

3.3.2. Thermal Energy Released: Solidification

The solidification process presents a different scenario since the volume of the PCM layer is reduced during this part of the cycle. First, we consider the internal energy change experienced by the layer of liquid phase that will not transform into its solid phase between an initial state at t a and a final state at t b . The thickness of this fraction of liquid is equal to ξ 1 ( t b ) as illustrated in Figure 7c. The internal energy change experienced by this mass of liquid is given by:
Δ U 1 ( s ) = ρ C 0 ξ 1 ( t b ) T 2 ( ) ( x , t b ) T 2 ( ) ( x , t a ) d x ,
where the superindex ( s ) represents the solidification stage and Δ U 1 ( s ) is the internal energy change experienced by this fraction of liquid phase. The thickness L ( t a ) ξ 1 ( t a ) of the solid phase at t = t a and illustrated in Figure 7c, remains constant during the phase change process. This fraction of solid is shifted to the left a distance equal to Δ L = L ( t a ) L ( t b ) , due to the shrinkage of the liquid layer in contact with ξ 1 ( t a ) , as shown in Figure 7c,d. The contraction of the PCM layer Δ L can be used to determine the thickness of the liquid layer after the phase change as illustrated in Figure 7d, as follows:
ξ 1 ( t a ) ξ 1 ( t b ) = L ( t a ) L ( t b ) .
Consequently, ξ 1 ( t b ) is given by
ξ 1 ( t b ) = ξ 1 ( t a ) + L ( t b ) L ( t a ) .
The internal energy change experienced by the initial mass of solid can be obtained as follows:
Δ U 2 ( s ) = ρ s C s ξ 1 ( t b ) L ( t b ) T 1 ( s ) ( x , t b ) d x ξ 1 ( t a ) L ( t a ) T 1 ( s ) ( x , t a ) d x ,
where ξ 1 ( t b ) is shown in Figure 7d and given by Equation (63). The fraction of liquid that will be transformed to its solid phase, will release sensible heat by changing from an initial state at t = t a to a saturated liquid state at some time t = t p t between t a and t b . The thickness of this fraction of liquid is ξ 1 ( t a ) ξ 1 ( t b ) and is illustrated in Figure 7c. The sensible heat released during this process can be obtained as follows
Δ U 3 ( s ) = ρ C ξ 1 ( t a ) ξ 1 ( t b ) T m ρ C ξ 1 ( t b ) ξ 1 ( t a ) T 2 ( ) ( x , t a ) d x .
Additionally, when this mass of liquid changes to its solid phase, it will release sensible heat from an initial state as saturated solid at t = t p t , to a final state at t = t b , where its temperature distribution is T 1 ( s ) ( x , t b ) . The sensible heat released by the newly formed solid phase is given by:
Δ U 4 ( s ) = ρ s C s ξ 1 ( t b ) ξ 1 ( t b ) T 1 ( s ) ( x , t b ) d x ρ s C s ξ 1 ( t b ) ξ 1 ( t b ) T m ,
where ξ 1 ( t b ) is given by Equation (63).
Finally, the latent heat released during the phase transition is given by
Δ Q f = ρ L f ξ 1 ( t b ) ξ 1 ( t a )

4. Numerical and Semi-Analytical Methods

Front tracking methods are applied to solve the model described in Section 2 with the volume corrections proposed in this work. The HBIM is used to find approximate semi-analytical solutions, and a FEM with first order Lagrange interpolating functions is used to verify the consistency of the semi-analytical solutions. The HBIM demands continuity and smoothness of the temperature profile in each phase. The isothermal boundary condition at the liquid-solid front, introduces a discontinuity in the spacial derivative of the temperature when two fronts collide and the entire PCM layer is in its solid phase. In this work, the HBIM is modified by introducing a local energy balance at the collision site and the solutions are verified through comparison with the FEM.

4.1. Heat Balance Integral Method

The HBIM is adapted to find approximate analytical solutions to the model described in Section 2. Continuous and smooth temperature profiles at each phase are required to apply the HBIM used in Refs. [14,30,31]. The HBIM consists on proposing a polynomial function for the temperature at each phase. The space integral of Equation (4) is performed at each phase, to obtain a set of ordinary differential equations (ODE) in time. Quadratic functions of x that satisfy the boundary conditions, are proposed in this work. The following temperature profiles are proposed to solve the two-front configuration scenario illustrated in Figure 2, as follows:
T 1 ( s ) ( x , t ) = a 1 ( t ) x ξ 1 ( t ) + b 1 ( t ) x ξ 1 ( t ) 2 + T m , for ξ 1 ( t ) x L ( t ) , T 2 ( ) ( x , t ) = a 2 ( t ) x ξ 1 ( t ) + b 2 ( t ) x ξ 1 ( t ) 2 + T m , for ξ 2 ( t ) x ξ 1 ( t ) , T 3 ( s ) ( x , t ) = a 3 ( t ) x ξ 2 ( t ) + b 3 ( t ) x ξ 2 ( t ) 2 + T m , for 0 x ξ 2 ( t ) .
Here, the time dependent coefficients a i ( t ) at each region i = 1 , 2 , 3 can be expressed in terms of b i ( t ) through the boundary conditions given by:
T 1 ( s ) L ( t ) , t = T C , T 2 ( ) ξ 2 ( t ) , t = T m , and T 3 ( s ) 0 , t = T 0 + δ sin ω t + ϕ .
The boundary condition at x = ξ 1 ( t ) is satisfied by the proposed temperature profiles. Applying the boundary conditions shown in the last equation the following time dependent coefficients are found:
a 1 ( t ) = b 1 ( t ) L ( t ) ξ 1 ( t ) + T m T C ξ 1 ( t ) L ( t ) , a 2 ( t ) = b 2 ( t ) ξ 1 ( t ) ξ 2 ( t ) , and a 3 ( t ) = b 3 ( t ) ξ 2 ( t ) + T m T 0 δ sin ω t + ϕ ξ 2 ( t )
The temperature profiles are substituted into Equation (4) and the result is integrated over the domain of each phase to obtain a set of three ODEs in time for the coefficients b i ( t ) . The two-front configuration of the PCM layer becomes a dynamic problem for the time dependent variables ξ 1 ( t ) , ξ 2 ( t ) , L ( t ) and b i ( t ) with i = 1 , 2 , 3 . The three ODEs in time for the coefficients b i ( t ) , along with Equations (13), (16) and (18), constitute a set of six ODEs that is solved through an explicit finite difference method with a forward first order approximation to the time derivatives.
The one-front configuration problem illustrated in Figure 1 is solved similarly. The space integral of Equation (4) is performed at each region shown in Figure 1. Two ODEs in time for the coefficients b i ( t ) with i = 1 , 2 is obtained through the spacial integrals. The PCM layer in the presence of one-front, becomes a dynamic problem for the time dependent variables ξ 1 ( t ) , L ( t ) and b i ( t ) with i = 1 , 2 [14].
The HBIM just described, requires continuity and smoothness of each temperature profile within its domain. The method must be slightly adapted to scenarios where the temperature distribution is not smooth as illustrated in Figure 5b. The situation depicted in Figure 5b results from the collision of two fronts and just after the saturated liquid is transformed into its solid phase. The two-front collision scenario shown in Figure 5b takes place at some time t a , when the temperature is approaching the daily minimum, and the entire PCM layer is in its solid state. The isothermal boundary condition at ξ 1 and ξ 2 is no longer required. The liquid-solid interfaces disappear and the thickness of the PCM layer becomes a constant of the motion in the absence of a phase change process. The problem consists of solving Equation (4) with the initial temperature profile illustrated in Figure 5b. A single polynomial function with the initial temperature field shown in Figure 5b is not possible, and instead a piecewise function that is consistent with the initial conditions is proposed as follows:
T 1 ( s ) ( x , t ) = a 1 ( t ) x ξ 1 + b 1 ( t ) x ξ 1 2 + T u ( t ) , for ξ 1 x L , T 2 ( s ) ( x , t ) = T u ( t ) , for ξ 2 x ξ 1 , and T 3 ( s ) ( x , t ) = a 3 ( t ) x ξ 2 + b 3 ( t ) x ξ 2 2 + T u ( t ) , for 0 x ξ 2 ,
where T u ( t ) is the temperature at the collision site of thickness Δ ξ = ξ 1 ξ 2 and initially, equal to the melting temperature of the PCM. Initially, the region of thickness Δ ξ represents a saturated solid. The HBIM can be applied to regions 1 and 3 through the space integral of Equation (4), as previously described. The coordinates x = ξ 1 and x = ξ 2 now represent locations in space within the solid, and define a thin solid layer of thickness Δ ξ L . In this work, instead of solving the heat equation in region 2, we assume a uniform temperature profile T u ( t ) in this thin layer that will change in time according to a local energy balance principle. Following the basic idea of energy balance, the thin solid layer of thickness Δ ξ , will release(absorb) thermal energy that comes from the net heat flux at x = ξ 1 and x = ξ 2 , as illustrated in Figure 5b. The thin solid layer Δ ξ , releases thermal energy when the temperature is distributed through the PCM layer, as illustrated in Figure 5b. The rate of energy released by the solid layer Δ ξ during a small time interval Δ t , is equal to the internal energy change Δ U s experienced by this layer, as follows
Δ U s Δ t = C s ρ s T u ( t ) T u ( t + Δ t ) Δ t .
The internal energy change shown through the last equation, results from the energy released to the solid in regions 1 and 3 as depicted in Figure 5b. The net rate of thermal energy transferred can be obtained as follows
d Q s ( 1 ) d t + d Q s ( 3 ) d t = k s T 1 ( s ) ( x , t ) x | x = ξ 1 + k s T 3 ( s ) ( x , t ) x | x = ξ 2 .
According to the last two equations, the local energy balance at region 2 in the limit Δ t 0 , is given by
C s ρ s d T u ( t ) d t = k s T 1 ( s ) ( x , t ) x | x = ξ 1 k s T 3 ( s ) ( x , t ) x | x = ξ 2 ,
which constitutes a differential equation for the temperature T u ( t ) in the thin solid layer between x = ξ 2 and x = ξ 1 . Additionally, the energy balance given by Equation (74) may also be applied when the thin solid layer absorbs thermal energy. Finally, the time evolution of T u ( t ) is determined through the solution of Equation (74), which is used to estimate the coefficients b i ( t ) with i = 1 , 3 through the classical HBIM, previously described.

4.2. Finite Element Method

The FEM was applied to find the spacial dependence of the temperature field in each phase [32]. Equation (4) was solved through the methodology described in Ref. [32], but using the linear Lagrange shape functions at each element as follows:
N 1 ( x ) = x x 2 x 1 x 2 , and N 2 ( x ) = x x 1 x 2 x 1 .
Here x 1 and x 2 represent the nodal coordinates of any given element. The temperature field within each element is given by:
T ˜ ( x , t ) = N 1 ( x ) T ^ 1 ( t ) + N 2 ( x ) T ^ 2 ( t ) ,
where, T ^ 1 ( t ) and T ^ 2 ( t ) is the time dependent part of the temperature at each node. The time evolution of the temperature at each node was obtained through a first order approximation to the time derivative of Equation (4). The implicit finite difference scheme was applied to estimate the nodal temperatures at the next time level.
Equations (13), (16), and (18) were solved for the dynamical variables ξ 1 ( t ) , ξ 2 ( t ) , and L ( t ) by using an explicit finite difference scheme with a first order approximation to the time derivatives. The position of each interface and the thickness of the PCM layer in the next time level were used to solve the heat equation through the FEM previously mentioned [32].

5. Results and Discussion

The PCM used as an insulating material is octadecane. The thermal performance of a PCM layer of octadecane is determined by estimating the time evolution of the dynamical variables, the thermal energy released (absorbed) and the thermal energy released by the interior surface. The thermodynamic properties of octadecane are assumed to be constant in the temperature range considered in this work, and equal to their values close to the saturation temperature of the PCM [27]. For the liquid(solid) phase k = 0.152 ( k s = 0.334 ) W / m · K , C = 1.921 ( C s = 2.230 ) kJ / kg · K and ρ = 776.860 ( ρ s = 867.914 ) kg / m 3 . The liquid-solid saturation properties are: L f = 236.98 kJ / kg and T m = 301.13 K [27].

5.1. One-Front Dynamics: Transient and Steady Periodic Regimes

The one-front dynamics scenario is present when the ambient temperature oscillates above the melting temperature of Octadecane during the complete cycle. In this case, the parameters for the ambient ambient temperature are: T 0 = 308.15 K , δ = 5.0 K and ϕ = 0.73779 rad . The phase angle ϕ represents a shift of the sine function, and the temperature at the inner surface is fixed at T C = 295.15 K .
Equations (2)–(4) were solved through the FEM and the HBIM described in the previous section. Solutions were found for different values of the Stefan number defined as:
SteNo = L f C T max T m ,
where T max = T 0 + δ is the maximum temperature on the exterior surface in this example. Figure 8, shows the transient and steady periodic parts of the solution for the interface position ξ 1 ( t ) and the thickness of the PCM layer L ( t ) . Approximate HBIM solutions are validated through FEM solutions in the one-front dynamics problem as shown in Figure 8. We have predicted lower and upper bounds with volume corrections, for the interface position and thickness of the PCM layer, as shown through Equations (5)–(8). Figure 8 shows solutions in the steady periodic regime for several values of SteNo . According to Equations (5)–(8), the interface position and thickness of the PCM layer are bounded in the steady periodic regime. The solutions are observed to oscillate within the predicted bounds for several values of SteNo .
The maxima and minima in the oscillations of ξ 1 ( t ) and L ( t ) close to the steady periodic regime, are tested by probing the solutions in the range SteNo = [ 008841 88.4120 ] . The Stefan number is modified by changing the magnitude of the latent heat. According to Equations (5)–(8), upper and lower bounds for ξ 1 and L close to the steady periodic regime, should not depend on L f when the system is subjected to non-homogeneous isothermal boundary conditions. The steady periodic solutions must approach asymptotically to the lower and upper bounds predicted through Equations (5)–(8) for low values of L f . The asymptotic behavior close to the steady periodic regime according to Equations (5)–(8), was captured by the HBIM and the FEM solutions. Maxima and minima in the oscillations of ξ 1 ( t ) and L ( t ) , were registered for several values of the Stefan number as shown in Figure 9. Melting (solidification) rates are expected to increase for low values of the latent heat as observed in Figure 9; therefore, for small Stefan numbers, the amplitude of the oscillations in the steady periodic regime is exactly bounded by the values predicted through Equations (5)–(8).

5.2. Two-Front Dynamics: Transient and Steady Periodic Regimes

The temperature on the exterior surface was extracted from weather data of 10 August 2021 at the city of Villahermosa, Tabasco in Mexico [28]. The temperature data was fitted to the periodic function shown in Equation (1), with the following fitting parameters: T 0 = 302.884 K , δ = 5.24632 K and ϕ = 0.73779 rad . The melting temperature of the PCM selected for the numerical and semi-analytical examples, and the ambient temperature of the selected region can give rise to the formation of several fronts. The ambient temperature oscillates around the melting temperature of octadecane, which is chosen as the PCM for the numerical and semi-analytical examples. The initial thickness of the PCM layer is L 0 = 3.0 cm , and a two-front configuration is observed during part of the cooling stage of the cycle. The temperature at the inner surface is T C = 295.15 . The inner temperature of T C = 295.15 has been chosen to observe the two-front configuration scenario with the longest possible duration, but with inner temperatures that produce some thermal comfort. Lowest values of T C are possible; however, volumetric effects would be less evident.
The model described through Equations (4), (13), (16), and (18) will be solved with the FEM and HBIM discussed in the previous section. The model is equivalent to Equations (4), (16), (18), and (22), where Equation (22) results from applying total mass conservation to Equation (13). Figure 10 shows the FEM and HBIM solutions for ξ 1 ( t ) , ξ 2 ( t ) , and L ( t ) in a Octadecane layer with an initial thickness of L 0 = 3.0 cm . Initially, the system is almost in its liquid phase, with a thin solid layer of thickness ξ 2 ( 0 ) = 1.0 mm close to the exterior surface. Additionally, a solid layer of L ( 0 ) ξ 1 ( 0 ) = 1.0 mm of thickness lies close to the interior surface. The system is close to a steady periodic regime after being exposed to the ambient temperature for 4–5 days.
Figure 10 contains the information on the front configuration during a complete cycle: two-front, one-front and single solid phase configurations. The transient part of the solution is observed for t 1 day , where the layer thickness is significantly reduced. Initially, the ambient temperature is decreasing towards the daily lows and a significant amount of thermal energy is released by the PCM layer.
Internal energy and latent heat changes were estimated during the transient and steady periodic regimes. The model described through Equations (4), (13), (16), and (18) is solved, and thermal energy changes during small time intervals of Δ t = 0.1 s are estimated. Figure 11 shows the sensible heat, latent heat, and total energy released (absorbed) by the PCM during the transient and steady periodic regimes. The thermal energy released (absorbed) can also be found by performing the time integral of the net thermal flux through the PCM layer [25]. The net rate of thermal energy change in the PCM can be obtained through the difference between the thermal flux on the exterior and inner surface as follows
d Q PCM d t = d Q ext d t d Q in d t .
The energy released (absorbed) by the PCM layer between t = 0 and t = t a is then, given by:
Δ Q PCM = 0 t a k i T j ( i ) ( x , t ) x | x = 0 + k s T 1 ( s ) ( x , t ) x | x = L ( t ) d t ,
where k i represents the thermal conductivity of phase i = , s and T j ( i ) ( x , t ) is the temperature distribution of phase i within region j = 2 , 3 . The last equation may be used instead of estimating the thermal energy released (absorbed) through the process described in Section 3. On the one hand, Equation (79) contains global information related with the volumetric effects on the sensible and latent heat released (absorbed) by the PCM, and only provides information on the total energy released (absorbed). On the other hand, the process described in Section 3 includes detailed information related with the manner in which sensible and latent heats are being released or absorbed. Additionally, Equation (79) only depends on the behavior of the temperature on the exterior and inner surfaces, and does not provide information related with the entire temperature field in the PCM. The process described in Section 3 is also preferred, since the time evolution of the temperature profiles integrated over the PCM domain, can be used to perform a more consistent comparison between the numerical and semi-analytical solutions used in this work.
Thermal energy is released during the solidification stages, which are observed when the ambient temperature is below T m (night hours) and during the cooling part of the day. The amount of sensible heat released by the PCM in the presence of two liquid–solid interfaces is obtained by substituting the solutions at two different time values t a and t b in Equations (39), (40), (45), and (46). The sensible heat released when the system is in its pure solid state and during the intervals with lowest ambient temperatures is estimated through Equation (52). Additionally, during the solidification process, and in the presence of one-front, the sensible heat released is obtained through Equations (61) and (64)–(66). Latent heat is released during the solidification processes, during the formation of a solid layer and at the instant of collision between ξ 1 and ξ 2 . The latent heat released in the presence of two fronts and one-front, was estimated through Equation (48) and (67), respectively. Latent heat released during the formation of a solid layer close to the exterior surface and at the time of collision was obtained through Equation (37) and (49), respectively. Similarly, the sensible heat absorbed during the melting process was estimated through Equations (54) and (57)–(59). Sensible heat is also absorbed when the PCM layer is in its solid state and during the intervals when the ambient temperature is increasing from its daily minimum and towards the melting temperature of the PCM. Finally, the latent heat absorbed during the formation of a liquid layer and during the melting process, is obtained through Equations (53) and (60), respectively.
Table 1 shows the highest relative percent difference (RPD) between the model discussed in this work and constant volume models commonly by other authors. The RPD was obtained for the sensible heat, latent heat and thermal energy released (absorbed) through each method used in this work. The RPD was estimated as follows:
RPD Δ E ( max , i ) = | Δ E p Δ E o | | Δ E p | × 100 % ,
where Δ E p and Δ E o represents the sensible heat, latent heat or thermal energy released(absorbed) by the PCM layer according to the proposed model in this work and the constant volume methods used by other authors, respectively. The maximum RPD with i = 1 and i = 2 corresponds to the RPD between the proposed model ( ρ s ρ ) and a constant volume model with ρ s = ρ , and the RPD between the proposed model ( ρ s ρ ) and a constant volume model with ρ s = ρ , respectively. The maximum RPD is found when the system is close to the steady periodic regime and is shown in Table 1.
Numerical and semi-analytical solutions to models that do not incorporate volumetric effects are also shown in Figure 11. Two scenarios are considered by assuming an octadecane sample of L = 3.0 cm with phases of equal densities. On the one hand, the density of the liquid and solid phase are equal to ρ s = 867.914 kg / m 3 . On the other hand, the density of both phases is equal to the liquid density ρ = 776.860 kg / m 3 . The solutions for Δ U , Δ Q f , Δ Q , and total mass are labeled as ρ = ρ s and ρ s = ρ in Figure 11. Volume changes expansion or shrinkage disappear by assuming a PCM with no density change during the phase transition. In absence of volume changes, total mass is implicitly conserved and an additional equation of motion for the thickness L of the PCM layer is not required. Thickness and total mass are constants of the motion when equal densities are assumed. Sensible and latent heats can be obtained by taking the corresponding limit of equal densities on each of the equations discussed in Section 3.
The effects of PCM expansion and shrinkage can be observed in Figure 11 and depend on the initial conditions of the problem. Initially, the system is practically in its liquid form, since ξ 2 ( 0 ) = 1.0 mm and ξ 1 ( 0 ) = 2.9 cm . Additionally, the temperature on the exterior surface is below T m and evolving towards the daily minimum. The initial conditions chosen on this example produce solidification of the liquid layer with an initial thickness of Δ ξ liq ( 0 ) = ξ 1 ( 0 ) ξ 2 ( 0 ) = 2.8 cm . The shrinkage of the system from its initial state to its pure solid form can be obtained through a simple mass balance where ρ Δ ξ liq ( 0 ) = ρ s Δ ξ sol . Here, Δ ξ sol is the thickness of the initial mass of liquid, but in its solid form. Consequently, after solidification of the initial liquid layer, the system must shrink an amount equal to: Δ ξ liq ( 0 ) Δ ξ sol = 2.938 mm ; therefore, the thickness of the PCM layer should be L ( t a ) = 2.706 cm at the time t a of the first collision, as observed in Figure 10b. Assuming for example, a PCM where the liquid phase has the same density as the PCM in its solid form ρ = ρ s , the mass of the PCM layer is higher, as shown in Figure 11d. The liquid phase has a larger mass in this case, and consequently the latent heat released is significantly higher as illustrated in Figure 11b. The latent heat initially released, is shifted along the energy axis as a consequence of the extra mass that results by assuming ρ = ρ s . The other scenario, when ρ s = ρ , does not produce the observed shift on the latent heat, since total mass does not change significantly.
Finally, Figure 12 shows the energy released by the solid in region 1 and to the interior of the room. The energy released by the interior surface, represents the amount of thermal energy that must be removed from the room to keep a constant temperature at the inner surface of the PCM layer. The energy released by the solid phase at region 1 can be obtained through the time integral of the thermal flux on the inner surface as follows
Δ Q in = k s 0 t a T 1 ( s ) ( x , t ) x | x = L ( t ) d t .
The energy released by the inner surface represents another way to evaluate the thermal efficiency of the PCM layer as a thermal barrier. The result is compared with the models of other authors, where volumetric effects are not considered.
The enhanced thermal performance of the PCM when equal densities are assumed, can be understood in terms of the initial state of the system and the thickness of the solid phase in region 1. On the one hand, the thickness of the PCM layer remains constant in this case, while it is reduced initially by almost 3.0 mm when ρ ρ s , as previously discussed. Close to the steady periodic regime, the thickness of the PCM layer will oscillate between 2.7 cm and 2.8 cm as a consequence of this initial shrinkage. On the other hand, the thickness of the PCM layer will remain constant and equal to its initial value of L = 3.0 cm when equal densities are assumed. The oscillations of the interface at x = ξ 1 ( t ) are very similar in both cases (equal and different densities); therefore, the effective thickness of the solid phase in region 1 and close to the steady periodic regime, will be smaller when volumetric effects are incorporated into the problem. Consequently, since the temperature gradient in region 1 is Δ T = T m T C in both cases, higher energy transfer rates at x = L ( t ) will be expected when ρ ρ s .

6. Conclusions

The main goal of this work consisted in estimating the effects of volume changes in the thermal performance of a PCM layer, when the external surface is exposed to temperature oscillations about the fusion point of the PCM. Despite the relatively small density variations in this material, the volumetric effects on the energy transferred by the PCM layer were significant. The initial conditions of the system produced an overall maximum volume shrinkage of 10% during the transient regime, which is observed when the two fronts collide for the first time. The latent heat, sensible heat, and thermal energy released in the transient regime are shifted along the energy axis as a result of the initial shrinkage of the PCM layer. The shift of the thermal energy initially released produces a significant relative percent difference between the model proposed in this work and the constant volume methods used by other authors. The largest relative percent difference is observed when the liquid density is equal to the solid density. The result is expected since assuming phases with densities equal to ρ s , the total mass of the system is increased. Additionally, it is found that the system does not recover its initial thickness when the PCM layer oscillates in the steady periodic regime. The thickness of the PCM layer is significantly reduced and oscillates around a smaller value in the steady periodic regime, constituting a less effective thermal barrier as a result of the initial system shrinkage. The results indicate that the initial state of the system has a significant impact on the thermal performance of the PCM, which is overestimated when neglecting volume changes. The situation may be reversed, and thermal performance could be enhanced by using an initial state with two liquid fronts that produce a system expansion from which the PCM layer may not recover. The result will depend on the thermodynamic properties of the PCM and the temperature oscillations on the exterior surface. The volumetric effects produced by a combination of different transient states has yet to be determined. Different transient states may appear in systems with non-sinusoidal temperature oscillations. The system is unable to reach a steady periodic regime in this case, and the volumetric effects predicted in this work may be more profound.
Additional contributions to volume changes produced by considering the temperature dependence of liquid and solid densities must be addressed as well. Thermal expansion effects on the problem of several front formation will depend on the type of PCM. Ambient temperature variations may also give rise to a different phase configuration than discussed in this work when considering the thermal expansion of each phase. Finally, to the authors’ knowledge, there is little experimental evidence in the literature concerning the problem of several front formation. The present work is still limited by the experimental validations that will allow estimation of the effects predicted through the proposed model.

Author Contributions

Conceptualization and discussion of proposed model was done by R.D.S.-A., E.M.H.-C., R.P.-Á. and J.A.O., the HBIM was developed by R.D.S.-A. and E.M.H.-C., the FEM was implemented by J.A.O., the manuscript was written by R.D.S.-A. and E.M.H.-C. and revised by J.A.O. and R.P.-Á., animations and figures were elaborated by J.A.O. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations and symbols are used in this manuscript:
PCMPhase change material
FEMFinite element method
HBIMHeat balance integral method
SteNoStefan number
k Thermal conductivity of the liquid
k s Thermal conductivity of the solid
C Specific heat capacity of the liquid
C s Specific heat capacity of the solid
ρ Liquid density
ρ s Solid density
T m Melting temperature
L f Latent heat of fusion
LThickness of PCM layer
ξ s p ( u ) Upper bound for interface position in the steady periodic regime
ξ s p ( ) Lower bound for interface position in the steady periodic regime
L s p ( u ) Upper bound for the layer thickness in the steady periodic regime
L s p ( ) Lower bound the layer thickness in the steady periodic regime
ξ i Position of the ith liquid-solid interface
T j ( i ) ( x , t ) Temperature profile of phase i in region j
T j ( i ) ( 0 , t ) Temperature of the exterior surface
T amb ( t ) Ambient temperature
T C Temperature of the inner surface
T 0 Daily average temperature
δ Amplitude of temperature oscillations
ω Angular frequency of temperature oscillations
Δ Q Thermal energy that penetrates the PCM layer
L f Latent heat of fusion
T m Melting temperature
Δ U Internal energy change
Δ Q f released(absorbed) latent heat
Δ Q released(absorbed) thermal energy
Δ Q in Thermal energy released by the interior surface

References

  1. Gil, A.; Medrano, M.; Martorell, I.; Lázaro, A.; Dolado, P.; Zalba, B.; Cabeza, L.F. State of the art on high temperature thermal energy storage for power generation. Part 1-Concepts, materials and modellization. Renew. Sustain. Energy Rev. 2010, 14, 31–55. [Google Scholar] [CrossRef]
  2. Liu, M.; Saman, W.; Bruno, F. Review on storage materials and thermal performance enhancement techniques for high temperature phase change thermal storage systems. Renew. Sustain. Energy Rev. 2012, 16, 2118–2132. [Google Scholar] [CrossRef]
  3. Zhou, D.; Zhao, C.Y.; Tian, Y. Review on thermal energy storage with phase change materials (PCMs) in building applications. Appl. Energy 2012, 92, 593–605. [Google Scholar] [CrossRef] [Green Version]
  4. Scalat, S.; Banu, D.; Hawes, D.; Parish, J.; Haghighata, F.; Feldman, D. Full scale thermal testing of latent heat storage in wallboard. Sol. Energy Mater. Sol. Cells 1996, 44, 49–61. [Google Scholar] [CrossRef]
  5. Ibrahim, N.I.; Al-Sulaiman, F.A.; Rahman, S.; Yilbas, B.S.; Sahin, A.Z. Heat transfer enhancement of phase change materials for thermal energy storage applications: A critical review. Renew. Sustain. Energy Rev. 2017, 74, 26–50. [Google Scholar] [CrossRef]
  6. Cabeza, L.F.; Castellon, C.; Nogues, M.; Medrano, M.; Leppers, R.; Zubillaga, O. Use of microencapsulated PCM in concrete walls for energy savings. Energy Build. 2007, 39, 113–119. [Google Scholar] [CrossRef]
  7. Saman, W.Y.; Belusko, M. Roof Integrated Unglazed Transpired Solar Air Heater. Ph.D. Thesis, ANZSES, Hobart, Australia, 1997. [Google Scholar]
  8. Vakilaltojjar, S.M.; Saman, W. Analysis and modelling of a phase change storage system for air conditioning applications. Appl. Therm. Eng. 2001, 21, 249–263. [Google Scholar] [CrossRef]
  9. Mathieu-Potvin, F.; Gosselin, L. Thermal shielding of multilayer walls with phase change materials under different transient boundary conditions. Int. J. Therm. Sci. 2009, 48, 1707–1717. [Google Scholar] [CrossRef]
  10. Madruga, S.; Curbelo, J. Effect of the inclination angle on the transient melting dynamics and heat transfer of a phase change material. Phys. Fluids 2021, 33, 055110. [Google Scholar] [CrossRef]
  11. Dallaire, J.; Gosselin, L. Various ways to take into account density change in solid-liquid phase change models: Formulation and consequences. Int. J. Heat Mass Transf. 2016, 103, 672. [Google Scholar] [CrossRef]
  12. Dallaire, J.; Gosselin, L. Numerical modeling of solid-liquid phase change in a closed 2D cavity with density change, elastic wall and natural convection. Int. J. Heat Mass Transf. 2017, 114, 903. [Google Scholar] [CrossRef]
  13. Hetmaniok, E.; Slota, D.; Zielonka, A. Solution of the Direct Alloy Solidification Problem including the Phenomenon of Material Shrinkage. Therm. Sci. 2017, 21, 105. [Google Scholar] [CrossRef] [Green Version]
  14. Santiago, R.D.; Hernández, E.M.; Otero, J.A. Constant mass model for the liquid-solid phase transition on a one-dimensional Stefan problem: Transient and steady state regimes. Int. J. Therm. Sci. 2017, 118, 40. [Google Scholar] [CrossRef]
  15. Lopez, J.; Caceres, G.; Del Barrio, E.P.; Jomaa, W. Confined melting in deformable porous media: A first attempt to explain the graphite/salt composites behaviour. Int. J. Heat Mass Transf. 2010, 53, 1195–1207. [Google Scholar] [CrossRef]
  16. Pitié, F.; Zhao, C.; Cáceres, G. Thermo-mechanical analysis of ceramic encapsulated phase-change-material (PCM) particles. Energy Environ. Sci. 2011, 4, 2117–2124. [Google Scholar] [CrossRef]
  17. Nomura, T.; Zhu, C.; Sheng, N.; Saito, G.; Akiyama, T. Microencapsulation of metal-based phase change material for high-temperature thermal energy storage. Sci. Rep. 2015, 5, 1–8. [Google Scholar] [CrossRef] [Green Version]
  18. Hernández, E.M.; Otero, J.A. Fundamental incorporation of the density change during melting of a confined phase change material. J. Appl. Phys. 2018, 123, 085105. [Google Scholar] [CrossRef]
  19. Lamberg, P. Approximate analytical model for two-phase solidification problem in a finned phase-change material storage. Appl. Energy 2004, 77, 131–152. [Google Scholar] [CrossRef]
  20. Mazzeo, D.; Oliveti, G.; De Simone, M.; Arcuri, N. Analytical model for solidification and melting in a finite PCM in steady periodic regime. Int. J. Heat Mass Transf. 2015, 88, 844–861. [Google Scholar] [CrossRef]
  21. Ho, C.; Chu, C. Periodic melting within a square enclosure with an oscillatory surface temperature. Int. J. Heat Mass Transf. 1993, 36, 725–733. [Google Scholar] [CrossRef]
  22. Casano, G.; Piva, S. Experimental and numerical investigation of the steady periodic solid–liquid phase-change heat transfer. Int. J. Heat Mass Transf. 2002, 45, 4181–4190. [Google Scholar] [CrossRef]
  23. Chang-Yong, C.; Hsieh, C. Solution of Stefan problems imposed with cyclic temperature and flux boundary conditions. Int. J. Heat Mass Transf. 1992, 35, 1181–1195. [Google Scholar] [CrossRef]
  24. Mazzeo, D.; Oliveti, G. Thermal field and heat storage in a cyclic phase change process caused by several moving melting and solidification interfaces in the layer. Int. J. Therm. Sci. 2018, 129, 462–488. [Google Scholar] [CrossRef]
  25. Mazzeo, D.; Oliveti, G.; Arcuri, N. Dynamic Parameters to Characterize the Thermal Behaviour of a Layer Subject to Periodic Phase Changes. Energy Procedia 2016, 101, 129–136. [Google Scholar] [CrossRef]
  26. Cabeza, L.F.; Zsembinszki, G.; Martín, M. Evaluation of volume change in phase change materials during their phase transition. J. Energy Storage 2020, 28, 101206. [Google Scholar] [CrossRef]
  27. Faden, M.; Höhlein, S.; Wanner, J.; König-Haagen, A.; Brüggemann, D. Review of thermophysical property data of octadecane for phase-change studies. Materials 2019, 12, 2974. [Google Scholar] [CrossRef] [Green Version]
  28. Estaciones Meteorológicas Automáticas (EMAS). Available online: https://smn.conagua.gob.mx (accessed on 26 September 2021).
  29. Van Miltenburg, J. Fitting the heat capacity of liquid n-alkanes: New measurements of n-heptadecane and n-octadecane. Thermochim. Acta 2000, 343, 57–62. [Google Scholar] [CrossRef]
  30. Mitchell, S. An accurate nodal heat balance integral method with spatial subdivision. Numer. Heat Transf. Part B Fundam. 2011, 60, 34–56. [Google Scholar] [CrossRef]
  31. Mitchell, S.L.; Myers, T.G. Application of heat balance integral methods to one-dimensional phase change problems. Int. J. Differ. Equ. 2012, 2012, 187902. [Google Scholar] [CrossRef] [Green Version]
  32. Rodríguez-Alemán, S.; Hernández-Cooper, E.M.; Pérez-Álvarez, R.; Otero, J.A. Effects of Total Thermal Balance on the Thermal Energy Absorbed or Released by a High-Temperature Phase Change Material. Molecules 2021, 26, 365. [Google Scholar] [CrossRef]
Figure 1. Schematic representation of the PCM layer in the presence of one liquid–solid interface or front with the temperature profiles in the liquid phase (region 2) and solid phase (region 1).
Figure 1. Schematic representation of the PCM layer in the presence of one liquid–solid interface or front with the temperature profiles in the liquid phase (region 2) and solid phase (region 1).
Molecules 27 02158 g001
Figure 2. Three phase configuration with a PCM layer in the presence of two fronts. The two-front configuration is observed when the temperature on the exterior surface is below T m . The liquid layer releases thermal energy to the solid phases at regions 1 and 3 (supercooling effects are not considered).
Figure 2. Three phase configuration with a PCM layer in the presence of two fronts. The two-front configuration is observed when the temperature on the exterior surface is below T m . The liquid layer releases thermal energy to the solid phases at regions 1 and 3 (supercooling effects are not considered).
Molecules 27 02158 g002
Figure 3. Solid phase formation close to the exterior surface produced by the change of the ambient temperature from T amb ( t a Δ t ) > T m to T amb ( t a ) < T m during a small time interval Δ t . (a) Thin liquid layer that will change to its solid phase. The layer is still in its liquid form at t = t a Δ t , when the temperature on the exterior surface is just above the melting temperature of the PCM. (b) Volume shrinkage Δ L = L ( t a Δ t ) L ( t a ) produced by the shift of the interior surface, after the liquid layer in contact with the exterior surface is transformed into its solid phase.
Figure 3. Solid phase formation close to the exterior surface produced by the change of the ambient temperature from T amb ( t a Δ t ) > T m to T amb ( t a ) < T m during a small time interval Δ t . (a) Thin liquid layer that will change to its solid phase. The layer is still in its liquid form at t = t a Δ t , when the temperature on the exterior surface is just above the melting temperature of the PCM. (b) Volume shrinkage Δ L = L ( t a Δ t ) L ( t a ) produced by the shift of the interior surface, after the liquid layer in contact with the exterior surface is transformed into its solid phase.
Molecules 27 02158 g003
Figure 4. Liquid phase formation close to the exterior surface produced by the change of the ambient temperature from T amb ( t a Δ t ) < T m to T amb ( t a ) > T m during a small time interval Δ t . (a) Thin solid layer that will melt close to the exterior surface. The layer is still in its solid form at t = t a Δ t , when the temperature on the exterior surface is just below T m . (b) Volume expansion Δ L = L ( t a ) L ( t a Δ t ) after the formation of the liquid layer. The exterior surface is pushed rightwards, due to the expansion of the solid layer after the transition to its liquid state.
Figure 4. Liquid phase formation close to the exterior surface produced by the change of the ambient temperature from T amb ( t a Δ t ) < T m to T amb ( t a ) > T m during a small time interval Δ t . (a) Thin solid layer that will melt close to the exterior surface. The layer is still in its solid form at t = t a Δ t , when the temperature on the exterior surface is just below T m . (b) Volume expansion Δ L = L ( t a ) L ( t a Δ t ) after the formation of the liquid layer. The exterior surface is pushed rightwards, due to the expansion of the solid layer after the transition to its liquid state.
Molecules 27 02158 g004
Figure 5. Solid phase formation during the collision of ξ 1 ( t ) and ξ 2 ( t ) . (a) Saturated liquid layer before the phase transition at some time t = t a . (b) PCM layer contraction Δ L = L ( t a ) L ( t a ) after the phase transition. The exterior surface is pulled to the left due to volume changes produced by the difference between solid and liquid densities.
Figure 5. Solid phase formation during the collision of ξ 1 ( t ) and ξ 2 ( t ) . (a) Saturated liquid layer before the phase transition at some time t = t a . (b) PCM layer contraction Δ L = L ( t a ) L ( t a ) after the phase transition. The exterior surface is pulled to the left due to volume changes produced by the difference between solid and liquid densities.
Molecules 27 02158 g005
Figure 6. Schematic illustration of liquid layers ξ 2 ( t p t ) ξ 2 ( t a ) and ξ 1 ( t a ) ξ 1 ( t p t ) , before and after the phase transition. (a) Saturated liquid layers before the transformation at some time t = t a . The thickness of each liquid layer ξ 2 ( t p t ) ξ 2 ( t a ) and ξ 1 ( t a ) ξ 1 ( t p t ) is obtained by applying mass conservation. (b) Transformed liquid layers after the phase transition at some time t = t b . The thickness of each solid layer ξ 2 ( t b ) ξ 2 ( t a ) and ξ 1 ( t b ) ξ 1 ( t b ) is determined through mass conservation.
Figure 6. Schematic illustration of liquid layers ξ 2 ( t p t ) ξ 2 ( t a ) and ξ 1 ( t a ) ξ 1 ( t p t ) , before and after the phase transition. (a) Saturated liquid layers before the transformation at some time t = t a . The thickness of each liquid layer ξ 2 ( t p t ) ξ 2 ( t a ) and ξ 1 ( t a ) ξ 1 ( t p t ) is obtained by applying mass conservation. (b) Transformed liquid layers after the phase transition at some time t = t b . The thickness of each solid layer ξ 2 ( t b ) ξ 2 ( t a ) and ξ 1 ( t b ) ξ 1 ( t b ) is determined through mass conservation.
Molecules 27 02158 g006
Figure 7. Volume changes during melting and solidification in a PCM layer with a two phase configuration (one-front). (a) Saturated solid layer in direct contact with ξ 1 ( t a ) that will transform into its liquid phase. (b) Transformed solid layer during the melting process between t a and t b . (c) Liquid layer at some time t = t a that will change to its solid phase. (d) Transformed liquid layer in contact with ξ 1 ( t b ) after the phase transition.
Figure 7. Volume changes during melting and solidification in a PCM layer with a two phase configuration (one-front). (a) Saturated solid layer in direct contact with ξ 1 ( t a ) that will transform into its liquid phase. (b) Transformed solid layer during the melting process between t a and t b . (c) Liquid layer at some time t = t a that will change to its solid phase. (d) Transformed liquid layer in contact with ξ 1 ( t b ) after the phase transition.
Molecules 27 02158 g007
Figure 8. Transient and steady periodic behavior of ξ 1 ( t ) and L ( t ) , and for several values of the Stefan number. Solid and dashed lines represent HBIM and FEM solutions, respectively for: (a) ξ 1 ( t ) and (b) layer thickness L ( t ) . The amplitude of the oscillations in ξ 1 ( t ) and L ( t ) are tuned through the Stefan number by changing the value of L f . Corresponding values of the Stefan number are indicated through black arrows. Lower and upper bounds in the steady periodic regime given by Equations (5)–(8) are shown through horizontal dotted lines.
Figure 8. Transient and steady periodic behavior of ξ 1 ( t ) and L ( t ) , and for several values of the Stefan number. Solid and dashed lines represent HBIM and FEM solutions, respectively for: (a) ξ 1 ( t ) and (b) layer thickness L ( t ) . The amplitude of the oscillations in ξ 1 ( t ) and L ( t ) are tuned through the Stefan number by changing the value of L f . Corresponding values of the Stefan number are indicated through black arrows. Lower and upper bounds in the steady periodic regime given by Equations (5)–(8) are shown through horizontal dotted lines.
Molecules 27 02158 g008
Figure 9. (a) HBIM and FEM steady periodic solutions as a function of the natural logarithm of the inverse Stefan number ln SteNo 1 for ln ξ 1 ( max ) ( t ) / L 0 and ln ξ 1 ( min ) ( t ) / L 0 , and (b) ln L ( max ) ( t ) / L 0 and ln L ( min ) ( t ) / L 0 . Maxima and minima in ξ 1 ( t ) and L ( t ) were registered close to the steady periodic regime at t = 9 days.
Figure 9. (a) HBIM and FEM steady periodic solutions as a function of the natural logarithm of the inverse Stefan number ln SteNo 1 for ln ξ 1 ( max ) ( t ) / L 0 and ln ξ 1 ( min ) ( t ) / L 0 , and (b) ln L ( max ) ( t ) / L 0 and ln L ( min ) ( t ) / L 0 . Maxima and minima in ξ 1 ( t ) and L ( t ) were registered close to the steady periodic regime at t = 9 days.
Molecules 27 02158 g009
Figure 10. Time evolution of: (a) ξ 1 ( t ) , ξ 2 ( t ) , and (b) L ( t ) , according to the HBIM and FEM solutions to Equations (4), (13), (16), and (18). Solid and dotted lines represent the FEM and HBIM solutions, respectively.
Figure 10. Time evolution of: (a) ξ 1 ( t ) , ξ 2 ( t ) , and (b) L ( t ) , according to the HBIM and FEM solutions to Equations (4), (13), (16), and (18). Solid and dotted lines represent the FEM and HBIM solutions, respectively.
Molecules 27 02158 g010
Figure 11. Time evolution of (a) internal energy change Δ U , (b) latent heat released (absorbed) Δ Q f , and (c) total energy changes Δ Q per unit area in kWh / m 2 according to the model proposed in this work ( ρ s ρ ) and the constant volume method used by other authors ( ρ s = ρ or ρ = ρ s ) [21,23,24,25]. (d) Time evolution of total mass is registered to verify that mass is not created or destroyed during the transient and steady periodic regimes. The effects of volume changes on the thermal energy released (absorbed) during solidification, melting, front formation and collision are discussed in Section 3.
Figure 11. Time evolution of (a) internal energy change Δ U , (b) latent heat released (absorbed) Δ Q f , and (c) total energy changes Δ Q per unit area in kWh / m 2 according to the model proposed in this work ( ρ s ρ ) and the constant volume method used by other authors ( ρ s = ρ or ρ = ρ s ) [21,23,24,25]. (d) Time evolution of total mass is registered to verify that mass is not created or destroyed during the transient and steady periodic regimes. The effects of volume changes on the thermal energy released (absorbed) during solidification, melting, front formation and collision are discussed in Section 3.
Molecules 27 02158 g011
Figure 12. Energy transferred to the interior of the room according to the numerical and semi-analytical solutions used in this work. Constant volume methods of Refs. [21,23,24,25] are compared with the result obtained through the proposed model in this work. According to the initial state of the system and the result shown in this figure, the thermal performance of the PCM layer is enhanced by assuming equal densities.
Figure 12. Energy transferred to the interior of the room according to the numerical and semi-analytical solutions used in this work. Constant volume methods of Refs. [21,23,24,25] are compared with the result obtained through the proposed model in this work. According to the initial state of the system and the result shown in this figure, the thermal performance of the PCM layer is enhanced by assuming equal densities.
Molecules 27 02158 g012
Table 1. Maximum RPD for the sensible heat, latent heat and thermal energy released (absorbed) between the model discussed in this work and the constant volume methods used by other authors [21,23,24,25].
Table 1. Maximum RPD for the sensible heat, latent heat and thermal energy released (absorbed) between the model discussed in this work and the constant volume methods used by other authors [21,23,24,25].
RPD Δ U ( max , 1 ) RPD Δ Q f ( max , 1 ) RPD Δ Q ( max , 1 ) RPD Δ U ( max , 2 ) RPD Δ Q f ( max , 2 ) RPD Δ Q ( max , 2 )
FEM3.73%2.50%2.50%13.01%11.72%11.78%
HBIM4.01%2.51%2.52%12.84%11.72%11.77%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Santiago-Acosta, R.D.; Hernández-Cooper, E.M.; Pérez-Álvarez, R.; Otero, J.A. Effects of Volume Changes on the Thermal Performance of PCM Layers Subjected to Oscillations of the Ambient Temperature: Transient and Steady Periodic Regimes. Molecules 2022, 27, 2158. https://doi.org/10.3390/molecules27072158

AMA Style

Santiago-Acosta RD, Hernández-Cooper EM, Pérez-Álvarez R, Otero JA. Effects of Volume Changes on the Thermal Performance of PCM Layers Subjected to Oscillations of the Ambient Temperature: Transient and Steady Periodic Regimes. Molecules. 2022; 27(7):2158. https://doi.org/10.3390/molecules27072158

Chicago/Turabian Style

Santiago-Acosta, Rubén D., Ernesto M. Hernández-Cooper, Rolando Pérez-Álvarez, and José A. Otero. 2022. "Effects of Volume Changes on the Thermal Performance of PCM Layers Subjected to Oscillations of the Ambient Temperature: Transient and Steady Periodic Regimes" Molecules 27, no. 7: 2158. https://doi.org/10.3390/molecules27072158

APA Style

Santiago-Acosta, R. D., Hernández-Cooper, E. M., Pérez-Álvarez, R., & Otero, J. A. (2022). Effects of Volume Changes on the Thermal Performance of PCM Layers Subjected to Oscillations of the Ambient Temperature: Transient and Steady Periodic Regimes. Molecules, 27(7), 2158. https://doi.org/10.3390/molecules27072158

Article Metrics

Back to TopTop