**3. Characterization of the heat exchanger**

The main objectives of this section are:


#### **3.1 Experimental set up to test PCM-air heat exchangers**

An experimental setup was designed to study different PCM-air heat exchangers (Dolado, 2011; Lazaro, 2009). A closed air loop setup was used to simulate indoor conditions. The setup design was based on the ANSI/ASHRAE STANDARD 94.1-2002 "Method of Testing Active Latent-Heat Storage Devices Based on Thermal Performance" (ANSI/ASHRAE, 2002). The setup is constituted of: 1) Inlet air conditioner allowing the simulation of different

selected because it has been a deeply studied geometry since London & Seban, 1943. It involves: 1) Easy-to-control PCM thickness, which is a crucial design factor as it allows regulating elapsed times of the melting and solidification; 2) Uniformity of the PCM thickness and, therefore, of the phase change process; 3) Simplicity of the manufacturing process (both small scale and large scale) and versatility of handling (transportation, installation ...); 4) Commercial accessibility in a wide variety of plate-shaped encapsulations

Fig. 5. PCM panels and air flow system (left); PCM-air heat exchanger (right).

Finally, the rigid metallic plate encapsulation has been selected to avoid both compatibility issues (Lazaro et al., 2006) as well as leakage problems detected previously (Lazaro, 2009)

The basics of the heat transfer in PCM are compiled by Zalba et al., 2003, and discussed in a very understanding way by Mehling & Cabeza, 2008. The authors describe the basics of the heat transfer by means of: 1) Analytical models; 2) Numerical models; 3) Modelling; 4)

An experimental setup was designed to study different PCM-air heat exchangers (Dolado, 2011; Lazaro, 2009). A closed air loop setup was used to simulate indoor conditions. The setup design was based on the ANSI/ASHRAE STANDARD 94.1-2002 "Method of Testing Active Latent-Heat Storage Devices Based on Thermal Performance" (ANSI/ASHRAE, 2002). The setup is constituted of: 1) Inlet air conditioner allowing the simulation of different

Comparison of models vs. experimental; 5) Methods to improve the heat transfer.

4. Importance of the uncertainties in measurements and their propagation.

in different materials, both metallic and plastic.

when using pouches.

**2.5 Heat transfer mechanisms: basics** 

The main objectives of this section are:

2. Gathering experimental results.

**3. Characterization of the heat exchanger** 

1. How to test a prototype of PCM-air heat exchanger.

**3.1 Experimental set up to test PCM-air heat exchangers** 

3. Analyzing data and obtaining empirical models.

operating modes (5 kW air chiller and 4.4 kW electrical resistance); 2) Air flow measurements; 3) Difference between inlet and outlet air temperature measurements (thermopile); 4) Inlet and outlet air temperature and humidity measurements; 5) PCM and air channels temperature measurements (31 thermocouples); 6) Data logger and data screening; 7) Air ducts and gates; 8) PID controller.

The energy balance of air between the prototype's inlet and outlet is utilized to evaluate the cooling (equation 4). As the main parameters are the air flow and the air temperature difference between the inlet and the outlet, the accuracy depends on the precision when measuring these parameters. The methods used are:

$$
\dot{Q} = \dot{m}\_{\text{air through HX}} \cdot \Delta \mathbf{l}\_{\text{air}} \approx \dot{m}\_{\text{air through HX}} \cdot \mathbf{c}\_{\text{p}\_{\text{air}}} \cdot \Delta T \tag{4}
$$


The reader can find more information on the experimental setup in Lazaro, 2009.

## **3.2 Two prototypes**

Two real-scale prototypes of PCM-air heat exchangers were constructed and incorporated into the experimental setup to characterize them. Initially tests were conducted with the equipment filled with bags of a hydrated salt PCM (prototype 1). Subsequently, the bags were replaced by plates of a paraffin based PCM, and the unit was tested filled with plates. These two geometries were arranged vertically and parallel to the airflow. The casing of the heat exchanger unit used in both cases was the same. PCM thickness was a critical parameter to obtain the required cooling rates (Dolado et al., 2007). Vertical position was a requirement; therefore a metallic grid was used to force PCM thickness below a maximum in vertical position. The experimental setup built to test this kind of heat exchangers is shown in Figure 6. Tests using a constant inlet air temperature setpoint were accomplished.

Figure 7 (left) shows the cooling power evolution in prototype 1. Results showed that the cooling rates were very low and the total melting times were double the melting design time (2h). Different air flow rates were tested. As it can be seen in figure 8 (left), the air flow influence on melting times and cooling rates were negligible (in the figure *HH:mm* denotes time, hours:minutes). Cooling power does not increase by a rise of air flow rates. Indicating that, contrary to what was at first designed, heat transfer by conduction inside the PCM resistance is dominant. The prototype was opened to confirm the diagnosis, and PCM

PCM-Air Heat Exchangers: Slab Geometry 437

Fig. 8. Cooling rate evolution with constant inlet air temperature in prototype 1 during tests with different air flow rates (left); in prototype 2 (right); Flow 1 is the mass flow, equals to

Prototype 2 was designed using aluminium panels filled with organic PCM. Configuration was also vertical. Air flows were parallel to the panels from top to bottom. Due to the fact that the PCM in this prototype was organic, it presented lower thermal conductivity than PCM in prototype 1. Furthermore, the total stored energy in prototype 2 was also lower than in prototype 1. The first tests were accomplished using constant inlet air temperature. As it can be seen in figure 7 (right), cooling rates were higher than those obtained with prototype

