Complex Fourier series expansion for the liquid-solid phase transition in PCM layers: transient and steady state periodic regimes

Semi-analytical solutions to the classical two phase Stefan problem are proposed. Time dependent solutions to the one-dimensional liquid-solid phase transition in a PCM wallboard subjected to isothermal and periodic Dirichlet boundary conditions are obtained. Transient and steady state solutions are found in finite size systems, and the semi-analytical solutions are validated through the asymptotic time limit behaviour of the phase transition. In this work, complex Fourier methods are proposed to find the solutions in the transient and steady state periodic regimes. Semi-analytical solutions based on the heat balance integral method (HBIM) are used to verify the consistency of the proposed method. The Fourier method can be pictured as a generalization of the phasors based method recently introduced by other authors. The proposed method incorporates a complete set of complex functions, which allows finding the transient and steady state response of the system. Finally, solutions for the time dependent interface position, liquid and solid temperature distributions and the thermal energy penetrating through the PCM wallboard, are shown. The solutions from the proposed method are found to be consistent when compared to the semi-analytical solutions estimated through the HBIM.


Introduction
The study of phase change processes through the analysis of moving boundary problems has a wide variety of applications. These applications rely on the use of phase change materials (PCMs) in concentrating solar power (CSP) plants for thermoelectric generation as described by Gil and Mathur [1][2]. Additionally, solar energy harvesting for domestic hot water systems has been studied by Cabeza and Xie [3][4]. High temperature phase change materials (HTPCMs) have been studied by Cabeza, Xie and Cá ceres [3][4][5]. HTPCMs are used as latent heat thermal energy storage (LHTES) units in backup systems to provide thermal energy when sunlight is not available. PCMs with mild fusion temperatures between 323 K and 333 K have been used by Cabeza [3] in hot water systems for domestic applications. Additionally, PCMs with lower melting temperatures may find applications for thermal shielding in the construction industry, as discussed by Hadjieva and Ye [6][7].
The usage of PCMs in thermal shielding applications are extremely appealing due to the isothermal nature of the phase transition and thermal comfort provided by the absorption(release) of thermal energy during the light(dark) period of the day. Current uses of PCMs are related to reducing electric energy consumption costs by usage of heating or cooling units in buildings and residences. Graphite composites with salt hydrates as PCMs have been used by Hadjieva [6] to determine the structural and thermo-physical properties of wallboards. Concrete composites with CaCl2 6H2O salt hydrate as the PCM, have been used by Ye [7], to estimate the thermal performance of a PCM based panel room with an artificial climatic chamber. Khan, et al. [8] have studied salt hydrates as PCMs for thermal energy storage and thermal shielding applications due to the high latent heat values and low to mild liquid-solid saturation temperatures at standard pressure. However, hydrate salts tend to experience segregation during the melting process according to Kahn [8]. Paraffins, on the other hand, do not experience segregation but have lower latent heat values, which may diminish their thermal performance in building applications according to a review given by Kahn, et al [8]. Buildings based on PCM wallboards in places with extreme day/night temperatures can provide thermal comfort by absorbing and releasing thermal energy during the day and night periods, reducing electric consumption costs associated with the air conditioning unit.
The nonlinear nature of the equation of motion for the liquid-solid interface has also been an appealing property associated with continuum models of first order phase change processes. Several authors are devoted on finding numerical solutions for the liquid-solid interface position, volume changes, amount of melted (solidified) mass through enthalpy and front tracking methods. Radiation and convection effects during the phase change process have been taken into account by Archibold [9] through enthalpy methods. Front tracking methods have also been used by Santiago [10][11] and Hetmaniok [12]. Santiago [10][11] used finite difference and heat balance integral methods to determine the volume changes experienced during an isobaric phase transition by using finite difference and heat balance integral methods. Density changes and pressure increments during melting of confined salts with high melting temperatures have also been addressed by Herná ndez [13] through front tracking methods. Incompressible solid phases were assumed by Lopez [14] and enthalpy-based methods were used to determine the fracture limit of spherically encapsulated salts. Incompressible liquid and solid phases were assumed by Dallaire [15][16] to apply mass accommodation methods during the freezing process of water in a deformable cavity.
Savovic [17] and Mazzeo [18][19] have used front tracking methods to determine the thermal performance of hydrate salts and paraffins during isobaric phase change processes in wallboards subjected to periodic boundary conditions. The work carried out by Savovic [17] proposes an explicit finite difference scheme with time dependent Dirichlet and Neumann boundary conditions. Recent works by Mazzeo [18][19] consider a two-layer system with daily periodic boundary conditions where the interface position and thermal flux are obtained for long time values. Mathematical techniques based on phasors were applied to find analytical solutions to the thermal diffusion equations at each phase and the equation of motion for the liquid-solid interface. The complexity of the problem relies on the non-linear nature of the resulting system of differential equations due to the motion of the liquid-solid interface.
In this work, a method based on a Fourier series to analyze the behavior of the interface, is established. Timedependent boundary conditions are considered by taking into account the periodic nature of the ambient temperature during a 24-hour time window. Periodic Dirichlet boundary conditions are applied at the surface of the wallboard in contact with the ambient and homogeneous Dirichlet (Isothermal) boundary conditions are applied at the inner surface of the wall. The result obtained by Mazzeo [18][19] is generalized through the proposed method, where analytical solutions to the transient part of the solution are also obtained. The first section is devoted to discuss the model and the set of base functions required to solve the thermal diffusion problem. The second section is dedicated to analyze the problem with periodic conditions at the wall boundaries in the long time limit. In the third section, the initial boundary conditions and periodic conditions at the boundaries of the system are discussed. Finally, the conclusions are presented in the last section of this work.

