A NUMERICAL SOLUTION THAT DETERMINES THE TEMPERATURE FIELD INSIDE PHASE CHANGE MATERIALS: APPLICATION IN BUILDINGS

The use of novel building materials that contain active thermal components would be a major advancement in achieving significant heating and cooling energy savings. In the last 40 years, Phase Change Materials or PCMs have been tested as thermal mass components in buildings, and most studies have found that PCMs enhance the building energy performance. The use of PCMs as an energy storage device is due to their relatively high fusion latent heat; during the melting and/or solidification phase, a PCM is capable of storing or releasing a large amount of energy. PCMs in a wall layer store solar energy during the warmer hours of the day and release it during the night, thereby decreasing and shifting forward in time the peak wall temperature. In this paper, an algorithm is presented based on the general Fourier differential equations that solve the heat transfer problem in multi-layer wall structures, such as sandwich panels, that includes a layer that can change phase. In detail, the equations are proposed and transformed into formulas useful in the FDM approach (finite difference method), which solves the system simultaneously for the temperature at each node. The equation set proposed is accurate, fast and easy to integrate into most building simulation tools in any programming language. The numerical solution was validated using a comparison with the Voller and Cross analytical test problem.


Introduction
Within the framework of directives set by the European Union to reduce green-house gas emissions, the actions to reduce the energy impact of building sectors are extremely important. In 2009, the total primary energy demand in Europe and the world was 1,654 and 12,132 Mtoe (million tons of oil equivalent), respectively. To adopt these new global policies and measures aimed at rationalising energy consumption, the world primary energy will increase from 12,132 Mtoe in 2009 to 16,961 Mtoe in 2035 ('40%); in the building sector, the global energy demand will increase by 34%, assuming that the sector's share of total final consumption will remain stable. Simulation models predict that the primary energy demand in the European Union will increase by 5% from the 2009 levels by 2035; however, natural gas demand will increase by 24% over this period, which will largely be due to the power sector and also the heating in buildings (OECD IEA 2011).
In scientific literature, a simplified thermal analysis of buildings has been approached several times due to its relevance for technicians who need simple and accurate methods to simulate the energy behaviour of buildings (Butera et al. 1991;Principi et al. 2004;Cellura et al. 2011).
It became clear that reducing the heat losses of buildings or in general, the total energy consumption of buildings, is a key objective for many nations characterised by high energy consumption due to the old housing stock that requires improvements in its thermophysical characteristics and the need to minimise the size of thermal plants (Cardona, Piacentino 2004;Beccali et al. 2005a, b). Indeed, approximately 80% of the total thermal energy consumption of buildings is consumed because of heat losses (Pupeikis et al. 2010). The costs attributed to the thermal energy needs in buildings is strongly influenced by the characteristics of the building envelope (Zavadskas et al. 2005(Zavadskas et al. , 2008Kaklauskas et al. 2012;Kanapeckiene et al. 2011).
The poor quality of thermal insulation of old residential and public buildings makes the renovation of these buildings an urgent problem. In these cases, renovation is vitally important not only because of the immediate consequences, i.e. reduced energy consumption and improved state of buildings, but also because of the positive external effects, such as an increased quality of life and reduced climate change (Tupenaite et al. 2010). Studies on the Lithuanian housing stock built during the Soviet year, which were characterised by high values of thermal transmittance, have shown that the highest economic effect can be obtained by insulating the external walls. In this case, the economic saving is approximately double compared with the savings realised if the windows are replaced (Ginevičius et al. 2008).
Phase change materials (PCMs) represent a possible solution that may reduce peak loads and thermal energy consumption in buildings due to their good insulation properties and their thermal inertia related to the phase change phenomenon (Tyagi, Buddhi 2007;Halawa et al. 2005;Liu, Awbi 2009). Even in hot climates, thermal insulation materials are extensively used in building structures to reduce the heat flow into the building's indoor space by providing an effective thermal resistance to the heat flow (Alawadhi 2008). Indeed, many authors have investigated the use of PCMs in building elements and structures (Ibáñ ez et al. 2005;de Gracia et al. 2011;Waqas, Kumar 2011;Weinlaeder et al. 2011).
In this work, authors present an algorithm based on the general Fourier differential equations that solves the heat transfer problem in multi-layer wall structures, such as sandwich panels, which includes a layer that changes its phase. The equations are proposed and transformed into formulas useful in the FDM approach (finite difference method).

Generalities of phase change materials
Phase change materials, commonly known as PCMs, are materials that are capable of storing large amounts of heat because they have a high latent heat of liquefaction and are solid at ambient temperature, i.e. PCMs absorb or release excessive thermal energy by undergoing a phase transition. There are different substances with a high latent heat property that allows the storage of large amounts of heat during a solidliquid phase transition. For pure substances, this phase change occurs at a constant temperature, and for certain materials, the process of melting and solidifying can be repeated over a large number of cycles with no change in their physical or chemical properties.
The main features of the PCMs materials are the following: Á high thermal energy storage capacity; Á no performance decay; Á no toxicity; Á chemical neutrality.
From the thermal point of view, the specific heat of PCMs varies as a function of its temperature, reaching its maximum value during the phase transition. For thermal problems involving phase changes, it is more suitable to use enthalpy to define specific heat at a constant pressure: Generally, during a generic phase transition, the isobaric phase change is isothermal; otherwise, for impure substances, the temperature variation remains within a very limited range.
To determine the thermal field in a dynamic regime domain, we can use the Fourier equation: where: ; l is the thermal conductivity; r is the density; c p is the specific heat at a constant pressure.
In particular cases, the analytical method may be used to obtain exact mathematical solutions for steady conduction problems. These solutions have been generated for an assortment of simple geometries and boundary conditions and are well documented in the literature. However, in most cases, thermal problems involve geometries and/or boundary conditions that preclude such solutions (Shilei et al. 2007).
The authors have chosen to use the FDM to solve the problem in a one-dimensional geometry. This method represents one of the simplest ways to search for a numerical solution of the heat transfer problem because it is possible to reduce a partial differential equation system into an algebraic equation system (Mottaghy, Dijkshorn 2012). The resolution algorithm can be coded, and the thermal analysis can be performed easily and quickly. Due to the phase change, as explained below, it is convenient to use an enthalpy formulation in the heat transfer problem.

PCMs for heating and cooling in buildings
By inserting phase change materials into the walls of the building envelope, the structure can store large amounts of energy while maintaining the indoor temperature within a small variation. Therefore, PCMs with a melting range near the desired room temperature should be chosen.
There are several examples that can describe the integration of a PCM layer in a building envelope, where the position of the PCM strongly influences the thermal behaviour of the wall (Izquierdo-Barrientos et al. 2012): Á A PCM layer can be placed within the wall near the external surface; in this case, during hot and sunny days, the PCM stores a great amount of thermal energy. Later, during the night, the energy is released, where it primarily flows to the outside of the building. This configuration is used in a building to diminish the thermal load during the summer; Á A PCM layer can be placed within the wall near the internal surface to increase the thermal mass of the buildings. In this way, the PCM layer stores energy when the indoor air temperature is higher than the melting temperature. Conversely, the PCM releases stored energy when the indoor air temperature is lower than the melting temperature. This configuration allows the building to retain more energy in the interior when the highest thermal load occurs during the winter; Á Furthermore, a PCM layer can be placed as an underfloor system or in a ceiling panel (Koschenz, Lehmann 2004). In this case, it is possible to reduce the peak temperature, which helps to maintain the internal thermal comfort conditions. The latent heat required to induce a PCM to change phase is much higher than the specific heat of conventional building materials, e.g., the latent heat capacity of a wallboard with 30% PCM by weight is approximately five times greater than the heat stored by conventional wallboards for a variation of 5.5 8C. Consequently, the PCM effectively increases the thermal mass of the building material when temperatures rise above or below the transition temperature (IEA 2005).
The impact of PCMs in building in terms of energy saving varies greatly between different building types. Climate also has a strong influence: in a moderate climate, where the outdoor temperature is often near the PCM transition temperature, will realise greater savings than climates where the temperatures during the cooling season remain predominantly above or below the transition temperature.
Simulations for a building located in Dayton, Ohio, where 10% of the 305 mm concrete walls contained PCM characterised by a transition temperature of 26 8C, resulted in a 13% reduction in annual cooling loads. The application to a steel roof with 25 mm thick panels with 28% of the same PCM decreased the annual cooling loads through that building component by 14%. Simulations of a wall frame found that a gypsum board with PCM (also at 28%) reduced the annual cooling loads through the wall by approximately 9% (Roth et al. 2007).
Another study on the effect of PCMs at reducing the heating load was performed for a case study in Helsinki, where a 13 mm thick wallboard was impregnated with 30% PCM by mass with a melting temperature of 21.4 8C. Numerical simulations for a 120 m 2 house with PCM-enhanced wallboards re-duced the supplementary heating energy by approximately 2 GJ/year, or approximately 6%, which is an annual benefit of $34 for the present cost of heat, approximately 17 $/GJ, which results in a payback time of 18 years, assuming a cost of 1.5 $/kg for industrial-grade fatty acids (Batens et al. 2010).
Regarding the Mediterranean climate, different types of PCMs have been experimentally tested in Spain. The results of the monitoring campaign showed a reduction in consumption related to a cooling load between 15 and 17% using a PCM with a transition temperature between 26 and 28 8C (Cabeza et al. 2010).

The finite difference method
The finite difference method is based on the approximation of the spatial and temporal derivatives of the heat diffusion equation with relations between finite differences of space and time. In a one-dimensional system, the heat diffusion equation for the finite difference method becomes: and a double discretisation (space and time) is necessary.
In contrast to an analytical solution, a numerical solution allows the determination of the temperature only at discrete points. The original domain must be subdivided into a finite number of parts, assigning each a reference point at its centre (node).
Each n node represents a certain region, and its temperature is a measure of the average temperature of the region. The value of this derivative at the n nodal point may be approximated as (Incropera, DeWitts 2001): The temperature gradients may be evaluated as a function of the nodal temperatures: Substituting Eqs (5) and (6) into (4) results are: For the time discretisation, dividing the observation time t of the phenomenon into a finite number of time steps p with amplitude^t results are: Therefore, the time derivative present in the second member of the heat diffusion equation can be replaced with the finite difference approximation using the ''central difference'' operator: where: p'1 is the current time step and p is the previous time step. Finally, replacing Eqs (7) and (9) into (3) results are: In other words, an explicit finite difference method can calculate the state of a system at a future time from the state of the system at the current time.
Beginning from the approach by Zivkovic and Fujii (2001), in which a simple, implicit computational model for an isothermal phase change was presented, the authors developed a modified algorithm based on an explicit finite difference formulation of the heat equation. To this aim, a forward difference at time t and a first-order central difference for the space derivative at position x was used. To use this method for building envelope walls, two sets of recursive equations were developed for two types of spatial domains: a boundary domain characterised by a length equal to Dx/2 with a representative node placed on the surface and an internal domain with a length equal to Dx with a representative node placed in the middle. The novelty concerning the definition of the boundary domain allows us to directly consider the radiative and convective heat transfer occurring on the surface in the recursive equation. Furthermore, the presence of a generic heat flow q is always considered to correctly simulate the presence of a possible active layer in the wall.
The explicit approach followed by the authors can occasionally result in stability problems. In these cases, the response is affected by fluctuations that disrupt and eventually make the solution completely unreliable. There are several criteria that ensure the stability of the method, and in our case, in the definition of the time step and for the size of the domain, we considered the following condition: where Fo is the discretised Fourier number. If we consider a building envelope coupled with a PCM system, the energy balance must take into account the presence of the phase change material. Figure 1 schematically shows the energy exchanges in a building envelope-PCM system. Due to the simple geometry, it was possible to choose a one-dimensional approach that considers only heat flow orthogonal to the wall plane.
Specifically, a plausible building envelope coupled with a PCM heat storage system is composed of a succession of layers of plaster, insulation, concrete and PCM. An optional layer of air, representing a possible imperfect contact, can be supposed between the wall structure and the heat storage system. Furthermore, the two plastic layers constituting the bag that contains the PCM (bag layer) must be taken into account.
Additionally, in the algorithm, it is possible to implement the thickness of the envelope containing the PCM and a possible layer of air due to the imperfect contact of the envelope of PCM with the wall.
The wall-PCM system is presented as a multilayer plate invested by solar radiation, where heat is exchanged between the outdoor and inner environments via convection and radiation. Depending on the properties of the PCM and the amount of captured solar energy, the PCM layer can partially or completely melt during maximum insolation and returns that amount of energy during the night, with the possibility of solidifying again. The first hypothesis is that the phase change is isothermal. In the case of a nonisothermal transition, when c p as a function of temperature is known, the problem is reduced to a simple heat conduction case. The occurrence of a phase change that maintains the temperature at a given value is a more interesting case and requires the determination of the PCM liquid fraction inside the discretised domain.
This hypothesis is not far from reality because many PCMs are characterised by an isothermal phase change, and several paraffin and eutectic mixtures have a small transition temperature. 4. One-dimensional thermal analysis of an isothermal phase change using the explicit finite difference approach To obtain the finite difference approximation of the heat equation in one-dimensional form, it is possible to use the ''central difference'' operator of equation (7), where the n index identifies the position along the x-axis of the nodal point under examination and the p index denotes the time dependence; p'1 denotes the present time, and p denotes the previous time interval.
However, in the case of a phase change of the medium, it is easier to use equations that depend on the enthalpy variation. If we assume that the total enthalpy is the sum of the sensible and the latent component, we can state: where: H is the total enthalpy; h is the sensible component; L is the latent fusion heat; and f is the liquid fraction. Furthermore: where T m is the temperature of the phase change and c is the specific heat. If we suppose that the phase change is isothermal, it is possible to define the liquid fraction as: Thus, rearranging the previous equations: Figure 2 represents the schema of the thermal exchange at the surface and shows the boundary node length, where a thickness of Dx/2 is used.

Balance equations at the boundaries nodes
In the case of border nodes, it is possible to write the energy balance of the domain that pertains to these nodes. Assuming that on the external surface, there is a generic heat flux q and a convective process and that T pþ1 0 À T p 0 > 0 (increasing temperature), then the discretised explicit form of the thermal balance will be: where: r is the density of the control volume; c is the specific heat; A is the area; Dx is the thickness of the domain; T pþ1 0 is the temperature of the superficial node at the present time; T p 0 is the temperature of the superficial node at the previous time; r is the convective heat transfer coefficient; T p 1 is the temperature of the air at the previous time; l is the thermal conductivity; T p 1 is the temperature of the first internal node following the superficial temperature at the previous time and q p 0 is a generic external heat flux at the previous time.
For the sake of simplicity, in the following section of the paper, we refer to the density, thermal conductivity and specific heat without manifestly indicating that they are functions of temperature at the previous instant: (20) Using the enthalpy notation and solving the equation for the control volume, the following is the enthalpy at a superficial node at time p'1: When there is no phase change, the time variation of the liquid fraction is null, and the equation is useful to calculate the superficial temperature: In the case when the phase change has just begun, it is possible to make several observations. First, if the phase change is isothermal, during this process, the temperature is constant at T m as long as the liquid fraction is 1. Clearly, during the process, the time variation of the sensible enthalpy is null. Furthermore, an additional term has to be evaluated to consider the sensible heat that the control volume has absorbed to overcome the previous temperature and to reach the phase change temperature in one time step.
When the phase change begins, T 0 0T m , and the liquid fraction is: When the phase changing is occurring, T 0 0T m , and the liquid fraction is: When the phase change is just ending, the superficial temperature is again free to change; however, we must compute an additional term that concerns the necessary latent heat to end the phase change. The time variation of the liquid fraction becomes null, and the superficial temperature can be calculated with the following equation: When the transition of the substance has ended and is fully liquid, the superficial temperature is:

Internal node
Even in case of an internal node, it is possible to write the energy balance of the associated control volume. If assuming only conductive heat flux coming from the previous and next nodes, the thermal balance in explicit form will be: Or in enthalpy notation: Solving the equation for the control volume, the enthalpy of the internal node at time p'1 is: When there is no phase change, the time variation of the liquid fraction is null, and the temperature of the internal node is: With similar considerations, it is possible to state that when the phase change has just begun, the temperature of the internal node will be T m , and the liquid fraction is: When the phase change is occurring, T i 0T m , and the liquid fraction is: Journal of Civil Engineering and Management, 2013, 19(4): 518Á528 When the phase change has just ended, the time variation of the liquid fraction is null, and the temperature of the internal node is: When the phase change has already ended and the substance is fully liquid, the temperature of the internal node is: This set of equations permits us to solve the onedimensional thermal analysis of an isothermal phase change using the explicit finite difference approach.
The explicit method is numerically stable and convergent if the limits imposed by Eqn (11) are respected. In contrast with the method of Zivkovic and Fujii (2001), our approach is generally less numerically intensive than the implicit method.
The above equations constitute an algorithm that can easily be implemented in a software procedure that uses a common programming language, such as VB.NET.
To better understand how our algorithm works, a logic flow chart is shown in Figure 3.
To test the validity and the performance of the presented algorithm, we compared the numerical results with the analytical solutions available for a specific problem of phase change: the Voller and Cross problem.

Validation by comparison with an analytical solution
To validate the numerical solution, it is possible to compare the results with those obtained from an exact solution. The verification of the mathematical model presented was made with the test problem of Voller and Cross (Voller, Cross 1981;Voller 1985Voller , 1990. This test concerns a specific case of phase change in a one-dimensional geometry, where a pure, liquid substance initially at 2 8C occupies a semiinfinite space. At time t00, the surface at x00 is characterised by a temperature of(4 8C, which is less than the phase change temperature of the substance fixed at T m 00 8C. For t0, the separation surface between the solid and liquid moves inside the half plane. The test consists of determining the position of the separation surface between the solid and liquid phase at different times, x 0X(t). The analytical solution to this classical problem (Rubestein 1971;Onyejekwe, O. O., Onye-jekwe, O. N. 2011) is described by the following equation: where a is the thermal diffusivity, and the constant b can be evaluated as: At each point, the temperature is equal to: ; for x > X t ð Þ: The thermal properties of the material are assumed to be constant and the same for both the solid and liquid phase and are: At these conditions, the phase change at x00.25 m occurs after approximately 453,600 seconds. The proposed model described earlier is able to correctly determine this time. For greater clarity, Fig. 4 shows the comparison between the evolutions of the temperature obtained with the analytical solution of Voller and Cross and the temperature trend from the proposed numerical solution.
In particular, the chart refers to the temperature at 25 cm from the surface of a semi-plan, where the time interval chosen for the discretisation was 30 seconds. With a time resolution of less than 30 seconds, the numerical solution is closer to the analytical solution; however, the computation time becomes excessive.

Comments and conclusions
The use of PCMs in the building envelope as thermal storage has attracted great interest in the research community for many years. Applying these systems that simultaneously have thermal inertia and good insulation characteristics requires appropriate calculation tools to evaluate the thermal fields. The latter is more difficult to determine in the case of isothermal phase change.
In Italy, the energy performance evaluation of buildings built in compliance with European regulations has been assigned as a calculation procedure implemented in certified software. The calculation procedure and current certified software are not able to add PCMs in the building envelope, due to the lack of a calculation algorithm. In the paper, we developed a finite difference model that can be used to predict the thermal behaviour of a wall coupled with phase change materials.
Inspecting the graphs in Fig. 4, it is possible to observe that the numerically calculated temperature trend is in good agreement with the analytical calculated temperatures, validating the reliability of the calculation model.
The results show that the proposed model is valid and can be used to determine the thermal behaviour of a multi-layer wall that includes a phase change material. The model was validated by a comparison with the analytical solution of the classic Voller and Cross problem. The good agreement between the analytical resolution and numerical prediction has shown that the algorithm, although simplified and for a one-dimensional geometry, can be used to determine the temperature trend of a multi-layer wall with a PCM.
The proposed equation set can be immediately usable and could be implemented in a simple and easy manner in other software tools. Thanks to its simple form, the energy balance can be easily modified and adapted to a specific problem. The method is extremely flexible and could be adopted to calculate the thermal field inside a multi-layer wall of a free-floating room. Knowing the thermo-physical characteristics of the materials that compose the wall, it is possible to simulate the dynamic profile of each inner surface temperature and additionally, the inner air temperature. The proposed algorithm and the equations are fully described, which allows maximum clarity. The explanation of all the steps of the calculations makes this method that can be applied even in different situations where a few parameters or relationships are varied.
The flexibility of the proposed method is a fundamental characteristic that can improve the performance of thermal building simulations.
To improve the accuracy of the simulation, the authors are working on a more sophisticated calculation scheme using the Crank-Nicolson approach. Although this method is more complex to implement in software, this approach should allow us to obtain a good simulation using less computation time. Furthermore, an actual scale test facility is under construction on the roof of the Energy Department. This experimental system consists of a single housing structure of approximately 35 m 2 , where a monitoring system for indoor and outdoor climate conditions will be installed. The temperature monitoring across the walls will allow us to validate the proposed algorithm experimentally. The results from experimental testing will be analysed after considering the energy benefits and economic and environmental costs.