Different air flow rates in prototype 2 were tested. It was observed that it had influence on the melting time and cooling power. Figure 8 (right) shows that for the lowest air flow rate, the heat rate curve changes its shape and is more similar to prototype 1. This indicates that heat transfer by conduction inside the PCM starts to be relevant when compared to air convection. The first results of prototype 2 were satisfactory, so more test were planned to evaluate its behaviour under real conditions. Two types of experiments were accomplished: constant rise of inlet air temperature and constant heating power. Temperature rise ramps were then set into the resistances controller: results showed that faster the rise, higher cooling power and lower melting time. For constant power tests, different heating powers of electrical resistances were fixed. Results showed that prototype 2 was able to maintain a cooling capacity over 3 kW for approximately 1h 30' or approximately 1 kW for more than 3 h. This result is useful to design the optimal operation mode depending on the application.

The total energy exchanged during melting and solidification, as well as the time elapsed until total melting/solidification were determined from the heat rate curves experimentally obtained. The influence of the inlet air temperature and air flow was studied, and results showed that the continuous thermal cycling of the unit is a repetitive process: running experiments with similar conditions led to the same thermal behaviour; no degradation in the PCM properties was noticed. Pressure drop was measured for different air flows. Depending on the inlet air temperature, full solidification of the PCM could be achieved in less than 3 h for an 8 ºC temperature difference between the inlet air and the average phase

1 and the melting times were half the melting times with prototype 1.

0.34 kg/s.

**3.3 Experimental results** 

Fig. 6. Experimental installation arrangement to test PCM-air heat exchangers prototypes.

Fig. 7. Cooling rate evolution in prototype 1 during tests with different inlet air temperatures (left); in prototype 2 (right).

leakages were found out. Some pouches were torn and the metallic grid was deformed by the pushing force of the solidification process of the PCM inside pouches. PCM thickness was twice higher than the designed, causing a higher and dominant heat transfer resistance by conduction inside the PCM. Therefore melting times were higher and flow rate had almost no influence. This prototype did not fulfil melting time requirements and was discarded.

Fig. 6. Experimental installation arrangement to test PCM-air heat exchangers prototypes.

Fig. 7. Cooling rate evolution in prototype 1 during tests with different inlet air

leakages were found out. Some pouches were torn and the metallic grid was deformed by the pushing force of the solidification process of the PCM inside pouches. PCM thickness was twice higher than the designed, causing a higher and dominant heat transfer resistance by conduction inside the PCM. Therefore melting times were higher and flow rate had almost no influence. This prototype did not fulfil melting time requirements and was

temperatures (left); in prototype 2 (right).

discarded.

Fig. 8. Cooling rate evolution with constant inlet air temperature in prototype 1 during tests with different air flow rates (left); in prototype 2 (right); Flow 1 is the mass flow, equals to 0.34 kg/s.

Prototype 2 was designed using aluminium panels filled with organic PCM. Configuration was also vertical. Air flows were parallel to the panels from top to bottom. Due to the fact that the PCM in this prototype was organic, it presented lower thermal conductivity than PCM in prototype 1. Furthermore, the total stored energy in prototype 2 was also lower than in prototype 1. The first tests were accomplished using constant inlet air temperature. As it can be seen in figure 7 (right), cooling rates were higher than those obtained with prototype 1 and the melting times were half the melting times with prototype 1.

Different air flow rates in prototype 2 were tested. It was observed that it had influence on the melting time and cooling power. Figure 8 (right) shows that for the lowest air flow rate, the heat rate curve changes its shape and is more similar to prototype 1. This indicates that heat transfer by conduction inside the PCM starts to be relevant when compared to air convection. The first results of prototype 2 were satisfactory, so more test were planned to evaluate its behaviour under real conditions. Two types of experiments were accomplished: constant rise of inlet air temperature and constant heating power. Temperature rise ramps were then set into the resistances controller: results showed that faster the rise, higher cooling power and lower melting time. For constant power tests, different heating powers of electrical resistances were fixed. Results showed that prototype 2 was able to maintain a cooling capacity over 3 kW for approximately 1h 30' or approximately 1 kW for more than 3 h. This result is useful to design the optimal operation mode depending on the application.

#### **3.3 Experimental results**

The total energy exchanged during melting and solidification, as well as the time elapsed until total melting/solidification were determined from the heat rate curves experimentally obtained. The influence of the inlet air temperature and air flow was studied, and results showed that the continuous thermal cycling of the unit is a repetitive process: running experiments with similar conditions led to the same thermal behaviour; no degradation in the PCM properties was noticed. Pressure drop was measured for different air flows. Depending on the inlet air temperature, full solidification of the PCM could be achieved in less than 3 h for an 8 ºC temperature difference between the inlet air and the average phase

PCM-Air Heat Exchangers: Slab Geometry 439

from all experimental outcomes under real conditions. The results to test prototype 2 with constant heating power were evaluated and an average temperature was obtained for the air temperature plateau ( Tplateau ). There was obtained a linear correlation (equation 5) between the average plateau temperatures and the heating power of the electrical resistances ( Qresistances ). The origin ordinate was the average phase change temperature of the PCM

The temperature on the surface of the PCM encapsulation was measured at 20 locations, distributed in such a way that the melting evolution can be studied and that any cold point can be detected. For each measurement, the average surface temperatures were obtained, as well as the plotting of the temperature difference between the inlet and the outlet (*ΔT*), the air temperature at the inlet of the TES unit (*Tinletair*) and the average surface temperatures ( surface T ). All locations come from experiments with the same airflow rate and heat exchange geometry. All measured locations were contained in a plane. A fitting tool was used to obtain the equation adjustment (equation 6). For the tested heat exchanger, expression A in figure 9 corresponds to equation 6. As it has been detailed, the heat exchanged with air was evaluated by an energy balance. The stored energy was then obtained for each measurement time step, as an accumulative result from the exchanged energy between air and PCM. The relationship between the stored energy and the average surface temperatures ( surface T ) corresponds to expression C in figure 9. Using the empirical model, different conditions were simulated, and the equations for the expressions shown in figure 9 were obtained, as well as various conclusions concerning PCM selection criteria

