A numerical solution that determines the temperature field inside phase change materials: application in buildings.

Ciulla, Giuseppina ; Brano, Valerio Lo ; Messineo, Antonio 等

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; Ciulla et al. 2010; 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 thermo-physical 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, 2008; Kaklauskas 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 (Ginevicius 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 (Ibanez 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).

1. 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 solid-liquid 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:

[c.sub.p] = [([partial derivative]h/[partial derivative]T).sub.p=const] (1)

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:

[[nabla].sup.2]T = [1/[alpha]] [partial derivative]T/[partial derivative]t, (2)

where: [[nabla].sup.2] = ([[partial derivative].sup.2]/ [partial derivative][x.sup.2] + [[partial derivative].sup.2]/ [partial derivative][y.sup.2] + [[partial derivative].sup.2]/ [partial derivative][z.sup.2]); T is temperature; t is time; [alpha] is the thermal diffusivity; [alpha] = [lambda]/[rho][c.sub.p][m/[s.sub.2]]; [lambda] is the thermal conductivity; [rho] is the density; [c.sub.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.

2. 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 [degrees]C. 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 [degrees]C, 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 [degrees]C. Numerical simulations for a 120 [m.sup.2] house with PCM-enhanced wallboards reduced 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 [degrees]C (Cabeza et al. 2010).

3. 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:

[[DELTA].sup.2]T/[DELTA][x.sup.2] = 1/[alpha] [DELTA]T/[DELTA]t, (3)

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):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

The temperature gradients may be evaluated as a function of the nodal temperatures:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

Substituting Eqs (5) and (6) into (4) results are:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

For the time discretisation, dividing the observation time t of the phenomenon into a finite number of time steps p with amplitude [DELTA} t results are:

t = p x [DELTA]t. (8)

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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

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:

[[T.sup.p.sub.n+1] + [T.sup.p.sub.n-1] - 2[T.sup.p.sub.n]]/[([increment of x]).sup.2] = [1/[alpha]] [[T.sup.p+1.sub.n] - [T.sup.p.sub.n]]/[DELTA]t. (10)

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 [increment of x]/2 with a representative node placed on the surface and an internal domain with a length equal to [increment of x] 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:

Fo [less than or equal to] 1/2 ; Fo = [alpha][DELTA]t/[([increment of x]).sup.2], (11)

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.

[FIGURE 1 OMITTED]

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 non-isothermal transition, when [c.sub.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:

H = h + L x f, (12)

where: H is the total enthalpy; h is the sensible component; L is the latent fusion heat; and fis the liquid fraction. Furthermore:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (13)

where [T.sub.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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)

Thus, rearranging the previous equations:

H = h + L x f [??][partial derivative]H/[partial derivative]t = [partial derivative]h/[partial derivative]t + L [partial derivative]f/[partial derivative]t; (15)

dh = c x dT; (16)

[partial derivative]T/[partial derivative]t = [alpha] [[partial derivative].sup.2]T/[partial derivative][x.sup.2]; (17)

[partial derivative]h/[partial derivative]t = [partial derivative]/ [partial derivative]x ([[lambda]/[rho]] [partial derivative]T/[partial derivative]x) - L [[partial derivative]f/[partial derivative]t]. (18)

Figure 2 represents the schema of the thermal exchange at the surface and shows the boundary node length, where a thickness of [increment of x]/2 is used.

5. 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.sup.p+1.sub.0] - [T.sup.p.sub.0] > 0 (increasing temperature), then the discretised explicit form of the thermal balance will be:

[[rho].sub.i][c.sub.i]A [[increment of x]/2] ([T.sup.p+1.sub.0] - [T.sup.psub.0]) = [??]A ([T.sup.p.sub.[infinity]] - [T.sup.p.sub.0]) + 2[[lambda].sub.i]A/[increment of x] ([T.sup.p.sub.1] - [T.sup.p.sub.0]) + [q.sup.P.sub.0]A, (19)

where: [rho] is the density of the control volume; c is the specific heat; A is the area; [DELTA]x is the thickness of the domain; [T.sup.p+1.sub.0] is the temperature of the superficial node at the present time; [T.sup.p.sub.0] is the temperature of the superficial node at the previous time; [??] is the convective heat transfer coefficient; [T.sup.p.sub.[infinity]] is the temperature of the air at the previous time; [lambda] is the thermal conductivity; [T.sup.p.sub.1] is the temperature of the first internal node following the superficial temperature at the previous time and [q.sup.p.sub.0] is a generic external heat flux at the previous time.

[FIGURE 2 OMITTED]

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:

[[rho].sub.i] = [rho]([T.sup.p.sub.i]); [[lambda].sub.i] = [lambda]([T.sup.p.sub.i]); [c.sub.i] = c([T.sup.p.sub.i]). (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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (21)

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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (22)

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.sub.m] as long as the liquid fraction is. 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.sub.0] = [T.sub.m], and the liquid fraction is:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (23)