Model
Consider a wallboard of length and cross section A with a PCM consisting of its two phases (1: solid and 2: liquid). The time dependent liquid-solid interface position is denoted as . The liquid phase lies between the left boundary and the interface, and the solid phase lies between the interface and the right boundary. The temperature at the left and right boundaries is homogeneously distributed, so the problem becomes onedimensional. The liquid-solid phase transition is held at constant pressure and buoyancy effects produced by temperature dependent densities within the liquid are not considered. Then, heat transfer through buoyancy and pressure induced natural convection is negligible. The equations that describe the behavior of the temperature are given by: Where ( − 1) ≤ ≤ (2 − ) + ( − 1) and the heat diffusion coefficients in phase are = ⁄ , which depend on the specific heat capacities , the densities and the thermal conductivities of each phase. In addition, the energy balance at the interface is given by: where is the latent heat of fusion. The expected solution is the piecewise function of the temperature given by: where the temperature at the interface position is equal to the melting temperature 1 ( , ) = 2 ( , ) = 0 , and the boundary conditions are: ( , 0) = ( ), (0, ) = 1 ( , ); ( , ) = 2 ( , ).
A complex Fourier series expansion is proposed as a solution within each phase ( , ) Where 2 = and 2 2 + = 0 is a necessary condition. Equation (5) is valid at any phase within the PCM. The linear part 0 + 1 is added since it is an expression that satisfies the heat equation and is linearly independent of the rest of the base functions used in this work [20]. Two possibilities of real functions can be considered. The first case, = − 2 2 where is a real value, the temperature is given by: or using sinusoidal functions The exponential decreasing time behavior of these temperature distributions may be neglected when the system is subjected to time-indepndent boundary conditions. Additionally, a second possibility occurs when = 0 + 1 is a complex number. The base functions in this case, are given by ( , ) = ( 0 + 1 ) / −( 0 + 1 ) 2 2 .
Equations (7) and (11) can be conveniently combined in order to obtain a more complex set of expressions. Additionally, the proposed solutions are not unique. For example, if is a constant then ( − , ) also satisfies the heat equation.

Homogeneous Dirichlet (Isothermal) boundary conditions
The case of a one-dimensional PCM layer with two phases, where the temperature is kept constant at the ends (0, ) = 1 and ( , ) = 2 , is considered in this section. The initial profile is given by Two set of functions are used (liquid and solid parts). Temperature is described by the equation (7) in liquid medium 0 , 1 , , coefficients, and the translated temperature in solid medium with 0 , 1 , , coefficients. At long times, the temperature is described only by the linear part. Applying the conditions at the extremes and considering ( , ) = ( , ) = 0 , the unknown coefficients are obtained: With 0 is the initial interface position, is the asymptotic time limit of the interface position and Equation (2) that is determined through local thermal balance at the liquid-solid interface is used to estimate the time evolution of the interface. Using the results previously discussed, the equation of motion for the interface takes the following form This equation is solved numerically with the RK4 method; therefore, the proposed method provides semi-analytical solutions to the dynamics of the phase transition. Fig. 1 shows the interface motion when the system is subjected to isothermal (Dirichlet) boundary conditions. The result is practically the same as obtained in Ref. [21].  The temperature field estimated with the Fourier method introduced in this work is shown in Fig. 2.  Fig. 2. Temperature field of a PCM wallboard subjected to isothermal (Dirichlet) boundary conditions, and according to the solutions estimated with the proposed method.