The guide EA-4/02 Expression of the Uncertainty of Measurement in Calibration, 1999, has been followed to estimate the uncertainty of measurements. Air flow was determined using an energy balance method that consists in applying an energy balance to the air flow that passes through the electrical resistances (equation 7). The air flow is measured with an accuracy of ± 0.026 kg/s (5.5% of measurement). Table 1 summarizes the expanded uncertainty estimation. The same procedure was followed to estimate the uncertainty of the cooling power determination. Equation 8 expresses the energy balance for the heat exchanger. Table 2 shows an example of uncertainty estimation in a cooling power measurement. In this case, the cooling power was measured with a ± 0.301 kW uncertainty

air m =Q c · air resistances p thermopile <sup>Δ</sup>T (7)

HX air p inlet-outlet air Q =m ·c · <sup>Δ</sup><sup>T</sup> (8)

T =26.6+1.58·Q plateau resistances (5)

surface air ΔT K =-1.4683-1.10943·T ºC +1.10706·T ºC inlet (6)

used.

(Lazaro et al., 2009b).

(9%).

**3.5 Uncertainties propagation** 

change of the PCM. Average heat rates of up to 4.5 kW and 3.5 kW for 1 h were obtained for melting and solidification stages, respectively (Dolado et al., 2011b; Dolado, 2011; Lazaro, 2009; Lazaro et al., 2009a).

#### **3.4 Empirical models**

From experimental results, the empirical models were built aimed at simulating the thermal behaviour in the tested heat exchanger in different cases. These simulations were used to evaluate the technical viability of its application. Since the thermal properties of PCM vary with temperature, a PCM-air heat exchanger works as a transitory system and therefore, its design must be based on transitory analysis. This section shows that PCM selection criteria must include the power demand. The conclusions obtained for the PCM-air heat exchange can be useful for selecting PCM for other heat exchanger applications that use the tested geometry as well as for applications that use such technology: green housing, curing and drying processes, industrial plant production, HVAC, free-cooling.

Fig. 9. Flow diagram to evaluate the technical feasibility of the system (Lazaro et al., 2009b).

The technical viability of TES with PCM systems depends on the capability to maintain the temperature below a maximum level (*Tob*), during a specific period of time (*Δtob*). Figure 9 shows a flow diagram of the relevant parameters to test the technical viability of a system built with a specific geometry and an established air flow, for the heat exchange between air and PCM. To simulate a room under various conditions, the internal, external, and ventilation loads of the room, as well as the transitory response of the TES system using PCM must be known. Such transitory response was obtained by an empirical model built

change of the PCM. Average heat rates of up to 4.5 kW and 3.5 kW for 1 h were obtained for melting and solidification stages, respectively (Dolado et al., 2011b; Dolado, 2011; Lazaro,

From experimental results, the empirical models were built aimed at simulating the thermal behaviour in the tested heat exchanger in different cases. These simulations were used to evaluate the technical viability of its application. Since the thermal properties of PCM vary with temperature, a PCM-air heat exchanger works as a transitory system and therefore, its design must be based on transitory analysis. This section shows that PCM selection criteria must include the power demand. The conclusions obtained for the PCM-air heat exchange can be useful for selecting PCM for other heat exchanger applications that use the tested geometry as well as for applications that use such technology: green housing, curing and

Fig. 9. Flow diagram to evaluate the technical feasibility of the system (Lazaro et al., 2009b).

The technical viability of TES with PCM systems depends on the capability to maintain the temperature below a maximum level (*Tob*), during a specific period of time (*Δtob*). Figure 9 shows a flow diagram of the relevant parameters to test the technical viability of a system built with a specific geometry and an established air flow, for the heat exchange between air and PCM. To simulate a room under various conditions, the internal, external, and ventilation loads of the room, as well as the transitory response of the TES system using PCM must be known. Such transitory response was obtained by an empirical model built

drying processes, industrial plant production, HVAC, free-cooling.

2009; Lazaro et al., 2009a).

**3.4 Empirical models** 

from all experimental outcomes under real conditions. The results to test prototype 2 with constant heating power were evaluated and an average temperature was obtained for the air temperature plateau ( Tplateau ). There was obtained a linear correlation (equation 5) between the average plateau temperatures and the heating power of the electrical resistances ( Qresistances ). The origin ordinate was the average phase change temperature of the PCM used.

$$
\overline{T}\_{\text{plateau}} = 26.6 + 1.58 \text{ \dot{Q}}\_{\text{resistance}} \tag{5}
$$

$$
\Delta\text{T[K]=-1.4683-1.10943 }\overline{\text{T}}^{\text{surface}}\left[\text{\textquotedbl{}C}\right] + 1.10706 \cdot \text{T}^{\text{air}}\_{\text{irlet}}\left[\text{\textquotedbl{}C}\right] \tag{6}
$$

The temperature on the surface of the PCM encapsulation was measured at 20 locations, distributed in such a way that the melting evolution can be studied and that any cold point can be detected. For each measurement, the average surface temperatures were obtained, as well as the plotting of the temperature difference between the inlet and the outlet (*ΔT*), the air temperature at the inlet of the TES unit (*Tinletair*) and the average surface temperatures ( surface T ). All locations come from experiments with the same airflow rate and heat exchange geometry. All measured locations were contained in a plane. A fitting tool was used to obtain the equation adjustment (equation 6). For the tested heat exchanger, expression A in figure 9 corresponds to equation 6. As it has been detailed, the heat exchanged with air was evaluated by an energy balance. The stored energy was then obtained for each measurement time step, as an accumulative result from the exchanged energy between air and PCM. The relationship between the stored energy and the average surface temperatures ( surface T ) corresponds to expression C in figure 9. Using the empirical model, different conditions were simulated, and the equations for the expressions shown in figure 9 were obtained, as well as various conclusions concerning PCM selection criteria (Lazaro et al., 2009b).