When the phase changing is occurring, [T.sub.0] = [T.sub.m], and the liquid fraction is:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (24)

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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (25)

When the transition of the substance has ended and is fully liquid, the superficial temperature is:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (26)

6. 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:

[[rho].sub.i][c.sub.i]A[increment of x]([T.sup.p+1.sub.i] - [T.sup.p.sub.i]) = [[lambda].sub.i]A/[increment of x]([T.sup.p.sub.i-1] - [T.sup.p.sub.i]) + [[lambda].sub.i]A/[increment of x] ([T.sup.pi+.sub.1] - [T.sup.p.sub.i]) + [q.sup.p.sub.i]A. (27)

Or in enthalpy notation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (28)

Solving the equation for the control volume, the enthalpy of the internal node at time p +1 is:

[h.sup.p+1.sub.i] - [h.sup.p.sub.i] + [[lambda].sub.i][DELTA]t/[[rho].sub.i][DELTA][x.sup.2]([T.sup.p.sub.i-1] - 2[T.sup.p.sub.i] + [T.sup.p.sub.i+1]) + [q.sup.p.sub.i][DELTA]t/[[rho].sub.i][increment of x] - L ([f.sup.p+1.sub.i] - [f.sup.p.sub.i]). (29)

When there is no phase change, the time variation of the liquid fraction is null, and the temperature of the internal node is:

[T.sup.p+1.sub.i] = [T.sup.p.sub.i] + [[lambda].sub.i][DELTA]t/[[rho].sub.i][c.sub.i][DELTA][x.sup.2]([T.sup.p.sub.i-1] - 2[T.sup.p.sub.i] + [T.sup.p.sub.i+1]) + [q.sup.p.sub.i][DELTA]t/[[rho].sub.i][c.sub.i][increment of x]. (30)

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.sub.m], and the liquid fraction is:

[f.sup.p+1.sub.i] = [f.sup.p.sub.i] + [[lambda].sub.i][DELTA]t/[[rho].sub.i]L[DELTA][x.sup.2] ([T.sup.p.sub.i-1] - 2[T.sup.p.sub.i] + [T.sup.p.sub.i+1]) [q.sup.p.sub.i][DELTA]t/[[rho].sub.i]L[increment of x] - [c.sub.i]/L ([T.sub.m] - [T.sup.p.sub.i]) (31)

When the phase change is occurring, [T.sub.i] = [T.sub.m], and the liquid fraction is:

[f.sup.p+1.sub.i] = [f.sup.p.sub.i] + 2 [[lambda].sub.i][DELTA]t/[[rho].sub.i]L[DELTA][x.sup.2] ([T.sup.p.sub.i-1] - 2[T.sup.p.sub.i] + [T.sup.p.sub.i+1]) + [q.sup.p.sub.i][DELTA]t/[[rho].sub.i]L[increment of x]. (32)

When the phase change has just ended, the time variation of the liquid fraction is null, and the temperature of the internal node is:

[T.sup.p+1.sub.i] = [T.sup.p.sub.i] + [[lambda].sub.i][DELTA]t/[[rho].sub.i][c.sub.i][DELTA][x.sup.2] ([T.sup.p.sub.i-1] - 2[T.sup.p.sub.i] + [T.sup.p.sub.i+1]) + [q.sup.p.sub.i][DELTA]t/[[rho].sub.i][c.sub.i][increment of x] - L(1 - [f.sup.p.sub.i])/[c.sup.i]. (33)

When the phase change has already ended and the substance is fully liquid, the temperature of the internal node is:

[T.sup.p+1.sub.i] = [T.sup.p.sub.i] + [[lambda].sub.i][DELTA]t/[[rho].sub.i][c.sub.i][DELTA][x.sup.2] ([T.sup.p.sub.i-1] - 2[T.sup.p.sub.i] + [T.sup.p.sub.i+1]) + [q.sup.p.sub.i][DELTA]t/[[rho].sub.i][c.sub.i][increment of x]. (34)

This set of equations permits us to solve the one-dimensional 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.

7. 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 1985, 1990). This test concerns a specific case of phase change in a one-dimensional geometry, where a pure, liquid substance initially at 2[degrees]C occupies a semi-infinite space. At time t = 0, the surface at x = 0 is characterised by a temperature of -4[degrees]C, which is less than the phase change temperature of the substance fixed at [T.sub.m] = 0[degrees]C.

For t >0, 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 = X(t). The analytical solution to this classical problem (Rubestein 1971; Onyejekwe, O. O., Onyejekwe, O. N. 2011) is described by the following equation:

X (t) = 2 [beta][([alpha]t).sup.1/2], (35)

where [alpha] is the thermal diffusivity, and the constant [beta] can be evaluated as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (36)

At each point, the temperature is equal to:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (37)

The thermal properties of the material are assumed to be constant and the same for both the solid and liquid phase and are:

[T.sub.m] = 0[degrees]C; T([infinity], t)= 2[degrees]C; T(x, 0) = 2[degrees]C;

T(0, t) = -4[degrees]C; [lambda] = 2 [W/mK];

c = 2.5 x [10.sup.6] J/kgK; [rho] = 1 kg/[m.sup.3]; L = 100 x [10.sup.6] J/kg.

At these conditions, the phase change at x = 0.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.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

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 ofa 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.sup.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.

doi: 10.3846/13923730.2013.778212

References

Alawadhi, E. M. 2008. Thermal analysis of a building brick containing phase change material, Energy and Buildings 40(3): 351-357. http://dx.doi.org/10.1016/j.enbuild.2007.03.001

Batens, R.; Jelle, B. P.; Gustavsen, A. 2010. Phase change materials for building applications: a state of art review, Energy and Buildings 42(9): 1361 -1368. http://dx.doi.org/10.1016/j.enbuild.2010.03.026

Butera, F.; Cannistraro, G.; Rizzo, G.; Yaghoubi, M. A. 1991. Simplified thermal analysis of naturally ventilated dwellings, Renewable Energy 1(5-6): 749-756. http://dx.doi.org/10.1016/0960-1481(91)90023-I

Beccali, G.; Cellura, M.; Lo Brano, V.; Orioli, A. 2005a. Is the transfer function method reliable in a European building context? A theoretical analysis and a case study in the south of Italy, Applied Thermal Engineering 25(2-3): 341 -357. http://dx.doi.org/10.1016/j.applthermaleng.2004.06.010

Beccali, G.; Cellura, M.; Lo Brano, V.; Orioli, A. 2005b. Single thermal zone balance solved by transfer function method, Energy and Buildings 37(12): 1268-1277. http://dx.doi.org/10.1016/j.enbuild.2005.02.010

Cabeza, L. F.; Castell, A.; Medrano, M.; Martorell, I.; Perez, G.; Fernandez, I. 2010. Experimental study on the performance of insulation materials in Mediterranean construction, Energy and Buildings 42(5): 630-636. http://dx.doi.org/10.1016/j.enbuild.2009.10.033