Steady state periodic regime
In this scenario, the temperature at the left boundary of the system, oscillates in time with a period 2 / , the temperature at the right boundary of the solid phase 2 is constant, and the temperature at the interface is equal to the melting temperature 0 of the PCM (19) Close to the periodic steady state regime, the set of functions used in the previous section is not required since these functions decay exponentially in time. In this regime, we use the set of functions shown in equation (10). The temperature distribution within the solid phase must satisfy: where the constants 0 , 1 are given by equation (14). In general, the Fourier series for the temperature profile at the left boundary when = 0 is given by When the temperature at the left boundary has the form 1 = 10 + 1 sin( ) , the final expression for the temperature in liquid phase is obtained by applying the boundary conditions as follows: ( , ) = 0 + 1 + 1 1 ( , ) + 2 2 ( , ) + 3 3 ( , ) + 4 4, ( , ), where the constants 0 , 1 are given by equation (13), and: with = √ 2 and = 4 − 2 2 cos(2 ) + 1.
The thermal balance at the liquid-solid interface is used in a similar way to that discussed in the previous section, in order to determine the interface dynamics. Fig. 3 shows the time behavior of the interface position for the salt hydrate 2 6 2 . In this case, only one functions set is used and a 5 K temperature amplitude is considered.

PCM layer under periodic boundary conditions: Transient and steady state regimes
Similarly, to the previous case, the boundary conditions are shown through equation (19). The temperature within the solid phase is shown by equation (20) with the constants given through equation (14). The temperature at the left boundary has a sinusoidal oscillation of the form 1 = 10 + 1 sin( ) ; therefore, the temperature in liquid phase is given by: where the constants 0 , 1 are given by (13), 1 , 2 , 3 and 4 by (23), and by (15). The equation of motion for the interface position is solved by using a fourth order Runge Kutta (RK4) method, in a similar fashion that in previous sections. Fig. 3 shows transient and steady state behavior of the liquid-solid interface position . The thermal performance of the PCM wallboard can be assessed by estimating the thermal energy that penetrates the PCM layer in a given time window. The thermal energy that penetrates the system, which is equal to the energy that has to be removed by an AC unit, can be obtained through the thermal flux at the right boundary = as follows: where ( ) is the total amount of thermal energy per unit area that has penetrated the wallboard up to certain time value . The thermal energy ( ) obtained through the Fourier method introduced in sections 2 and 4 is shown in Fig. 5. The thermal energy is the amount of heat that penetrates through a 10 cm thick PCM layer with a cross section of = 1 m 2 .

Discussion
The solutions estimated through the proposed method, which is based on a complex Fourier series expansion of the temperature profiles, are consistent with those obtained with the previously established HBIM. Both methods provide approximate solutions for the interface position in the transient regimes as shown in Fig. 1 and Fig. 3. Finite size effects are illustrated with the steady state values, which are well reproduced through both methods. Fig. 3 shows that the proposed method in this work can predict the steady state periodic regime in the same way as the phasors based method introduced in Ref. [18]. The time dependent behavior of the interface position in the transient regime according to the complex Fourier method introduced in this work is also obtained. Then, the predictions determined with the proposed method are more general. Daily temperature variations are cyclic during a 24-hour period, but they depend on weather conditions. Sine or cosine temperature variations are not realistic and proposed methods need to incorporate the transient response of the system. The thermal performance of the PCM based wallboard is shown in Fig.  5 when the system is subjected to the periodic boundary conditions at the left boundary 1 = 10 + 1 sin( ) , and the homogeneous Dirichlet boundary conditions (0, ) = 1 and ( , ) = 2 . The thermal performance is determined by the solution to the interface position as shown through equation (24). Fig. 5 shows the approximate solutions according to the proposed method for the amount of thermal energy that penetrates through the PCM wallboard. Small differences are found between the solutions with periodic boundary conditions and homogeneous Dirichlet (isothermal) boundary conditions. The small relative difference is due to the low amplitude temperature oscillations during the time window of the simulation and the type of PCM. Increasing the amplitude of the temperature oscillations requires to consider the formation of several fronts. Depending on the melting temperature of the PCM, the formation of several liquidsolid fronts can take place when the ambient temperature oscillates about 0 . The formation of several fronts will also depend on the amplitude of temperature variations. The proposed method needs to be generalized in order to incorporate the effects of multiple front formation, on the thermal performance of PCM wallboards.

Conclusions
In this work, Fourier series are used to analyze the temperature behavior of a one-dimensional sample with a liquid-solid phase transition. Several results are found that essentially coincide with the results obtained with the HBIM and MFD methods discussed in Ref. [21]. A basis of Fourier functions has been described that allows analyzing the behavior of the temperature distribution in a two-phase sample where one of the extremes has periodic conditions in temperature. The solution is only valid for climates where the temperature day and night is always higher than the melting temperature. The nature of the boundary conditions imposed on the sample implies an asymptotic behavior that can be predicted, allowing to find the amount of liquid and solid that will be oscillating in the sample. The results obtained in [2][3] are generalized to consider the transient part of the solution.
The method used to solve the energy balance equation at the interface does not take into consideration volumetric effects due to the density difference between the liquid and solid phases. Additionally the volume changes produced by thermal expansion are not considered in this work. The present method will be generalized to incorporate volume changes during the phase transition by including total mass as a constant of the motion and analyze the impact of these effects on the thermal performance of PCMs.