#### **3.5 Uncertainties propagation**

The guide EA-4/02 Expression of the Uncertainty of Measurement in Calibration, 1999, has been followed to estimate the uncertainty of measurements. Air flow was determined using an energy balance method that consists in applying an energy balance to the air flow that passes through the electrical resistances (equation 7). The air flow is measured with an accuracy of ± 0.026 kg/s (5.5% of measurement). Table 1 summarizes the expanded uncertainty estimation. The same procedure was followed to estimate the uncertainty of the cooling power determination. Equation 8 expresses the energy balance for the heat exchanger. Table 2 shows an example of uncertainty estimation in a cooling power measurement. In this case, the cooling power was measured with a ± 0.301 kW uncertainty (9%).

$$\dot{\mathbf{m}}\_{\text{air}} = \dot{\mathbf{Q}}\_{\text{resistance}} \sqrt{\left(\mathbf{c}\_{p\_{\text{air}}} \cdot \Delta \mathbf{T}\_{\text{thermropic}}\right)} \tag{7}$$

$$
\dot{\mathbf{Q}}\_{\rm HX} = \dot{\mathbf{m}}\_{\rm air} \mathbf{v}\_{\mathbf{p}\_{air}} \Delta T\_{\rm inlet-outlet} \tag{8}
$$

PCM-Air Heat Exchangers: Slab Geometry 441

are rigorous and flexible to simulate heat exchangers of air and PCMs. When developing a model, the trade-off between rigour and computational cost is crucial. There are many options reported in scientific literature to face the mathematical problem of phase change as well as to solve specific particularities such as hysteresis phenomena (Bony & Citherlet, 2007) or sub-cooling (Günther et al., 2007). In the review by Zalba et al., 2003, the authors presented a comprehensive compilation of TES with PCM. The authors remarked that although there is a huge amount of published articles dealing with the heat transfer analysis of the phase change, the modelling of latent heat TES systems still remains a challenging task. When working with commercially available PCM (or mixtures or impure materials), the phase change takes place over a temperature range and therefore a two-phase zone appears between the solid and liquid phases. In these cases, it is appropriate to consider the energy equation in terms of enthalpy (Zukowsky, 2007b). When the advective movements

> 

The solution of this equation requires knowledge of the h-T functional dependency and the λ-T curve. The advantage of this methodology is that the equation is applicable to every phase; the temperature is determined at each point and the value of the thermo-physical properties can then be evaluated. In thermal simulations of PCM, the accuracy of the results relies on the material properties' data (Arkar & Medved, 2005). In the geometry studied in this work, the main properties were enthalpy and thermal conductivity, but notice that the rate of melting/solidification can also depend on other material properties such as viscosity or density (Hamdan & Elwerr, 1996). In the models developed here the variation of

A PCM plate model was developed with finite differences, one-dimensional, implicit formulation. Implicit formulation was selected because of its unconditional stability. The basis model assumed only conduction heat transfer inside the PCM plate, in a normal direction to the air flow. The model analyzed the temperature of the airflow in a onedimensional way. Due to its symmetry, the analyzed system was a division of the prototype. In the present work, the model was implemented in Matlab R2008b. The software implements direct methods, variants of Gaussian elimination, through the matrix division

The study system is the PCM-air TES unit (figure 5, right). The air inlet was located on the upper side of the TES unit. The air flowed downwards in the TES unit, circulating parallel to the PCM slabs, exchanging energy with the PCM, and eventually was blown outside the TES unit by a centrifugal fan. The system was studied from the point of view of a single slab. As the PCM zone of the TES unit was insulated as well as due to the distribution of the slabs inside the TES unit, some symmetry relationships could be considered so only the dotted domain in figure 5 (left) was modelled. The nodal distribution of the mathematical model is shown in figure 10. Depending on whether the encapsulation is considered or not, two more nodes have to be included between the PCM surface and the airflow. In the experimental study, the heat transfer processes that take place inside the TES unit between the air flowing through the slabs and the PCM inside the slabs were: forced convection in

· · *ht T* (9)

within the liquid are negligible, the energy equation is expressed as follows:

thermophysical properties with temperature in all phases was considered.

operators, which can be used to solve linear systems.

**4.2 Development of a 1D finite differences equations model for PCM plates** 


Table 1. Air flow uncertainty of measurement estimation.


Table 2. Example of cooling power uncertainty estimation.

The air temperature difference between the inlet and the outlet of the heat exchanger was measured using a thermopile and two Pt-100 in the centre of the air ducts. Measurements were compared during stationary periods in order to confirm the fact that a thermopile was more appropriate. Standard deviations of thermopile measurements are lower than the ones for Pt-100 differences and the mean values are all higher for the Pt-100 differences. This is due to the fact that Pt-100 are located at a specific point in the centre of the air duct whereas the thermopiles are distributed in the air duct cross surface.

## **4. Study of the heat transfer**

In this section a theoretical model has been developed to perform the computer simulation of the thermal behaviour of a PCM-air heat exchanger, validating the theoretical model with the results obtained from the prototype in the experimental facility built for this purpose. In the archival literature, the approach of the solid-liquid phase change problem appears with different configurations, this section is focused on the case of macroencapsulated PCM, plate shape. Among the different numerical methods for solving the problem, in this section the energy equation is considered in terms of enthalpy as the governing equation can be applied at any stage, the temperature can be determined at each point, and therefore the values of the thermophysical properties can be evaluated. The PCM simulated are commercially available so the simulation involves, among other problems the non-linearity. The finite difference method for discretization of the governing equations is used. The models are based on 1D conduction analysis, using the thermo-physical data of the PCM measured in the laboratory. The models can take into account the hysteresis of the enthalpy curve and the convection inside the PCM, using effective conductivity when necessary.