Cardona, E.; Piacentino, A. 2004. A validation methodology for a combined heating cooling and power (CHCP) pilot plant, Journal of Energy Resources Technology ASME 126(4): 285-292. http://dx.doi.org/10.1115/L1803849

Cellura, M.; Campanella, L.; Ciulla, G.; Guarino, F.; Lo Brano, V.; Cesarini, D. N.; Orioli, A. 2011. The redesign of an Italian building to reach net zero energy performances: a case study of the SHC Task 40 ECBCS Annex 52 (Report), ASHRAE Transactions 117, Part 2, 331-339.

Ciulla, G.; Lo Brano, V.; Orioli, A. 2010. A criterion for the assessment of the reliability of ASHRAE conduction transfer function coefficients, Energy and Buildings 42(9): 1426-1436. http://dx.doi.org/10.1016/j.enbuild.2010.03.012

de Gracia, A.; Barreneche, C.; Farid, M. M.; Cabeza, L. F. 2011. New equipment for testing steady and transient thermal performance of multilayered building envelopes with PCM, Energy and Buildings 43(12): 3704-3709. http://dx.doi.org/10.1016/j.enbuild.2011.10.010

Ginevicius, R.; Podvezko, V.; Raslanas, S. 2008. Evaluating the alternative solutions of wall insulation by multi-criteria methods, Journal of Civil Engineering and Management 14(4): 217-226. http://dx.doi.org/10.3846/1392-3730.2008.14.20

Halawa, E.; Bruno, F.; Saman, W. 2005. Numerical analysis of a PCM thermal storage system with varying wall temperature, Energy Conversion and Management 46(15-16): 2592-2604. http://dx.doi.org/10.1016/j.enconman.2004.11.003

Ibanez, M.; Lazaro, A.; Zalba, B.; Cabeza, L. F. 2005. An approach to the simulation of PCMs in building applications using TRNSYS, Applied Thermal Engineering 25(11-12): 1796-1807.

IEA. 2005. Annex 17: Advanced thermal energy storage through phase change materials and chemical reactions -feasibility studies and demonstration projects. Final Report. International Energy Agency Implementing Agreement on Energy Conservation through Energy Storage, 2005.

Incropera, F.; DeWitts, D. 2001. Fundamentals of heat mass transfer. 5th ed. Wiley. 1048 p.

Izquierdo-Barrientos, M. A.; Belmonte, J. F.; Rodriguez-Sanchez, D.; Molina, A. E.; Almendros-Ibanez, J. A. 2012. A numerical study of external building walls containing phase change materials (PCM), Applied Thermal Engineering 47: 73-85. http://dx.doi.org/10.1016/j.applthermaleng.2012.02.038

Kaklauskas, A.; Rute, J.; Zavadskas, E. K.; Daniunas, A.; Pruskus, V.; Bivainis, J.; Gudauskas, R.; Plakys, V. 2012. Passive house model for quantitative and qualitative analyses and its intelligent system, Energy and Buildings 50: 7-18. http://dx.doi.org/10.1016/j.enbuild.2012.03.008

Kanapeckiene, L.; Kaklauskas, A.; Zavadskas, E. K.; Raslanas, S. 2011. Method and system for multi attribute market value assessment in analysis of construction and retrofit projects, Expert Systems with Applications 38(11): 14196-14207.

Koschenz, M.; Lehmann, B. 2004. Development of a thermally activated ceiling panel with PCM for application in lightweight and retrofitted buildings, Energy and Buildings 36(6): 567-578. http://dx.doi.org/10.1016/j.enbuild.2004.01.029

Liu, H.; Awbi, H. B. 2009. Performance of phase change material boards under natural convection, Building and Environment 44(9): 1788-1793. http://dx.doi.org/10.1016/j.buildenv.2008.12.002

Mottaghy, D.; Dijkshorn, L. 2012. Implementing an effective finite difference formulation for borehole heat exchangers into a heat and mass transport code, Renewable Energy 45: 59-71. http://dx.doi.org/10.1016/j.renene.2012.02.013

