Prediction of yearly energy requirements of indoor ice rinks

Energy and Buildings 41 (2009) 500–511 Contents lists available at ScienceDirect Energy and Buildings journal homepage: www.elsevier.com/locate/enbuild Prediction of yearly energy requirements of indoor ice rinks Lotfi Seghouani, Ahmed Daoud, Nicolas Galanis * Département de Génie Mécanique, Université de Sherbrooke, Sherbrooke, Qc, Canada J1K 2R1 A R T I C L E I N F O A B S T R A C T Article history: Received 24 April 2008 Received in revised form 20 October 2008 Accepted 22 November 2008 A model of the transient heat transfer between the ground under and around the foundations of an indoor ice rink and the brine circulating in pipes embedded in the concrete slab under the ice has been coupled with a previously developed model calculating heat fluxes towards the ice by convection, radiation and phase changes. Subroutines calculating the energy consumption for heating and humidifying (or cooling and reheating) the ventilation air have also been added to the model. The resulting simulation tool has been used to calculate monthly refrigeration loads and energy consumption by the ventilation system, the lights, the brine pump, the radiant heating system of the stands and the underground electric heating used to prevent freezing and heaving for four North American cities with very different climates. Correlations expressing the energy consumption of the ventilation air stream in terms of the sol-air temperature are formulated. ß 2008 Elsevier B.V. All rights reserved. Keywords: Zonal method Ground conduction Radiation exchanges Convection Condensation Refrigeration load 1. Introduction Indoor ice rinks are large buildings without internal partitions and with high-energy consumption. They have a complex energy system in which a large ice sheet is cooled and maintained at a low temperature by a refrigeration system, while the stands are heated (or cooled) to ensure comfortable conditions for the spectators. Also, the building is ventilated to ensure good air quality. The movement of the ventilation air through these wide-open areas and the simultaneous operation of heating and cooling equipments increase energy consumption and greenhouse gas (GHG) emissions. A study by Lavoie et al. [1] shows that the potential for energy savings in a typical ice rink in Quebec is roughly 620 MWh/year and the potential GHG emission reduction is 146 tons-equivalent CO2/year. Since there are 435 indoor ice rinks in Quebec and several thousand in North America, it would be interesting to improve their energy efficiency while preserving good ice quality and comfort for the spectators. To achieve this objective precise methods for the calculation of the corresponding loads are necessary. Three different methods are commonly used for the thermal modeling of buildings: the nodal method, computational fluid dynamics (CFD) and the zonal method. The first one is the simplest * Corresponding author. E-mail address: [email protected] (N. Galanis). Abbreviations: AIM, above ice model; BIM, below ice model; IST, ice surface temperature. 0378-7788/$ – see front matter ß 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.enbuild.2008.11.014 and is implemented by representing the inside volume of the entire building, or of large parts thereof, by a single node. Therefore the nodal method does not necessitate an important computing capacity but, on the other hand, it does not provide a detailed description of the indoor conditions. The application of such a model to ice rinks (or other large buildings without internal partitions such as supermarkets or gymnasia) can lead to very imprecise results because the mass fluxes between different parts of the inside volume are extremely difficult to estimate. On the other hand, the modeling of an ice rink for CFD calculations is very complex due firstly to their size and geometry, and secondly to the variety of the heat and mass transfer mechanisms which take place therein. Thus, the model must take into account heat transfer through the envelope and heat gains from the ground, air motion within the building due to forced and natural convection, vapour diffusion and condensation on the ice sheet, heat transfer by radiation between all internal surfaces, conduction in the ice and floor as well as heat generation by the lights, the resurfacing operations, the refrigeration system, etc. Hence, the literature review revealed few CFD studies for large buildings such as ice rinks. Jones and Whittle [2] described the status and capabilities of CFD for building air flow prediction while Jian and Chen [3] as well as Yang et al. [4] used a CFD code to evaluate air quality in large ventilated enclosures. However, these studies did not calculate heating and refrigeration loads and ignored the interaction between the indoor and outdoor environments. More recently Bellache et al. [5,6] have carried out numerical simulations in 2D and steady state conditions using a CFD code which predicts velocity, temperature and absolute humidity distributions in an indoor ice rink with ventilation and heating. L. Seghouani et al. / Energy and Buildings 41 (2009) 500–511 Nomenclature A Cp Cd g h k ṁi; j ṁB M P qcd qcv qrd qcond qrs QCool QHeat QHumid QRe heat QH QIce Qnode R Snode t T Tb Tgr Tnode Tsol-air Unode Ẇg z1, z2 area (m2) specific heat (J/kg K) discharge coefficient acceleration of gravity (m/s2) height (m) thermal conductivity (W/(m K) airflow between zones i and j (kg/s) brine flow rate (kg/s) mass (kg) static pressure (Pa) conductive flux (W/m2) convective flux (W/m2) radiative flux (W/m2) condensation flux (W/m2) heat flux due to resurfacing (W/m2) cooling rate (W) heating rate (W) energy rate due to humidification (W) energy rate due to reheating (W) electrical power in sand layer (W) heat rate into node 1 calculated by AIM (W) lateral heat transfer to node n (W) thermal resistance (m2 K/W) surface for lateral heat transfer at node n (m2) time (h) temperature (K) brine temperature (K) temperature of ground surface (K) temperature at node n (K) sol-air temperature (8C) average conductance for lateral heat transfer at node n (W/m2 K) moisture source term (kg/s) top and bottom depths of ground segment (m) Subscripts and superscripts i i, j o p p+1 cell i or surface i between surface or cell i and j outside present time step next time step Greek symbols Dt time step (s) eij constant depending on flow direction (1) r air density (kg/m3) v absolute humidity (kgmoisture/kgdry air) The CFD code also calculates the heat fluxes toward the ice due to convection from the air, to condensation of vapour and to radiation from the walls and ceiling. However, these calculations did not take into account the contributions of ice resurfacing, system pump work and ground heat to the refrigeration load. This 2D CFD model was later improved by Bellache et al. [7] by including transient phenomena, heat transfer through the ground and energy gains from lights as well as the effects of resurfacing and dissipation of pump work in the coolant pipes. The ground at a 501 depth of 2 m was assumed to have a constant temperature while the horizontal plane through the centers of the brine pipes was assumed to be an isothermal surface with temperature equal to the average of the supply and return brine temperatures. Ouzzane et al. [8] contributed preliminary experimental measurements for a Canadian indoor ice rink which provide a better understanding of its thermal and energy behaviour. These measured values were also used for the verification and calibration of the numerical model developed by Bellache et al. [7]. The main drawback of the CFD approach is that it requires considerable computer memory and CPU time for the simulations. Thus, the transient 2D model by Bellache et al. [7] requires approximately 24 h of calculations on a modern desktop computer to simulate the response of an ice rink over a period of 1 day. An alternative method to CFD, which requires less calculation time and computer memory, was developed by Daoud et al. [9–11]. It combines a zonal airflow model, a radiation model, a humidity transport and condensation model and takes into account resurfacing and occupation. This above ice model (AIM) predicts the heat fluxes through the envelope as well as the temperature and absolute humidity distributions for a 3D transient regime during an entire typical meteorological year. In particular it calculates the heat fluxes into the ice sheet by convection, radiation and condensation. The temperature below the ice sheet was assumed uniform and constant. The results show a satisfactory agreement with corresponding measurements and CFD calculations. The present article describes a second part in the development of a global 3D transient model of an ice rink. The below ice model (BIM) was developed using an implicit unidirectional electrical analogy, taking into account the secondary loop and brine movement and the heat gain from the ground (with changing meteorological conditions). The BIM was coupled successfully with the previously mentioned AIM. The combined model eliminates the assumption of constant temperature below the ice sheet used in AIM. Instead the temperature of the brine entering the pipes below the ice sheet must be specified. The combined model evaluates the return brine temperature, the total refrigeration load, the ice surface temperature (IST), the heat gain from ground, as well as the energy consumption of the ventilation system and of the radiant heaters. Parametric studies were undertaken in order to evaluate the impact of the climate, brine inlet temperature, ice thickness and other parameters on the calculated results and their results are presented in the last part of the present paper. 2. Description and modeling 2.1. Ice rink description Figs. 1 and 2 show a schematic representation of the studied ice rink ‘‘Camilien Houde’’ located in Montreal (Canada). The building is 64.2-m long, 41.5-m wide and its height is 9.2 m. The ice surface is 61-m long, 25.9-m wide and is surrounded by a narrow corridor. The space above the stands is heated by eight radiant heaters (22 kW  8) which are controlled by a thermostat. Seven inlets supply a stream of ventilation air. Its flow rate is 4270 L/s except during resurfacing of the ice when it is increased to 10,384 L/s to evacuate the combustion gases of the resurfacing vehicle. The air exits through four outlets on the walls. Heat gains from lighting are 10 W/m2 above the ice and 5 W/m2 above the stands; those due to the presence of the audience are also taken into account while the number of spectators is specified according to a weekly schedule [7,11]. The ice resurfacing takes place several times per day, lasts 12 min and is modeled as a 1 mm film of hot water at 60 8C. Its frequency, specified in the schedule mentioned above, is higher in the evenings and weekends. The stands, corridors and boards are also modeled in the AIM. 502 L. Seghouani et al. / Energy and Buildings 41 (2009) 500–511 Fig. 1. Schematic section of the ice rink and the different layers under the ice (not to scale). Fig. 2. Top view of the ice showing the different zones and the flow of the brine. The ground structure beneath the ice rink is represented in Fig. 1 and comprises horizontal layers of ice (50 mm), concrete (150 mm), thermal insulation (100 mm), sand (200 mm) and, finally, soil. The total depth of this structure included in the calculation domain is 4 m. The secondary coolant used to maintain the ice at the desired temperature is calcium chloride brine. It is supplied from a header located at the west end of the ice sheet and circulates in the concrete slab at a depth 57.5 mm below the ice surface within 74 uniformly distributed, four-pass polyethylene tubes (25 mm ID). The spacing between tubes is 87.5 mm. The main collector has an internal diameter of 150 mm. The flow rate of the pump is 28.5 L/s. An electrical heater of 8 kW is activated in the sand layer when the ground temperature at a depth of 4 m is below 4 8C to prevent freezing under the concrete slab which can cause damage to the underground structure and ice. 2.2. Model of air movement and heat exchanges above the ice (AIM) The air movement and heat exchanges occurring in the rink above the ice surface are simulated using a 3D transient model with 64 zones [9–11]. It consists of six coupled submodels solved with the ‘‘onion’’ method. The first submodel is the energy model which uses the Multizone Building Model (type 56) of TRNSYS [12]. It is based on two relations. The first one expresses energy conservation inside each thermal zone i: ðM C p Þi X dT i ṁ j ! i C p T j ¼ ðqcv AÞi þ dt j (1) while the second expresses energy conservation for each internal surface in contact with the air in the building: qcd ¼ qcv þ qrd þ qcond þ qrs (2) The conductive flux through the wall is evaluated using the transfer functions method while the convection flux between the wall surface and the air inside the building is calculated using a constant heat transfer coefficient (3 W/m2 K). The radiation flux between internal surfaces of the building is provided by a submodel (radiative transfer submodel) based on the Gebhart method [13]. The condensation flux is attributed to the ice surface when its temperature is below the dew point of the air above it. It is provided by a submodel (humidity transport submodel) which calculates the absolute humidity of the air in every thermal zone inside the building using the following conservation equation: M air;i dvi X ¼ ṁi; j ðv j  vi Þ þ Ẇg;i dt i; j (3) Finally, qrs corresponds to the heat flux occurring when the resurfacing operation deposits approximately 0.5 m3 of water at L. Seghouani et al. / Energy and Buildings 41 (2009) 500–511 503 Fig. 3. Schematic representation of the ventilation system. 60 8C on the ice surface. It is calculated using the equation recommended by ASHRAE [14]. The airflow mi,j between thermal zones used in Eq. (1) is provided by the zonal airflow submodel. The formulation used expresses the mass flow crossing the common surface between two zones i and j in terms of their pressure difference. Thus, in the case of a vertical interface ṁi; j ¼ ei; j C d pffiffiffiffiffiffiffiffi 2ri Ai; j jP j  P i j1=2 (4) While in the case of a horizontal interface ṁi; j ¼ ei; j C d  1=2 pffiffiffiffiffiffiffiffi   1 2ri Ai; j P j  P i  ðri ghi þ r j gh j Þ 2 (5) The coefficient eij is equal to +1 when flow is from zone i to zone j and equal to 1 for flow from zone j to zone i. A new submodel not described in our previous publications [9,10] is used to simulate the behaviour of the ventilation system. It consists of two units (see Fig. 3). The first one is used for cooling, dehumidifying and reheating the external air when its temperature is above 23 8C while the second unit is for heating and humidifying it when its temperature is below 15 8C. When the temperature of the entering air is between 23 8C and 15 8C none of the units is in operation unless the humidity level is too high or too low. The humidity controls are such that the relative humidity of the air entering the ice rink is maintained between 20% and 33%. The equations modeling the operation of these two units are the mass and energy conservation equations for a gas–vapour mixture in a steady state, steady flow process. It should be noted that the results presented here assume that only outdoor air is handled by the ventilation system (no recirculation). The final submodel evaluates the ventilation effectiveness by calculating the age of the air in every zone of the ice rink. The data exchange between these six submodels is represented in Fig. 4. It takes place several times at every timestep until the outputs of each submodel vary by less than 103. 2.3. Below ice modeling (BIM) The modeling of the ice rink ground structure shown in Figs. 1 and 2 is of a substantial nature. Indeed, the concrete slab is one of the most important parts of an ice rink. It forms and maintains the ice by removing heat from it, using a secondary coolant (brine) circulating in the embedded network of pipes. The BIM is based on the transient one-dimensional conduction equation and the electrical circuit analogy. For that, the ice surface is divided in eight equal square surfaces (one to eight) shown in Fig. 2, which correspond to those used by the zonal submodel of the AIM. Each of these surfaces is subdivided into two parts (A and B). The brine flows from west to east under part A of the eight surfaces and in the opposite direction under part B (see Fig. 2). Thus, the brine enters at 1A or 2A with a constant temperature Tb,in and exits at 1B or 2B with a temperature Tb,out which is not the same for the North and South brine loops. Heat transfer between the tubes is neglected since the brine temperature increases by less than 2 8C between inlet and outlet. This justifies the hypothesis of one-dimensional vertical heat transfer between the different layers of the ground structure. However, the BIM takes into consideration heat exchanges between the horizontal layers and the outdoor or the indoor appropriate ground surface as explained below. Therefore, it takes into account three-dimensional phenomena. Fig. 5 shows the equivalent electrical circuit under one-half of the ice. At each of the 16 subdivisions of the ice surface shown in Fig. 2 the heat flux calculated by AIM enters the ice at node 1. The electrical heating in the sand layer is added to nodes 5 and 6. The temperature of node 7 is considered as given by the following correlation based on the deep soil (4 m) temperature data for Montreal [15] T 7 ¼ 3:34  4:78  103 t  1:97  107 t 2 þ 1:15  109 t 3  2:66  1013 t 4 þ 2:15  1017 t 5  5:91  1022 t 6 þ 273:15 (6) where t is in h and T7 in K. The heat exchanges between the horizontal layers and the outdoor or indoor ground surface are added at nodes 2–6. They are evaluated using the ASHRAE method [16] and are divided in two halves attributed to the circuits under parts A and B. This method uses the formula: Q node ¼ U node Snode ðT o  T node Þ where the below-grade average U-factor is given by:      2ksoil 2k R 2k R U node ¼  ln z1 þ soil  ln z2 þ soil pðz2  z1 Þ p p (7) (8) For nodes 4, 5 and 6 on the north, east and west sides of the ice rink To is the temperature of the ground surface outside the ice rink. It is calculated at each time step by the following correlation based on the surface soil temperature data for Montreal [15] T gr ¼ 2:56  1:31  102 t  6:88  106 t 2  1:79  1010 t 3  2:02  1013 t 4 þ 2:42  1017 t 5  7:93  1022 t 6 þ 273:15 Fig. 4. Information flow in the above-ice model (AIM). where t is in h and Tgr in K. (9) 504 L. Seghouani et al. / Energy and Buildings 41 (2009) 500–511 Fig. 5. Thermal circuit model for below ice structure (1A, 1B. . . etc. identify zones in Fig. 2). On the south side To is taken as constant (equal to the temperature in the changing rooms). On the three other sides, for nodes 2, 3 and 4, To is taken as the temperature of the concrete surface in the corridors surrounding the ice surface which is calculated by the AIM. The brine enters each subdivision and absorbs the heat coming from nodes 2 and 4. Its temperature therefore increases and it enters the next subdivision where the process is repeated. Each of the 16 subdivisions (parts A and B of surfaces 1–8 in Fig. 2) is modeled in the same manner and only the entering brine temperature and the lateral heat exchanges calculated by Eq. (7) vary from one to another. Thus there are 7 unknown temperatures under each of the 16 subdivisions (at nodes 1–6 and at the brine outlet) or 118 unknown temperatures altogether. The implicit discretisation of the transient energy balance for each node under each subdivision leads to a system of linear equations which can be represented by the matrix equation AT ¼ B (10) where T is the vector of the seven unknown temperatures. The expressions of the 7  7 matrix A and of the vector B are given in the appendix. This system of linear equations is solved by inversing matrix A since this direct method is suitable for small-size systems. The solution starts under surface 1A and then continues to 3A, 5A, etc. following the direction of flow. The temperature of all the nodes is thus obtained for each of the 16 subdivisions. Then, the temperature of each level below the ice is obtained by averaging the 16 corresponding node temperatures. The final outlet temperature of the brine is calculated by assuming that the two streams from 1B and 2B are mixed adiabatically. The heat rates from the ice to the brine and from the soil to the brine are also evaluated and therefore the total refrigeration load is calculated. A FORTRAN subroutine of this model was incorporated as a new type in TRNSYS [12]. 2.4. Coupling of the AIM and BIM The coupling of the BIM with the AIM was realised using the ‘‘Onion’’ method. Fig. 6 shows a schematic representation of the coupling method. During one time step, the BIM calculates the temperature T2 between the ice and the concrete for each of the eight zones. It provides them as inputs to the AIM in which they are used as boundary conditions. This model calculates the eight total heat fluxes towards the ice and returns them as inputs to the BIM. Several such data exchanges take place until outputs of each model vary by less than 103. Then time is incremented and this procedure is repeated. For the simulations performed in the present project hourly average weather data for a typical meteorological year (temperatures of the air and the ground at 0 m and 4 m depth as well as solar L. Seghouani et al. / Energy and Buildings 41 (2009) 500–511 Fig. 6. Schematic representation of the onion coupling between AIM and BIM. radiation [15]) is used. In order to ensure the periodicity of the results (values at the end of the 365th day must be identical to those at the beginning of the first one) we carried out simulations over 17 months by repeating the meteorological data for January to May and did not consider the results of the initial 5 months to be sure that the effects of the arbitrary initial values are eliminated. The choice of time step is based on four considerations. The first is that TRNSYS accepts only time steps of the form 1/N where N is an integer. The second is the duration of resurfacing (12 min) which means that the time step should not exceed 0.2 h. The third is the need to minimize the time required for the calculation of the results over the 17-month period which implies that the time step must be large. Finally, a small time step is necessary to capture more transient details. As a compromise between these considerations we have opted for a time step equal to 0.1 h. With this choice the entire yearly simulation requires about 80 h on a personal computer (with an Intel Core 2 Duo 6400 2.13 GHz processor and 2 Go of RAM). This is a considerable improvement over an equivalent CFD two-dimensional simulation which requires 24 h to calculate the results for a single day [7]. 2.5. Model validation The AIM was successfully validated in previous studies [9–11]. The predictions of the more complete numerical model presented here are validated by comparison with measurements recorded in the Camilien Houde ice rink over relatively short periods in 2005 and 2006 [8]. Table 1 presents such a comparison of the measured and calculated brine outlet temperature and ice surface temperature. Measured values are the averages of recorded temperatures while calculated values are monthly averages for the typical meteorological year. Although these conditions are not identical, the small seasonal variation of these results makes this comparison acceptable. The agreement between measured and calculated values of these temperatures shows that the proposed model predicts satisfactorily these values. The maximum relative difference is less than 10%, which is acceptable since the simulations did not take into account the refrigeration system which regulates the inlet brine temperature. This difference is principally due to the boundary condition used for the simulation (constant brine inlet temperature Tb,in = 9 8C) and to imprecision of the measurements. Similar agreement has also been established between the measured and calculated values of the heat flux into the ice. The former was obtained by integrating the readings of four heat flux sensors installed under the ice sheet. The corresponding daily mean values for 1 October 2005 are: 94.9 W/m2, 56.6 W/m2, 90.6 W/m2 and 113.0 W/m2 [8]. Disregarding the second sensor which gives significantly lower values than the other three, the mean value of the experimentally measured heat flux is 99.5 W/ m2. On the other hand, the model predicts an average total heat flux of 108 W/m2 for October which is about 8% higher than the measurements. In view of the fact that these values do not correspond to identical climate conditions, their agreement is judged to be acceptable. In view of these results we consider that the proposed model can be used with confidence for the calculation of typical yearly refrigeration loads and for parametric studies which aim to establish the effect of design and operation conditions on these loads. 3. Parametric analysis The analysis in this section starts from a ‘‘base case scenario’’ which corresponds to meteorological conditions for a typical year in Montreal (latitude N 458470 , longitude W 738750 ), a constant brine inlet temperature to the slab equal to 9 8C, an ice thickness of 5.08 cm, an under slab insulation thickness equal to 10 cm, a setpoint of 15 8C with a nocturnal set back of 7 8C and an hysteresis of 0.2 8C for the electronic thermostat, while the underground heating (8 kW) is in operation when the ground temperature at a depth of a 4 m is below 4 8C. The results of this transient simulation are compared with corresponding results obtained by varying the parameters defining the base case scenario one at a time. 3.1. Effects of the climate Results have been calculated for typical meteorological years for the cities of Edmonton (latitude N 538310 , longitude W 114850 ), Houston (latitude N 298580 , longitude W 958220 ) and Pittsburgh (latitude N 408300 , longitude W 808130 ) for which correlations similar to those for Montreal (Eqs. (6) and (9)) have been obtained from the meteorological data [15]. These results are compared with those for the base case at Montreal. As indicated by the values in Table 2 these climatic conditions vary from very warm and humid during the summer in Houston to very cold and dry in winter in Edmonton. Figs. 7 and 8 illustrate the effects of the climate on the average daily energy consumption of the ventilation system which, as mentioned before, is constituted of two units (cf. Fig. 3). In particular, Fig. 7 shows that the energy consumption of the first unit, which dehumidifies the ventilation air by cooling and reheating it, is greatest for Houston where it is significant throughout the year. On the other hand, for the other three cities this quantity is essentially zero during the winter months but becomes important during the summer. The seasonal variation of this quantity as well as its relative magnitude between the four cities under consideration is consistent with the meteorological data of Table 2. The same is true for the results in Fig. 8 which shows the energy consumption of the second unit of the ventilation system. This one heats and humidifies the incoming ventilation air and is therefore high in winter and essentially zero in summer, largest for Edmonton and lowest for Houston. It is Table 2 Comparison of climatic conditions in the cities under consideration. Table 1 Comparison between measured and calculated temperatures. Return brine temperature (8C) Ice surface temperature (8C) Measured Calculated 7.3 From 5.5 to 6.3 From 7.8 to 8.0 From 5.2 to 6.0 505 Edmonton Houston Montreal Pittsburgh Heating dry bulb (99%) Cooling dry bulb (2%) Mean coincident wet bulb 30.5 8C 0.4 8C 21.8 8C 14.1 8C 24.0 8C 33.5 8C 25.8 8C 28.7 8C 15.7 8C 24.9 8C 19.5 8C 20.6 8C 506 L. Seghouani et al. / Energy and Buildings 41 (2009) 500–511 Fig. 7. Effect of the climate on the cooling and reheating energy consumption of the ventilation system. Fig. 8. Effect of the climate on the heating and humidification energy consumption of the ventilation system. interesting to note that the peak energy consumption for cooling– reheating is for every city higher than the one for heating– humidification. The predicted total annual energy consumption of the ventilation system for Edmonton, Houston, Montreal and Pittsburgh is 1278 MWh, 4033 MWh, 1435 MWh and 2228 MWh, respectively. The distribution of this total among the four processes is shown in Fig. 9. Cooling uses the largest part in Houston, Pittsburgh and, surprisingly, Montreal while in Edmonton the largest part is used for heating. Humidification requires a very small part of the total energy consumption in all cases. The impact of the climate on the temperature of the ice surface is quite small (1 8C) since this variable is essentially determined by the temperature of the brine which for these simulations is the same throughout the year for all four cities. Therefore the values of IST are not presented here. However, it should be noted that these values are somewhat influenced by the temperature of the air above the ice which depends on whether the radiant heating is on and whether the ventilation air is heated or cooled. Thus the IST for all cities is slightly higher in winter, when heating is required, than in summer, when the ventilation air requires cooling. Similarly, the IST during the summer is a little bit higher in Edmonton than in Houston since the former city requires much less cooling of the ventilation air and necessitates heating of the stands as shown in Fig. 10. This last figure also shows that the radiant heating of the stands exhibits the expected seasonal behaviour, that it is significantly smaller in Houston and quite important in Edmonton even during the summer. It should be noted that the daily energy consumption of the radiant heaters is for all cities significantly lower than the corresponding values of the ventilation system (sum of consumption shown in Figs. 7 and 8). Fig. 11 shows the influence of the climate on the refrigeration load, i.e. on the sum of the heat reaching the brine from the ice and from the ground. This quantity is greatest for Houston where the seasonal variation is also the least pronounced. This result is consistent with the climatic conditions which signify that the Fig. 9. Effect of the climate on the annual energy consumption of the ventilation system. L. Seghouani et al. / Energy and Buildings 41 (2009) 500–511 507 Fig. 10. Effect of the climate on the energy consumption of the radiant heating. Fig. 11. Effect of the climate on the refrigeration load. envelope and indoor air are warmer in Houston and that, therefore, heat fluxes towards the ice by both radiation and convection are also higher in this city. Furthermore, because the humidity is higher in Houston the corresponding heat flux due to condensation of water vapour on the ice is also higher in this city. The combination of all these effects explains the results shown in this figure. Table 3 compares the annual refrigeration load and the corresponding energy consumption of the different systems for each city. It confirms the fact that the refrigeration load is greatest for Houston and shows that this quantity is always smaller than the energy consumption by the ventilation system. It is important to note that the refrigeration load does not vary significantly (less than 7.5% of the average value) between these four locations despite their very different climates. This is attributed to the fact that the indoor conditions are quite similar due to the controls on the ventilation system and the radiant heaters. However, it is expected that the energy consumption of the refrigeration system will be higher in locations with warm climates since the condensation temperature and pressure will be higher than in locations with cold climates. Table 3 also shows that the energy consumed by the lights and the brine pump does not depend on the climate in accordance with the modeling assumptions. Lighting consumes more energy than the radiant heaters for each of the four locations under consideration while the same is true for the energy used by the brine pump in three of the four cities. The greatest and smallest contributors to the total energy consumption are the ventilation system and the underground heater respectively for each of the four cities under consideration. Finally, the total annual energy consumption is highest, by a considerable margin, in Houston (where summer operation requires considerable quantities of energy for cooling and dehumidification) and lowest in Edmonton. Table 4 shows the energy consumption by the underground heater which is activated to avoid freezing and heaving of the ground under the concrete slab. Its operation is not at all necessary in the case of the warmer climates (Houston and Pittsburgh) while in Edmonton and Montreal it is in operation for approximately 7 and 5 months, respectively. 3.2. Effects of the hysteresis of the thermostat The influence of this parameter is established by comparing the results for an electronic thermostat with a hysteresis of 0.2 8C with those for a conventional bimetallic thermostat with a hysteresis of 1.5 8C. The set point is 15 8C in both cases. Fig. 12 shows that the refrigeration load is always higher in the case of the electronic thermostat. The difference varies from Table 3 Annual energy consumption of different systems and annual refrigeration load for each city (MWh). Annual refrigeration load Edmonton Houston Montreal Pittsburgh 1023 1110 979 1018 Annual energy consumption of different systems Ventilation system Radiant heating Underground heating Lighting Brine pump Total 1278 4033 1435 2228 133 30 86 74 40.5 0 30.1 0 148.7 148.7 148.7 148.7 98.1 98.1 98.1 98.1 1698.3 4309.8 1797.9 2548.8 Table 4 Monthly underground electrical heating for each city (MWh). L. Seghouani et al. / Energy and Buildings 41 (2009) 500–511 508 Fig. 12. Effect of the thermostat hysteresis on the refrigeration load. Fig. 14. Effect of the brine inlet temperature on the refrigeration load. approximately 125 kWh/day in summer to 200 kWh/day in winter. This result can be explained by the fact that the radiant heating is turned on more frequently with an electronic thermostat in order to maintain the temperature in the zone occupied by the spectators at (15  0.2) 8C; therefore the corresponding average air and envelope temperatures are slightly higher and result in increased convective and radiative fluxes towards the ice. However, the energy savings associated with the conventional thermostat are achieved at the expense of the spectators comfort since in that case the air temperature above the stands oscillates between 13.5 8C and 16.5 8C. The smaller difference in summer energy consumption between the two cases under consideration is due to the fact that the radiant heating load is greatly reduced during these warm months. It should also be noted that the type of thermostat influences the ice temperature, i.e. its quality. Thus the temperature of the ice surface is higher for the base case (electronic thermostat) by approximately 0.25 8C in winter and 0.15 8C in summer. elements is set to 15 8C during the day and 7 8C during the unoccupied night hours, while in second case the thermostat is set to 15 8C throughout the day. Fig. 13 shows that the use of a thermostat with nocturnal set back (base case) reduces the refrigeration load. However, this reduction is very small (it varies from 25 kWh/day during the winter to 50 kWh/day during the summer). It is due to the decrease of the operation time of the radiant heating elements and the corresponding reduction of the air and envelope temperatures which in turn lower the convective and radiative fluxes towards the ice. The corresponding effect on the IST is insignificant (reduction of 0.1 8C when nocturnal set back is used) since these flux reductions are small compared to their respective values. Finally, it is important to note that the heat flux from the ground to the brine is totally unaffected by the use, or not, of the nocturnal set back since the conditions below the concrete slab are independent of those prevailing within the ice rink. 3.3. Effect of the nocturnal set back 3.4. Effect of the brine inlet temperature Two cases are compared in this section. The first one is the base case in which the thermostat controlling the radiant heating Three different constant brine inlet temperature were used for this study (Tb,in = 8 8C, 9 8C and 10 8C). Fig. 14 shows that the Fig. 13. Effect of the nocturnal set back on the refrigeration load. Fig. 15. Effect of the ice thickness on the refrigeration load. L. Seghouani et al. / Energy and Buildings 41 (2009) 500–511 509 Table 5 Effect of the ice thickness on the January temperatures of the ice and the brine. Ice thickness (cm) 2.54 5.08 7.62 Temperatures (8C) (January) Ice surface Brine (inlet) Difference 6.20 5.97 5.74 9.0 9.0 9.0 2.8 3.03 3.26 Table 6 Effect of the ice thickness on the July temperatures of the ice and the brine. Ice thickness (cm) 2.54 5.08 7.62 Temperatures (8C) (July) Ice surface Brine (inlet) Difference 5.46 5.13 4.83 9.0 9.0 9.0 3.54 3.87 4.17 refrigeration load increases when the inlet brine temperature decreases. The increase is about 100 kWh/day for each decrease of Tb,in by 1 8C. This behaviour is qualitatively logical since the brine acts as a heat sink while the conditions of the heat sources above and below the concrete slab remain the same for these three simulations. A 1 8C increase or decrease of the inlet brine temperature causes, respectively, an augmentation or diminution of 0.8 8C of the ice surface temperature, for all months of the year. An analogous observation can be made for the outlet brine temperature which increases or decreases by the same amount as the inlet brine temperature. Fig. 16. Effect of the insulation thickness on the heating rate from the soil. 3.5. Effect of the ice thickness Fig. 15 shows that the refrigeration load decreases slightly with the increase of the ice thickness. It is reduced by 8–10% when the ice thickness triples (from 2.54 cm to 7.62 cm). This is due to the corresponding increase of the thermal resistance of the ice sheet which also causes a small increase of the IST. The increase of the latter is greater in summer due to the augmentation of the convection and radiation fluxes during this warm season. Furthermore, Tables 5 and 6 show that, when the ice thickness is increased by 5 cm, the temperature difference between the ice surface and the brine increases by 14% (from 2.8 8C to 3.26 8C) for the month of January and by 18% (from 3.54 8C to 4.17 8C) for the month of July. These results should not be generalised. They have been calculated by fixing the brine inlet temperature and do not necessarily apply to other control strategies (such as, for example, cases where the ice surface temperature is kept constant). 3.6. Effect of the insulation thickness Fig. 16 shows that an increase of the underground insulation thickness from 10 cm to 30 cm reduces the heat rate transferred to the brine pipes from the soil by approximately 25% for all months of the year. It should also be noted that the change of insulation thickness has an insignificant effect on the heat rate from the ice to the brine pipes. However, the effect of this reduction of the heat rate from the soil on the total refrigeration load is not significant since the former is less than 10% of the latter (cf. Fig. 15). The most important impact of added underground insulation is to reduce the risk of freezing under the concrete slab which can cause heaving and damage to the ice. Fig. 17 which shows the average temperature of node 5 situated at the sand–insulation interface clearly illustrates this effect. Furthermore, as this insulation thickness is increased the energy consumption of the electrical element in the sand can be reduced without increasing the danger of freezing. Fig. 17. Effect of the insulation thickness on the temperature of the insulation–sand interface. 4. Non-linear correlations of energy loads This section presents four correlations between the monthly mean value of the energy consumption (in kWh/day) of the different processes taking place in the ventilation system and the corresponding sol-air temperature (in 8C) calculated with the expression and parameters proposed by ASHRAE [17]. They are based on numerous simulations carried out at the four selected cities (Montreal, Edmonton, Houston and Pittsburgh). Cooling load correlation: 2 Q Cool ¼ 283:96 þ 189:1 T sol-air þ 4:85 Tsol -air (14a) Reheating load correlation: 2 Q Re heat ¼ 8:63 þ 192:58 T sol-air  2:217 Tsol -air (14b) Heating load correlation: 2 Q Heat ¼ 1943:1  140:23 T sol-air þ 1:89 Tsol -air (14c) Humidification load correlation: 2 Q Humid ¼ 44:9  10:95 T sol-air þ 0:79 Tsol -air (14d) Figs. 18 and 19 show that these fairly simple correlations agree quite well with the calculated values. These results also show that 510 L. Seghouani et al. / Energy and Buildings 41 (2009) 500–511  The annual refrigeration load does not vary significantly (less than 7.5% from the mean value) between these four locations despite their very different climates.  The annual energy consumption by the ventilation system is significantly influenced by the meteorological conditions; it is highest in hot and humid locations and is always greater than the corresponding refrigeration load.  The annual energy consumptions by the radiant heaters and by the underground electric heater are very small compared to the refrigeration load even in the coldest of the four cities under consideration. The main results of the other parametric studies are the following: Fig. 18. Energy consumption of ventilation system (cooling/reheat) vs. sol-air temperature.  The use of an electronic thermostat with a low hysterisis (0.2 8C) rather than a bimetallic one with a high hysteresis (1.5 8C) increases the refrigeration load by more than 10% if the same set point is used in both cases; however the energy savings associated with the bimetallic thermostat are achieved at the expense of the spectator’s comfort.  The use of a thermostat with nocturnal set back results in a small reduction of the refrigeration load.  A reduction of the brine temperature causes an increase of the refrigeration load and an almost equal reduction of the ice surface temperature.  An increase of the ice thickness causes a decrease of both the refrigeration load and the ice surface temperature.  An increase of the underground insulation thickness causes a slight decrease of the refrigeration load and reduces significantly the danger of ground freezing and heaving. Finally, four correlations expressing the energy consumptions for heating and humidifying, or cooling and reheating, the ventilation air in terms of the sol-air temperature have been proposed. Acknowledgements Fig. 19. Energy consumption of the ventilation system (heating/humidification) vs. sol-air temperature. the energy consumption by the heating and cooling processes depend heavily on Tsol-air while those by the reheat and, in particular, humidification processes are less dependent on this parameter. 5. Conclusion A transient model which calculates heat transfer through the ground towards the brine pipes imbedded in the concrete slab under the ice of an indoor ice rink has been formulated and coupled with a previously developed one which calculates heat fluxes in the building by convection, radiation and phase changes. The resulting simulation tool has been enriched with subroutines which calculate the energy consumption for heating and humidifying, or cooling and reheating, the ventilation air. After validation with experimental results, this tool was used to evaluate the refrigeration load as well as the energy consumed by the radiant heaters of the stands, each process of the ventilation system, the lights, the brine pump and the underground electric heater over a typical year at four North American locations (Edmonton, Houston, Montreal, Pittsburgh) with very different meteorological conditions. The results of this analysis show that: This study was financed by the Natural Sciences and Engineering Research Council (NSERC) of Canada through the Strategic Project Grant STPGP 306792 (Title: Development of design tools and of operation guidelines for the heating, ventilation, air conditioning and refrigeration systems of ice rinks). Appendix A The expressions of the 7 T 7 matrix A and of the vectors B and T in Eq. (10) are: 1 1 0 0 0 0 0 C B Að1; 1Þ  R Ice C B C B 1 1 C B Að2; 2Þ 0 0 0 0 C B R Rc1 B Ice C B C 1 1 B 0 C Að3; 3Þ 0 0  ṁ C B PB B C R R c1 c2 B C C 1 1 A¼B B 0 C 0 Að4; 4Þ 0 0 B C R R Ins c2 B C B C 1 1 C B 0 Að5; 5Þ 0 0 0 C B RIns RSand C B C B 1 C B 0 0 0 0 Að6; 6Þ 0 A @ RSand 0 0 2 0 0 0 1 0 (A.1) L. Seghouani et al. / Energy and Buildings 41 (2009) 500–511 1 ðM C p Þeq1 p Q Ice þ T1 B C B C Dt B C ðM C Þ p eq2 p 1 B C T2  U 2 S2 T gr  B C B C 2 D t B C ðM C Þ p 1 B eq3 p C B ṁB C pB T In  U 3 S3 T gr  T3 C B C 2 Dt B C ðM C p Þeq4 p B¼B C 1 B C U  S T  T gr 4 4 4 B C 2 D t B C ðM C p Þeq5 p B C 1 B C T5 Q H  U 5 S5 T gr  B C 2 D t B C B C ðM C p Þeq6 p 1 1 B Q  U S T gr  T7 C T6  6 6 H @ A 2 RSoil Dt T in 0 ðM C p Þeq1 Að2; 2Þ ¼  Að3; 3Þ ¼  Að4; 4Þ ¼  Að5; 5Þ ¼  Að6; 6Þ ¼  1 RIce ! ðM C p Þeq2 1 1 1 þ þ U 2 S2 þ RIce Rc1 2 Dt ! ðM C p Þeq3 1 1 1 þ þ U 3 S3 þ Rc1 Rc2 2 Dt ! ðM C p Þeq4 1 1 1 þ þ U 4 S4 þ RIns Rc2 2 Dt ! ðM C p Þeq5 1 1 1 þ þ U 5 S5 þ RSand RIns 2 Dt ! ðM C p Þeq6 1 1 1 þ þ U 6 S6 þ RSoil RSand 2 Dt Dt 1 1 MIns C PIns þ M Sand C PSand 2 2 (A.5e) ðM C p Þeq6 ¼ 1 1 M CP þ M CP 2 Sand Sand 2 Soil Soil (A.5f) References (A.3) Here the coefficients in Eq. (A.1) are Að1; 1Þ ¼ ðM C p Þeq5 ¼ (A.2) and Pþ1 T ¼ T1pþ1 T2pþ1 T3pþ1 T4pþ1 T5pþ1 T6pþ1 TOut 511 þ (A.4) And the equivalent thermal capacity in Eqs. (A.2) and (A.4) is calculated as follows: ðM C p Þeq1 ¼ 1 M Ice C PIce 2 (A.5a) ðM C p Þeq2 ¼ 1 1 M Ice C PIce þ Mc1 C PC 2 2 (A.5b) ðM C p Þeq3 ¼ 1 M B C PB 2 (A.5c) ðM C p Þeq4 ¼ 1 1 M c2 C PC þ M Ins C PIns 2 2 (A.5d) [1] M. Lavoie, R. Sunyé, D. Giguére, Potentiel d’économies d’énergie en réfrigération dans les arénas du Québec, Report prepared by the CANMET Energy Technology Center, 2000. [2] P.J. Jones, G.E. Whittle, Computational fluid dynamics for building air flow prediction—current status and capabilities, Building and Environment 27 (3) (1992) 321–338. [3] Z. Jian, Q. Chen, Airflow and air quality in a large enclosure, ASME Journal of Solar Energy Engineering 117 (2) (1995) 114–122. [4] C. Yang, P. Demokritou, Q. Chen, J.D. Spengler, A. Parsons, Ventilation and air quality in indoor ice skating arenas, ASHRAE Transactions 106 (2) (2000) 338– 346. [5] O. Bellache, M. Ouzzane, N. Galanis, Coupled conduction, convection, radiation heat transfer with simultaneous mass transfer in ice rinks, Numerical Heat Transfer Part A 48 (2005) 219–238. [6] O. Bellache, M. Ouzzane, N. Galanis, Numerical prediction of ventilation and thermal processes in ice rinks, Building and Environment 40 (3) (2005) 417– 426. [7] O. Bellache, N. Galanis, M. Ouzzane, R. Sunyé, D. Giguére, Two-dimensional transient model of airflow and heat transfer in ice rinks (1289-RP), ASHRAE Transactions 112 (2) (2006) 706–716. [8] M. Ouzzane, R. Sunyé, R. Zmeureanu, D. Giguére, J. Scott, O. Bellache, Cooling load and environmental measurements in a Canadian indoor ice rink (1289-RP), ASHRAE Transactions 112 (2) (2006) 538–545. [9] A. Daoud, N. Galanis, Calculation of the thermal loads of an ice rink using zonal model and building energy simulation software, ASHRAE Transactions 112 (2) (2006) 526–537. [10] A. Daoud, N. Galanis, O. Bellache, Calculation of refrigeration loads by convection, radiation and condensation in ice rinks using a transient 3D zonal model, Applied Thermal Engineering (2007), http://dx.doi.org/10.1016/j.applthermaleng.2007. 11.011. [11] A. Daoud, Analyse des transferts de chaleur et de masse transitoires dans un aréna à l’aide de la méthode zonale, Thèse de doctorat en sciences appliquée, Université de Sherbrooke (Québec), Canada, 2007. [12] TRNSYS, A transient system simulation program, Software Manual, Solar Energy Laboratory, University of Wisconsin-Madison, 2000. [13] B. Gebhart, Heat Transfer, Second edition, McGraw-Hill, 1971. [14] ASHRAE, Ice Rinks, ASHRAE Handbook-Refrigeration, American Society of Heating, Refrigeration, and Air-Conditioning Engineers, Inc., Atlanta, 2002. [15] ENERGYPLUS, http://www.eere.energy.gov/buildings/energyplus/. [16] ASHRAE, Residential cooling and Heating load Calculations ASHRAE Handbookfundamentals, American Society of Heating, Refrigeration, and Air-Conditioning Engineers, Inc., Atlanta, 2005, pp. 29.11–29.12. [17] ASHRAE, Non-residential cooling and Heating load Calculations ASHRAE Handbook-fundamentals, American Society of Heating, Refrigeration, and Air-Conditioning Engineers, Inc., Atlanta, 2001, pp. 29.14–29.17.