#### **4.1 Modelling the solid-liquid phase change**

Modelling is a useful tool in a viability analysis of applications that involve TES by solid– liquid PCM. Therefore, there is a necessity to develop experimentally validated models that

ΔTthermopile[K] 0.51 0.255 9.3 1 0.0008

ΔTthermopile[K] 0.51 0.255 7 1 0.0013 mair [kg/s] 0.026 0.013 0.47 1 0.0008

> Estimated value

The air temperature difference between the inlet and the outlet of the heat exchanger was measured using a thermopile and two Pt-100 in the centre of the air ducts. Measurements were compared during stationary periods in order to confirm the fact that a thermopile was more appropriate. Standard deviations of thermopile measurements are lower than the ones for Pt-100 differences and the mean values are all higher for the Pt-100 differences. This is due to the fact that Pt-100 are located at a specific point in the centre of the air duct whereas

In this section a theoretical model has been developed to perform the computer simulation of the thermal behaviour of a PCM-air heat exchanger, validating the theoretical model with the results obtained from the prototype in the experimental facility built for this purpose. In the archival literature, the approach of the solid-liquid phase change problem appears with different configurations, this section is focused on the case of macroencapsulated PCM, plate shape. Among the different numerical methods for solving the problem, in this section the energy equation is considered in terms of enthalpy as the governing equation can be applied at any stage, the temperature can be determined at each point, and therefore the values of the thermophysical properties can be evaluated. The PCM simulated are commercially available so the simulation involves, among other problems the non-linearity. The finite difference method for discretization of the governing equations is used. The models are based on 1D conduction analysis, using the thermo-physical data of the PCM measured in the laboratory. The models can take into account the hysteresis of the enthalpy curve and

Modelling is a useful tool in a viability analysis of applications that involve TES by solid– liquid PCM. Therefore, there is a necessity to develop experimentally validated models that

the convection inside the PCM, using effective conductivity when necessary.

**4.1 Modelling the solid-liquid phase change** 

Estimated value

0.044 0.022 4.4 1 0.0000

Standard uncertainty

Estimated value

Sensibility coefficient

Sensibility coefficient

Expanded uncertainty

Contribution to uncertainty

Contribution to uncertainty

Standard uncertainty

Standard uncertainty

QHX [kW] 0.089 3.29 0.151 0.301

Expanded uncertainty

Table 1. Air flow uncertainty of measurement estimation.

Sum of contributions

Table 2. Example of cooling power uncertainty estimation.

the thermopiles are distributed in the air duct cross surface.

**4. Study of the heat transfer** 

Expanded uncertainty

Electrical power consumption of resistances [kW] are rigorous and flexible to simulate heat exchangers of air and PCMs. When developing a model, the trade-off between rigour and computational cost is crucial. There are many options reported in scientific literature to face the mathematical problem of phase change as well as to solve specific particularities such as hysteresis phenomena (Bony & Citherlet, 2007) or sub-cooling (Günther et al., 2007). In the review by Zalba et al., 2003, the authors presented a comprehensive compilation of TES with PCM. The authors remarked that although there is a huge amount of published articles dealing with the heat transfer analysis of the phase change, the modelling of latent heat TES systems still remains a challenging task. When working with commercially available PCM (or mixtures or impure materials), the phase change takes place over a temperature range and therefore a two-phase zone appears between the solid and liquid phases. In these cases, it is appropriate to consider the energy equation in terms of enthalpy (Zukowsky, 2007b). When the advective movements within the liquid are negligible, the energy equation is expressed as follows:

$$
\rho \mathbf{v} \,\hat{\imath} \mathbf{t} / \hat{\imath} \mathbf{t} = \nabla \left( \mathbb{X} \nabla T \right) \tag{9}
$$

The solution of this equation requires knowledge of the h-T functional dependency and the λ-T curve. The advantage of this methodology is that the equation is applicable to every phase; the temperature is determined at each point and the value of the thermo-physical properties can then be evaluated. In thermal simulations of PCM, the accuracy of the results relies on the material properties' data (Arkar & Medved, 2005). In the geometry studied in this work, the main properties were enthalpy and thermal conductivity, but notice that the rate of melting/solidification can also depend on other material properties such as viscosity or density (Hamdan & Elwerr, 1996). In the models developed here the variation of thermophysical properties with temperature in all phases was considered.

#### **4.2 Development of a 1D finite differences equations model for PCM plates**

A PCM plate model was developed with finite differences, one-dimensional, implicit formulation. Implicit formulation was selected because of its unconditional stability. The basis model assumed only conduction heat transfer inside the PCM plate, in a normal direction to the air flow. The model analyzed the temperature of the airflow in a onedimensional way. Due to its symmetry, the analyzed system was a division of the prototype. In the present work, the model was implemented in Matlab R2008b. The software implements direct methods, variants of Gaussian elimination, through the matrix division operators, which can be used to solve linear systems.

The study system is the PCM-air TES unit (figure 5, right). The air inlet was located on the upper side of the TES unit. The air flowed downwards in the TES unit, circulating parallel to the PCM slabs, exchanging energy with the PCM, and eventually was blown outside the TES unit by a centrifugal fan. The system was studied from the point of view of a single slab. As the PCM zone of the TES unit was insulated as well as due to the distribution of the slabs inside the TES unit, some symmetry relationships could be considered so only the dotted domain in figure 5 (left) was modelled. The nodal distribution of the mathematical model is shown in figure 10. Depending on whether the encapsulation is considered or not, two more nodes have to be included between the PCM surface and the airflow. In the experimental study, the heat transfer processes that take place inside the TES unit between the air flowing through the slabs and the PCM inside the slabs were: forced convection in