OECD IEA. 2011. World energy outlook 2011. Paris: OECD International Energy Agency. 696 p.

Onyejekwe, O. O.; Onyejekwe, O. N. 2011. Numerical solutions of the one-phase classical Stefan problem using an enthalpy green element formulation, Advances in Engineering Software 42(10): 743-749. http://dx.doi.org/10.1016/j.advengsoft.2011.05.012

Principi, P.; Di Perna, C.; Carbonari, A.; Ferrini, M. 2004. Validation of numeric finite element calculation model through the laboratory experimentation of artificial programmable inertia walls, in Proc. of the XXXII IAHS Word Congress on Housing, 21-25 September, 2004, Trento, Italy.

Pupeikis, D.; Stankevicius, V.; Burlingis, A. 2010. The effect of the Fourier number on calculation of an unsteady heat transfer of building walls, Journal of Civil Engineering and Management 16(2): 298-305. http://dx.doi.org/10.3846/jcem.2010.34

Roth, K.; Westphalen, D.; Brodrick, J. 2007. PCM technology for building materials, ASHRAE Journal 49(7): 129-131.

Rubestein, L. 1971. The Stefan problem (translations of mathematical monographs), Vol. 27. American Mathematical Society. 419 p.

Shilei, L.; Guohui, F.; Neng, Z.; Li, D. 2007. Experimental study and evaluation of latent heat storage in phase change materials wallboards, Energy and Buildings 39(10): 1088-1091. http://dx.doi.org/10.1016/j.enbuild.2006.11.012

Tyagi, V. V.; Buddhi, D. 2007. PCM thermal storage in buildings: a state of art, Renewable and Sustainable Energy Reviews 11(6): 1146-1166. http://dx.doi.org/10.1016/j.rser.2005.10.002

Tupenaite, L.; Zavadskas, E. K.; Kaklauskas, A.; Turskis, Z.; Seniut, M. 2010. Multiple criteria assessment of alternatives for built and human environment renovation, Journal of Civil Engineering and Management 16(2): 257-266. http://dx.doi.org/10.3846/jcem.2010.30

Voller, V. R.; Cross, M. 1981. Accurate solutions of moving boundary problems using enthalpy method, International Journal of Heat and Mass Transfer 24(3): 545-556. http://dx.doi.org/10.1016/0017-9310(81)90062-4

Voller, V. R. 1987. A heat balance integral method based on an enthalpy formulation, International Journal of Heat and Mass Transfer 30(3): 604-607. http://dx.doi.org/10.1016/0017-9310(87)90275-4

Voller, V. R. 1990. Fast implicit finite-difference method for the analysis of phase change problem, Numerical Heat Transfer, Part B: Fundamentals 17(2): 155-169. http://dx.doi.org/10.1080/10407799008961737

Waqas, A.; Kumar, S. 2011. Thermal performance of latent heat storage for free cooling of buildings in a dry and hot climate: an experimental study, Energy and Buildings 43(10): 2621-2630. http://dx.doi.org/10.1016/j.enbuild.2011.06.015

Weinlaeder, H.; Koerner, W.; Heidenfelder, M. 2011. Monitoring results of an interior sun protection system with integrated latent heat storage, Energy and Buildings 43(9): 2468-2475. http://dx.doi.org/10.1016/j.enbuild.2011.06.007

Zavadskas, E. K.; Kaklauskas, A.; Turskis, Z.; Tamosaitiene, J. 2008. Selection of the effective dwelling house walls by applying attributes values determined at intervals, Journal of Civil Engineering and Management 14(2): 85-93. http://dx.doi.org/10.3846/1392-3730.2008.14.3

Zavadskas, E. K.; Ustinovicius, L.; Turskis, Z.; Ambrasas, G.; Kutut, V. 2005. Estimation of external walls decisions of multistorey residential buildings applying methods of multicriteria analysis, Technological and Economic Development of Economy 11(1): 59-68.