PCM-Air Heat Exchangers: Slab Geometry 443

The dominant resistance of the process could be convection on the air side and not always conduction–convection in the PCM. In this case the thermal resistance of the encapsulation was very low, and therefore it was not necessary to consider encapsulation in the node system, as the heat transfer process was controlled by the convection on the air side and/or by the conduction–natural convection in the PCM. However, in other cases it is not always possible to disregard the thermal influence of the encapsulation, and therefore two models were developed: the first model did not take into account the thermal behaviour of the encapsulation and the second model did, and it was developed in order to be used for general purposes. The node equations of the two models are summarized in tables 3 and 4. Important aspects to consider when dealing with the simulation of this type of heat exchanger are as follows: friction factor (rugosity of the encapsulation surface), convection coefficient, thermophysical properties of the PCM (as functions of temperature), hysteresis, natural convection inside the PCM, thermal losses/gains through the TES casing, etc. More

detailed information can be found in Dolado et al., 2011a, and Dolado, 2011.

variables; 7) To determine the uncertainty of the theoretical model results.

Air conditions at the inlet of the TES unit (temperature and airflow);

were considered and classified into three groups:

conductivity curve, λ-T);

For the current study, the following parameters that introduce uncertainty in the results,

Material properties (parameterized enthalpy-temperature curve, h-T, and thermal

 Geometric parameters (PCM plate thickness and width of the air gap between plates). The parameterization and the range of uncertainty assumed for all these parameters are detailed in Dolado, 2011. The confidence level is 97.5%. Furthermore, as the probability distribution of the different parameters is unknown, a normal distribution is taken in all

**model** 

**4.3 Experimental validation: applying the uncertainties propagation approach to the** 

The validation stage of a theoretical model has become a fundamental objective to evaluate the precision, accuracy and reliability of computer simulations used in design. Uncertainties can be associated with the own theoretical model, as well as with the measurement systems used to characterize the process of interest or even with the manufacturing process of the equipment. Therefore, assessing the validity of an approximation of a theoretical model must be carried out based on stochastic measurements to ensure the trust of designers in the use of the model. This improved knowledge of the theoretical model helps to know what are the most critical factors as model inputs and, therefore, indicates what should be more controlled in its determination or measurement. It also allows establishing an uncertainty band set around the solution bringing more rigor to the model simulations. The methodology of uncertainties propagation is an external method used to analyze the system by means of the input-output analysis, instead of the traditional equation of uncertainties propagation applied to a known function. The whole methodology followed in is summarized in the next steps: 1) To select the variables under study; 2) To allocate the probability distributions of each variable; 3) To generate samples for the different runs of the theoretical model (by means of Latin hypercube sampling); 4) To run the program once per sample; 5) To analyze the relations between the inputs and the outputs; 6) To classify the

the air, conduction in the shell of the aluminium slab, and conduction and natural convection in the PCM itself.

Fig. 10. Nodal distribution in the 1D plate system.


Table 3. Node temperature equations not considering encapsulation.


Table 4. Node temperature equations considering encapsulation.

the air, conduction in the shell of the aluminium slab, and conduction and natural

convection in the PCM itself.

Surface Encapsulation

Encapsulation

Fig. 10. Nodal distribution in the 1D plate system.

Nodes Equations Air Flow T =T -NTU T -T air air-1 air air-1 surface

PCM Surface t-1 T = T +FoT +FoBiT 1+Fo+FoBi surface PCM PCM+1 air

 

enc-PCM PCM-enc

enc enc enc

surface enc enc enc enc air

T +2 Fo T +Fo Bi T

1+2 Fo Bi +Fo 

enc enc-PCM surface PCM-enc PCM

 

T +2 Fo T +2Fo T

1+2 Fo +2Fo

PCM Inner t-1 T = T +Fo T +T 1+2Fo PCM PCM PCM-1 PCM+1

PCM Central t-1 T = T +2FoT 1+2Fo PCM PCM PCM-1

Air Flow T =T -NTU T -T air air-1 air air-1 surface

Air Surface

t-1

PCM t-1 T = T +Fo 2T +T 1+3Fo PCM PCM enc PCM+1

PCM Inner t-1 T = T +Fo T +T 1+2Fo PCM PCM PCM-1 PCM+1

PCM Central t-1 T = T +2FoT 1+2Fo PCM PCM PCM-1

t-1

Table 3. Node temperature equations not considering encapsulation.

Nodes Equations

surface

enc

Table 4. Node temperature equations considering encapsulation.

T =

T =

The dominant resistance of the process could be convection on the air side and not always conduction–convection in the PCM. In this case the thermal resistance of the encapsulation was very low, and therefore it was not necessary to consider encapsulation in the node system, as the heat transfer process was controlled by the convection on the air side and/or by the conduction–natural convection in the PCM. However, in other cases it is not always possible to disregard the thermal influence of the encapsulation, and therefore two models were developed: the first model did not take into account the thermal behaviour of the encapsulation and the second model did, and it was developed in order to be used for general purposes. The node equations of the two models are summarized in tables 3 and 4. Important aspects to consider when dealing with the simulation of this type of heat exchanger are as follows: friction factor (rugosity of the encapsulation surface), convection coefficient, thermophysical properties of the PCM (as functions of temperature), hysteresis, natural convection inside the PCM, thermal losses/gains through the TES casing, etc. More detailed information can be found in Dolado et al., 2011a, and Dolado, 2011.

#### **4.3 Experimental validation: applying the uncertainties propagation approach to the model**

The validation stage of a theoretical model has become a fundamental objective to evaluate the precision, accuracy and reliability of computer simulations used in design. Uncertainties can be associated with the own theoretical model, as well as with the measurement systems used to characterize the process of interest or even with the manufacturing process of the equipment. Therefore, assessing the validity of an approximation of a theoretical model must be carried out based on stochastic measurements to ensure the trust of designers in the use of the model. This improved knowledge of the theoretical model helps to know what are the most critical factors as model inputs and, therefore, indicates what should be more controlled in its determination or measurement. It also allows establishing an uncertainty band set around the solution bringing more rigor to the model simulations. The methodology of uncertainties propagation is an external method used to analyze the system by means of the input-output analysis, instead of the traditional equation of uncertainties propagation applied to a known function. The whole methodology followed in is summarized in the next steps: 1) To select the variables under study; 2) To allocate the probability distributions of each variable; 3) To generate samples for the different runs of the theoretical model (by means of Latin hypercube sampling); 4) To run the program once per sample; 5) To analyze the relations between the inputs and the outputs; 6) To classify the variables; 7) To determine the uncertainty of the theoretical model results.

For the current study, the following parameters that introduce uncertainty in the results, were considered and classified into three groups:


The parameterization and the range of uncertainty assumed for all these parameters are detailed in Dolado, 2011. The confidence level is 97.5%. Furthermore, as the probability distribution of the different parameters is unknown, a normal distribution is taken in all

PCM-Air Heat Exchangers: Slab Geometry 445

evolution of the exchanged heat rate in the melting process with the associated uncertainty interval (97.5%) and also the results of the relative error in heat rate. It is observed that the relative error is below 10 % until the process is approaching the end of the melting stage. Is then when the absolute values of heat rate are smaller and the relative error increases until the melting ends, as expected. This result is analogous in the solidification stage. As stated by Dolado et al., 2011b, and Lazaro et al., 2009a, considering the instrumentation used in the experimental setup, an uncertainty of 9 % in terms of the heat rate during the first hour of a typical test process is obtained. The calculations to determine the uncertainty in the measurement of the heat rate are shown in table 7. The uncertainty is estimated for each measurement so that, according to the EA Guide 4/02 Expression of the Uncertainty of Measurement in Calibration, 1999, a band of uncertainty associated with the experimental

0 1800 3600 5400 7200 9000 10800

**time (s)**

Heat transfer rate, reference case Upper error curve Lower error curve Relative error Fig. 11. Simulated heat rate, uncertainty bands and relative error (melting stage).

> Standard uncertainty

Sum of contributions Estimated

7+(0.255/ΔTi)2] Q*<sup>i</sup>*

ΔTthermopile[K] 0.51 0.255 ΔTi 1 (0.255/ ΔTi)2 m [kg/s] 0.026 0.013 0.36 1 0.0013 cp [J/(kg·K)] 2 1 1007 1 9.86·10-7

As expected, the relative errors grow as the absolute value of the heat rate decreases (figures 11 and 12). This is because the expanded uncertainty associated with the measure of the thermopile is ± 0.51 °C (Lazaro, 2009) so that the error increases as the temperature difference between the air at the inlet and at the outlet of the TES unit decreases. Figure 12

Estimated value

value

0% 10% 20% 30% 40% 50% 60% 70% 80%

Sensibility coefficient

Standard uncertainty

w ·Q <sup>Q</sup> <sup>i</sup>

Contribution to uncertainty

Expanded uncertainty

<sup>Q</sup> <sup>i</sup> 2·w ·Q

**Relative error**

heat rate curve can be determined (figure 12).

0

Expanded uncertainty

<sup>Q</sup> [W] SQRT [(0.0013+9.86·10-

Table 7. Uncertainties determination of the experimental heat rate.

1000

2000

3000

**Heat transfer rate (W)**

4000

5000

cases. The study of uncertainties propagation is performed by numerical simulation of sets of input values of those parameters. Traditionally the random sampling technique has been used, which was followed by an improved version as the stratified sampling and subsequently the Latin hypercube (McKay et al., 1979) which is a valid technique used here. The variation range of the studied variable was analyzed. This result provides an estimate of the uncertainty range that would result as the TES unit is designed. This analysis sets the interval for the output variable of interest. For every particular time of the simulation the error was calculated as: 1) Histogram to the corresponding time of the results set of all simulations (the difference between the value in the current simulation and the reference case); 2) Function of the cumulative probability distribution; 3) Error in the given time at a certain confidence level; 4) Graphical representation of the reference case and of the error.


Table 5. Parameters setting of the enthalpy-temperature curve.


Table 6. Parameters setting of the inlet air.

The corresponding value for a cumulative relative frequency of 97.5% was taken as the uncertainty range in this approach. Running the theoretical model with the reference case conditions (tables 5 and 6) yielded to the results shown in figure 11. The figure shows that the heat rate has an initial plateau of about 4500 W, with a duration of 40 minutes, which reduces while reaching the complete melting of the PCM. The full melting takes place two hours after the start of the process. From these graphical results, the responses of interest were obtained. Among the responses of interest provided by the theoretical model, the analysis was focused on the following responses: average heat rate in the first hour ( Q *average,1h*), and time until the air reaches 32 ºC at the outlet (*tTair,out=32ºC*), as both responses pose a greater interest in practical operation of the TES unit. Immediately afterwards, applying the Latin hypercube sampling to run the numerical simulations, instead of getting a single result of Q *average,1h*, a complete distribution was obtained. Figure 11 shows the