Zivkovic, B.; Fujii, I. 2001. Analysis of isothermal phase change of phase change material within rectangular and cylindrical containers, Solar Energy 70(1): 51-61. http://dx.doi.org/10.1016/S0038-092X(00)00112-2

Giuseppina Ciulla (a), Valeric. Lo Brano (a), Antonio Messineo (b), Giorgia Peri (a)

(a) Dipartimento di Energia, ingegneria dell'Informazione e Modelli matematici, Universita degli Studi di Palermo, viale delle scienze Edificio 9, 90128, Palermo, Italy

(b) Universita Kore di Enna, cittadella Universitaria, 94100, Enna, Italy

Received 13 Jun. 2012; accepted 03 Oct. 2012

Corresponding author: Valerio Lo Brano

E-mail: lobrano@dream.unipa.it

Giuseppina CIULLA. Graduated in Engineering at the University of Palermo in 2006; PhD in Energy (2010). She has worked in the educational activities of the department and is an active member of several national and international research groups (Net Zero Energy Buildings--IEA Task 40; and Compact Thermal Energy Storage-Task42-IEA). Author of over 20 papers published in journals and national and international conferences. Other scientific interests are: processing of meteorological data, monitoring of solar systems, analysis of energy performance of PV systems. During the a.a. 2010/2011 was appointed aggregate professor of thermo-physics of buildings at the Faculty of Engineering and of Environmental Technology Control at the Faculty of Architecture, University of Palermo. The research activities were primarily intended to study: thermal dynamic simulation of buildings, thermo-physics and energy certification of buildings and use of renewable energy sources in urban areas.

Valerio LO BRANO. M.S. (Nuclear Engineering) Palermo University, Palermo, Italy, 1998. Thesis topic: safety analysis of AP600 nuclear reactor with Fault Tree approach, PhD (Environmental Applied Physics) Palermo University, Palermo, Italy, 2003. Application of neural network to heat transfer problems. Since 2011 Associate Professor; From 2004 Assistant Professor teaching Renewable Energy Sources and Civil Use of Energy, Faculty of Engineering, Palermo University, Palermo, Italy. From 2003 to 2004 Project Specialist, Department of Energy and Environmental Researches, Palermo University. Responsibilities included research and writing on the technical and energy aspects of building materials; analysis of buildings-plants coupling; application of neural techniques to energy demands forecasting. Preparation of requests for proposals; network administration and management of grants in the field of Information Technology. Author of more than 70 publications in international journals and conferences. Research interests: energy planning, lifecycle assessment, heat transfer in building structure, thermal dynamic simulation of buildings, fuzzy cognitive maps, urban weather monitoring.

Antonio MESSINEO. Graduated cum laude in Mechanical Engineering (major in Energy Conversion) in 2001 at the University of Palermo. PhD in Energetics (2007), he is Assistant Professor of Environmental Applied Physics at Faculty of Engineering and Architecture of the Kore University of Enna. He's Referee of different international journals and B2 Commission (Refrigerating Equipment) Member of International Institute of Refrigeration (IIFIIR), Paris. Over 60 papers published on national and international journals and/or presented at scientific congresses.

Giorgia PERI. MSc Degree in Environmental Engineering with full marks cum laude, gained at the University of Palermo in 2008. Research period at the Sustainable Engineering Department of the Technische Universitat Berlin in Berlin (Germany) from 2010 till 2011 (1 year, 2 months). PhD in Environmental Applied Physics gained at the University of Palermo in 2012 along with the "Doctor Europeaus" certification. Other scientific research interests include: the evaluation and assessment of the overall environmental performance of buildings and their rating; the using of environmental quality brands for sustainable tourism. Over 20 papers published on national and international journals and/or presented at scientific congresses. The main scientific research interest regards green roofs (planted roofs or rooftop gardens) with a specific focus on both the analysis of the green roof thermal properties and investigation of its energy performance, and the assessment of the green roof life cycle from an environmental, economic and social perspective.