cases. The study of uncertainties propagation is performed by numerical simulation of sets of input values of those parameters. Traditionally the random sampling technique has been used, which was followed by an improved version as the stratified sampling and subsequently the Latin hypercube (McKay et al., 1979) which is a valid technique used here. The variation range of the studied variable was analyzed. This result provides an estimate of the uncertainty range that would result as the TES unit is designed. This analysis sets the interval for the output variable of interest. For every particular time of the simulation the error was calculated as: 1) Histogram to the corresponding time of the results set of all simulations (the difference between the value in the current simulation and the reference case); 2) Function of the cumulative probability distribution; 3) Error in the given time at a certain confidence level; 4) Graphical representation of the reference case and of the error.

Parameter b ΔT hl Tsl

Thermal

value 3 kJ/(kg·K) 0.9 ºC 170 kJ/kg 27 ºC Uncertainty ± 15 % ± 0.2 ºC ± 20 kJ/kg ± 1 ºC

**Y**

yo

yt

window Latent heat

**X**

**Y**

yo

yt

**X**

xo

xo

Variable Reference value Uncertainty *V* 1400 m3/h ± 86 m3/h *T* experimental curve of Tair,in, ºC ± 0.6 ºC

The corresponding value for a cumulative relative frequency of 97.5% was taken as the uncertainty range in this approach. Running the theoretical model with the reference case conditions (tables 5 and 6) yielded to the results shown in figure 11. The figure shows that the heat rate has an initial plateau of about 4500 W, with a duration of 40 minutes, which reduces while reaching the complete melting of the PCM. The full melting takes place two hours after the start of the process. From these graphical results, the responses of interest were obtained. Among the responses of interest provided by the theoretical model, the analysis was focused on the following responses: average heat rate in the first hour ( Q *average,1h*), and time until the air reaches 32 ºC at the outlet (*tTair,out=32ºC*), as both responses pose a greater interest in practical operation of the TES unit. Immediately afterwards, applying the Latin hypercube sampling to run the numerical simulations, instead of getting a single result of Q *average,1h*, a complete distribution was obtained. Figure 11 shows the

Average temperature of phase change

**Y**

yo

yt

**X**

xo

Related with

Reference

Effect on h-T curve

Slope of the sensible heat

**X**

Table 5. Parameters setting of the enthalpy-temperature curve.

xo

**Y**

yo

Table 6. Parameters setting of the inlet air.

yt

evolution of the exchanged heat rate in the melting process with the associated uncertainty interval (97.5%) and also the results of the relative error in heat rate. It is observed that the relative error is below 10 % until the process is approaching the end of the melting stage. Is then when the absolute values of heat rate are smaller and the relative error increases until the melting ends, as expected. This result is analogous in the solidification stage. As stated by Dolado et al., 2011b, and Lazaro et al., 2009a, considering the instrumentation used in the experimental setup, an uncertainty of 9 % in terms of the heat rate during the first hour of a typical test process is obtained. The calculations to determine the uncertainty in the measurement of the heat rate are shown in table 7. The uncertainty is estimated for each measurement so that, according to the EA Guide 4/02 Expression of the Uncertainty of Measurement in Calibration, 1999, a band of uncertainty associated with the experimental heat rate curve can be determined (figure 12).

Fig. 11. Simulated heat rate, uncertainty bands and relative error (melting stage).


Table 7. Uncertainties determination of the experimental heat rate.

As expected, the relative errors grow as the absolute value of the heat rate decreases (figures 11 and 12). This is because the expanded uncertainty associated with the measure of the thermopile is ± 0.51 °C (Lazaro, 2009) so that the error increases as the temperature difference between the air at the inlet and at the outlet of the TES unit decreases. Figure 12

PCM-Air Heat Exchangers: Slab Geometry 447

not a design tool, this methodology is proposed to size the equipment. This technique greatly reduces the time spent in performing the simulations required to find the optimal equipment (Del Coz Díaz et al., 2010) as well as and a potential cost saving on the experimental (Del Coz Díaz et al., 2010; Gunasegaram et al., 2009) if the prototype-model similarity relations are met. Moreover, contrary to a sequential analysis, it is reasonable to use a mathematical and statistical methodology that allows planning the sequence of

An empirical model was built from the experimental results described in the previous section. The aim was to simulate the thermal behaviour of the tested heat exchanger in different cases. These simulations were used to evaluate the technical viability of application. The model describes the temperature evolution of a room with an internal cooling demand ( Q *demand*), where the PCM-air heat exchanger is operating and there is a ventilation system. The enclosure temperature was considered to be the average between the outside temperature and the room temperature. A diagram of the room is shown in figure 13. Expression D in figure 9 is equivalent to equation 9, expressing the energy balance

experiments on the philosophy of maximum information with minimum effort.

**5.1 Empirical model: simulations of a case study and modular design** 

Fig. 13. Schematic diagram of the room in which the temperature is evaluated.

 

ventilation p,air outside room demand air HX p,air

i i-1

T +T T +T

); and the surface T at instant *i* is obtained from the stored

i i-1

room outside room outside i i-1

p,enclosure air p,air room room

where *ΔT* is obtained at each instant as a function of surface T and the inlet air temperature,

The real-scale PCM-air heat exchanger tested was constituted of 18 parallel modules (*#modules* denotes de number of PCM modules in the heat exchanger). A module is constituted by a

m·c · - =ρ ·V·c · T -T 2 2

= m ·c · T -T +Q -m ·c ·ΔT ·Δt i,i-1

(9)

applied to the air inside the room.

*Tinletair* (at instant *i* equal to *Troomi*

energy evolution.

shows the overlap between the experimental curve (including the lower and upper limits associated with its uncertainty) and the simulation (including the uncertainty of the response heat rate calculated applying the reported technique). The agreement is significant in most of the process, finding the more relevant discrepancies as the curve reaches the end of the corresponding stage of the cycle (i.e. as the heat rate values are smaller).

Fig. 12. Comparison of experimental and simulated results (including their corresponding uncertainty bands) for the melting (up) and solidification (down) of a thermal cycle.
