Two-dimensional transient model of airflow and heat transfer in ice rinks.

ABSTRACT

A two-dimensional transient model of air movement as well as heat and mass transfer in an ice rink was developed and tested by comparing its predictions with measured values in a Montreal municipal ice rink. The model was then used to predict the daily heat flux profiles into the ice by convection, radiation, condensation, and resurfacing operations for a representative configuration and for average, steady, and periodic meteorological conditions. The heat loads from the ground and dissipation in the floor pipes were also calculated. The effects of climate and some design parameters or operating conditions (ceiling emissivity, below-floor insulation thickness, and average brine temperature) on the daily loads were evaluated and analyzed.

INTRODUCTION

Background

The several thousand indoor ice rinks of North America are used for spectator events such as ice sports and ice shows. They are different sizes (ASHRAE 2002) but are characterized by the following common features: large dimensions with few or no dividing walls; significant refrigeration load; simultaneous need of heating (or air conditioning), refrigeration, and ventilation; high energy consumption; and, finally, considerable emissions of greenhouse gases (GHG) due to energy consumption and synthetic refrigerant losses. The state (temperature, humidity, etc.) and motion of the ventilation air through these large buildings often lead to uncomfortable thermal conditions for the spectators and possible formation of fog as well as an increase of the refrigeration load and a deterioration of ice conditions. A recent study (Lavoie et al. 2000) indicates that in Quebec alone ice rinks offer a potential energy conservation of 1300 GWh (4436 x [10.sup.6] Btu) annually and a possible GHG reduction of 0.5 MTons of C[O.sub.2] equivalent per year. According to the study's authors, the current performance of ice rinks is poor due to the lack of design and operational guidelines taking into consideration the interaction between the ice-making, ventilation, and heating systems.

For its part, ASHRAE Technical Committee 10.2, Automatic Icemaking Plants and Skating Rinks, in a work statement (MacCracken and Scott 2003) stated that "the available information from ASHRAE is of limited usefulness for the design and operation of ice sheet cooling systems by the general HVAC & R design community or arena owner/operators." It therefore issued a request for proposals to "develop and verify methods for determining ice sheet cooling loads," which was eventually awarded to the CANMET Energy Technology Centre-Varennes. This paper presents the adopted methodology and initial results obtained while performing task 1 ("Develop a two-dimensional transient numerical model of ice sheet heat transfer") and task 2 ("Numerical model validation") of the project.

General Description of the Adopted Approach

The adopted model uses a detailed description of the skating rink and typical steady-periodic ambient conditions (ASHRAE 2005) to calculate the corresponding profiles of the heat loads (conductive, convective, and radiant) over a 24-hour period. A major difficulty in modeling the phenomena occurring in such buildings is due to their complexity: the model must take into account heat transfer through the envelope and heat gain from the ground, air motion within the building due to forced and natural convection, vapor diffusion and condensation on the ice sheet, heat transfer by radiation between all internal surfaces, heat transfer by convection between the indoor air and all internal surfaces, sensible and latent heat gains from the spectators, conduction in the ice and floor as well as heat generation by the lights, the resurfacing operations, and the secondary coolant pump work. A second difficulty in modeling the interrelated phenomena occurring in this type of building is due to their size, irregular form, and lack of internal partitions. These buildings cannot be readily subdivided into distinct, well-mixed zones surrounded by surfaces of uniform temperature and uniform long- and shortwave radiation (these are the major assumptions of the heat balance and radiant time series methods described in chapter 30 of the 2005 ASHRAE Handbook--Fundamentals). Furthermore, such zonal methods require values of the heat and mass transfer coefficients that are not easy to specify accurately, especially for zones that are not separated by solid partitions. The use of a computational fluid dynamics (CFD) code (Chen and Srebric 2002) for the simultaneous calculation of all unknown variables eliminates the uncertainty resulting from the assumptions necessary for the application of zonal methods. Such code does not require values of the heat and mass transfer coefficients and provides more detailed descriptions of the flow field. Specifically, CFD codes determine the spatial distribution of air velocity, air temperature, and absolute humidity over the ice sheet. They also provide the spatial distribution of the ice, wall, and ceiling temperatures as well as the corresponding distributions of the heat fluxes into the ice sheet and into the secondary coolant.

Nielsen et al. (1979) and Davidson (1989) were among the first to apply CFD methods for airflow studies in ventilated rooms. However, they ignored the interaction between the indoor and outdoor environments and neglected heat transfer by radiation between the interior surfaces of the room. Chen and Van der Kooi (1988) proposed a methodology using a CFD code to calculate time-independent airflow patterns and a separate load calculation code, which solved the algebraic room energy balance equation using transient temperature differences. Li (1992) extended the application of CFD codes in the field of HVAC by including heat conduction through the walls and radiation transfer between their internal surfaces. This formulation was then applied (Li et al, 1993) to study the radiative effects on flow and temperature fields in a room ventilated by displacement. Yuan et al. (1992) and Yang et al. (2000) used a commercial CFD code to study ventilation and air quality in gymnasia and skating rinks, respectively. Very recently, Bellache et al. (2005a) used the same commercial CFD code to analyze steady-state simultaneous heat and mass transfer in ice rinks. Their model was improved for this study 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 pipes under the ice.

[FIGURE 1 OMITTED]

MODELING OF THE PROBLEM

The characteristics of eight ice rinks in the Montreal area were analyzed in order to select one that would serve as the basis for the project. Among the criteria used for this selection were the following: closeness to the two-dimensional configuration requested by ASHRAE; agreement with the owner/operator for the installation of the instrumentation necessary for the validation of the model; representativeness of the ice rink architecture and its HVAC & R systems and operation/occupation schedule.

This analysis led to the selection of the Camilien Houde ice rink, which is approximately 64 m (210 ft) long. A schematic representation of its cross section is shown in Figure 1. Seven rows of stands run the length of the building on one side, and a narrow corridor encircles the ice surface. The dasher boards are extended with protective transparent barriers separating the corridor and stands from the ice surface. Natural gas infrared heaters provide spot heating of the players' and scorekeepers' benches, which are situated on the west side of the building symmetrically with respect to the plane at Z = 32 m (105 ft). The spectator area's infrared heaters are distributed regularly over the length of the building. They are activated by motion sensors and deactivated by a thermostatic control when the air temperature over the stands reaches 15[degrees]C (59[degrees]F). As the heaters are direct-fired, the combustion products are expelled through roof-mounted exhaust fans interlocked with the burner controls. Outside air mixed with return air provides ventilation through regularly spaced diffusers above the spectator stands when the propane-fueled resurfacing machine is operating. Return air grilles are situated at the north and south ends of the stands. Humidity control is provided by two mechanical dehumidifiers located in the roof-truss spaces at each end of the building. A 24-hour operating schedule, which specifies occupancy, lighting, and resurfacing as well as possible stoppages of the refrigeration and heating operations during unoccupied periods (e.g., between midnight and 5 a.m. depending on the next morning's activities), is selected and contributes to the transient character of the conditions in the building (see appendix). Calcium chloride brine is used as the secondary coolant to maintain the ice surface temperature at the desired value. It is supplied from a header located at the west end of the ice sheet and circulates in the concrete floor under the ice in 146 two-pass polyethylene 32 mm (1.25 in.) tubes on 89 mm (3.5 in.) spacing.

It follows from this description that the building and its heating, cooling, and ventilation systems are symmetrical with respect to the mid-plane at Z = 32 m (105 ft). Assuming that spectators are uniformly distributed throughout the length of the stands and that the differences of outdoor conditions due to solar radiation and wind direction on the north/south walls are negligible, it follows that in the x-y plane midway through the length of the building (Z = 32 m [105 ft]) there will be no velocity components and no heat or mass transfer in the Z direction. Therefore, the conditions in this plane approximate very closely the two-dimensionality required by ASHRAE.

The airflow within ice-rink buildings is turbulent, as evidenced by the studies mentioned earlier. Several turbulence models exist and several were tested during an early phase of the project. The detailed results of these tests are not presented here due to lack of space. The principal conclusion is that the standard k-[epsilon] model (Launder and Spalding 1974) was very stable and gave satisfactory predictions compared to experimental and numerical results from the literature. Chen (1995) has also compared the performance of five turbulence models in predicting natural, forced, and mixed convection in rooms. He used experimental data for their validation and concluded that the performance of the standard k-[epsilon] model is good. Therefore, this model was selected for the present investigation.

The partial differential equations modeling the air movement, heat transfer, and vapor diffusion for transient conditions can be written in the following general form:

[[partial derivative]([rho][phi])/[partial derivative]t] + div([rho][[??].V][phi] - [[GAMMA].sub.[phi]][nable][phi]) = [S.sub.[phi]] (1)

The appropriate values or expressions for [phi], [[GAMMA].sub.[phi]], and [S.sub.[phi]] in the region occupied by air are defined in Table 1. It should be noted that in this region, these seven partial differential equations are elliptic, nonlinear, and coupled, since the air density varies with both the temperature and the absolute humidity:

[rho] = [[rho].sub.o][??]1 + [beta](T - [T.sub.o]) + [beta]*(C - [C.sub.o])[??] (2)

On the other hand, within the solids in the calculation domain (ceiling, walls, protective barriers, ice, concrete floor, insulation, and ground) the velocity is zero. For the purpose of this study, it is also assumed that these solids are impermeable to vapor and do not include any heat sources. Therefore, the temperature is the only unknown quantity and can be calculated from the two-dimensional transient conduction equation.

The boundary conditions for this problem are as follows.

The external surfaces of the building are subject to variable, site-dependant meteorological conditions. Since the objective of the project is to obtain "typical" daily profiles of refrigeration load components rather than corresponding maximum values, it was decided to use average meteorological conditions. The hourly data for three North American sites with very different climatic conditions (Edmonton, Houston, and Pittsburgh) provided by the US Department of Energy (DOE 2005) were used to calculate average daily profiles of sol-air temperatures (ASHRAE 2005) for the roof and east and west walls for winter (December-February), spring (March-May), summer (June-August), and autumn (September-November). Selection of the data used for the calculations was based on the average daily dry-bulb and horizontal sol-air temperatures. Summer conditions in Houston were excluded since operation of ice rinks under such extreme conditions is unlikely. From the remaining 11 site/season combinations, spring in Pittsburgh was chosen as the base case because it constitutes the median of these combinations. We also retained the two extreme cases (winter in Edmonton and summer in Pittsburgh) for sensitivity studies. Figure 2 shows the daily profiles of the dry-bulb and sol-air temperatures on the roof for the base case (spring in Pittsburgh), which have been used for the calculation of all the results presented in this paper, except during the validation of the model. By combining radiation and convection fluxes in a single term, the energy balance for any element of the external surfaces is

[q.sub.cd, o] = [h.sub.o]A([T.sub.sol-air] - [T.sub.s]). (3)

For the purposes of this study, the value of [h.sub.o] was considered constant and equal to 17 W/[m.sup.2] x K (3 Btu/h x [ft.sup.2] x [degrees]F), as recommended by ASHRAE (2005).

The dry-bulb temperature and humidity of the ventilation air entering the building when resurfacing takes place are the corresponding outdoor values calculated as explained in the previous paragraph. The horizontal and vertical components of the ventilation air velocity at the supply diffusers are both equal to 0.27 m/s. The relative pressure at the return air grilles is assumed to be zero.

[FIGURE 2 OMITTED]

The ground at a depth of 2 m (6.5 ft) is assumed to have a constant temperature (4.9[degrees]C [40.8[degrees]F] for the base case). Heat fluxes across the underground surfaces of the computational domain are evaluated as proposed by ASHRAE (2005).

On the inside surfaces of the ice rink (walls, ceiling, stands, and ice) the velocity components, concentration gradient, and turbulent kinetic energy are equal to zero. The corresponding temperature is related to that of the air within the flow domain by an energy balance that takes into account conduction through the solid, convection to the inside air, and net radiation fluxes between the inside surfaces. Thus, the conductive heat flux for each element of the inside surfaces is

[q.sub.cd, i] = [q.sub.cv, i] + [q.sub.r, i] + [q.sub.s, i] + [q.sub.l, i]. (4)

The last two terms denote the heat flux due to resurfacing operations and the latent heat due to vapor condensation (mass transfer to ice sheet) and are only taken into account for the ice surface. The conductive and convective heat fluxes [q.sub.cd] and [q.sub.cv] are calculated by the CFD code, while the net radiative heat flux [q.sub.r] between surface j and all the N inside surfaces of the envelope is calculated using Gebhart's method. Its implementation in the CFD code has been described in a paper by Bellache et al. (2005b).

The total loads due to vapor condensation, resurfacing, and brine pump work dissipation are calculated assuming quasi-steady conditions. In particular, we have used the formulation from ASHRAE (2002) to calculate the resurfacing heat flow on the ice surface. Thus,

[q.sub.l, i] = [dot.m.sub.in]([C.sub.in] - [C.sub.out])[h.sub.ig] (5)

[q.sub.s, i] = [rho][V.sub.f][4.2([T.sub.f] - 0) + 334 + 2(0 - [T.sub.ice])] (6)

[q.sub.p, i] = [rho]gQ[H.sub.sys] (7)

[q.sub.l, i] and [q.sub.s, i] are attributed to the ice surface (cf. Equation 4), while [q.sub.p, i] is attributed to the brine. The heat produced by a resting occupant is about 100 W (ASHRAE 2005).

CALCULATION PROCEDURE

The numerical solution of the system of coupled partial differential equations has been carried out using the finite volume method, a staggered grid, and the hybrid scheme for discretization of the convection terms. The SIMPLE algorithm was used for the pressure correction, and the solution was obtained iteratively using the numerical code PHOENICS.

The calculation domain was first subdivided into zones corresponding to the physical discontinuities shown in Figure 1. In the present case there are 27 zones in the X direction and 30 zones in the Y direction. The number of grid points, or finite volumes, within each of the resulting 810 subdivisions of the domain has been varied to ensure that the results are independent of their number. The discretization grid was always non-uniform with higher density near the air inlets and the solid surfaces where gradients are higher. Numerical tests were carried out with 100 x 50 and 200 x 100 grid points in the X and Y directions, respectively. The number of iterations was also varied. Table 2 shows that all the results are essentially identical. Therefore, the results reported in the present paper were calculated with 100 x 50 grid points and 10,000 iterations. Analogous tests were conducted with different timesteps. The results presented here were all calculated with a timestep of 600 s. Finally, the simulation was run for two days with identical 24-hour meteorological conditions in order to ensure that the calculated "typical" daily profiles of the refrigeration loads were periodic. Under-relaxation was often necessary to achieve convergence, which was declared when the cumulative residuals for each of the conservation equations were less than [10.sup.-6].

The horizontal surface through the centerlines of the floor pipes is assumed to be isothermal. This assumption is based on the small center-to-center spacing of the floor pipes as well as on the small increase of brine temperature (~1[degrees]C [2[degrees]F]) between the supply and return headers. For the purposes of the present study, the temperature of this surface is assumed to be equal to the average brine temperature in the floor pipes.

When resurfacing takes place, the total heat reaching the brine generally exceeds the capacity of the refrigeration plant. During these time intervals, the energy balance for the brine can be approximated as follows:

M[C.sub.p][d[bar.T]/dt] = [summation over i]([q.sub.cd, i] + [q.sub.p, i] (8a)

The summation on the right-hand side includes all the conductive fluxes reaching the isothermal surface from above (ice side) and below (ground side). The heat capacity term on the left-hand side includes contributions from the brine and the concrete in the appropriate finite volumes.

When resurfacing is not taking place, the capacity of the refrigeration plant is greater than the total heat transferred to the brine. During these time intervals, the space average of the brine temperature is assumed to be constant. Therefore, the corresponding energy balance is

[dot.m][C.sub.p][DELTA]T = [summation over i]([q.sub.cd, i] + [q.sub.p, i]). (8b)

VALIDATION

The adopted model and calculation procedure have been validated in several steps. Initially, the selected k-[epsilon] two-equation model was successfully validated for turbulent natural and mixed convection in small enclosures by comparing steady-state velocity and temperature profiles with corresponding measured values (Bellache et al. 2005a). Secondly, predicted air, ice, and wall temperatures for an ice rink were successfully compared with corresponding measured values in a previous project (Bellache et al. 2005b). For the purposes of the present project, long-term transient measurements that serve for the validation of the prediction model and for the better understanding of physical processes affecting the cooling load of the ice sheet and of occupant comfort are being conducted in the selected ice rink. The details of this experimental project are presented in a separate paper (Ouzzane et al. 2006). Furthermore, additional measurements of air and surface temperatures have been collected using RTD sensors and an infrared temperature sensor. In the following paragraphs we compare the experimental values for February 24, 2005, corresponding to the weekday schedule of the operating system, with the corresponding numerical predictions.

Figure 3 shows the outdoor temperature evolution over this 24-hour period and compares measured and calculated air temperature evolutions at two positions in the ice rink. The predicted time evolution of the air temperatures for each position follows very closely the corresponding measured profile. This is a first indication that the model satisfactorily represents the transient behavior of the phenomena taking place in the rink. The small, systematic underestimation of these temperatures is attributed to the inaccuracy of the assumed initial conditions used for the simulation. Indeed, since a complete three-dimensional map of air temperatures is not available, the simulation was executed with a uniform initial air temperature throughout the rink. Second, the model correctly predicts that the air temperature over the stands (X = 28.2 m) is approximately 4[degrees]C (7[degrees]F) warmer than the corresponding value near the west wall (X = 0.6 m) due to the operation of the radiant heaters. However, it is worth noting that the measured temperature over the stands rarely reaches the setpoint (15[degrees]C or 59[degrees]F) at Y = 7.7 m.

Figure 4 shows the vertical profiles of calculated and measured air temperatures near the center of the ice surface (X = 14 m, Y < 5 m [10 ft]) at 10 a.m. Hot wire measurements indicate that hardly any air movement takes place in this zone. In the first 20 cm (8 in.) above the ice, the agreement between measured and calculated temperatures is excellent. In this region, the temperature increases linearly with height and the corresponding gradient is quite high. This first linear region is followed by an anomalous behavior (for Y between 0.2 and 0.5 m) and a second region of linear increase. The predictions underestimate both the extent of this second region of linear increase and the corresponding gradient. This discrepancy may be due to the initial conditions or to the assumption of a two-dimensional flow field for the calculations. Above this second linear region, the air temperature is independent of the height. These constant values are in good agreement with the corresponding temperatures in Figure 3. This is significant since the values in Figure 3 were measured with the instrumentation used to establish the long-term transient performance, while those in Figure 4 where obtained with the RTD sensors. Finally, it is important to note that temperature profiles exhibiting the same characteristics as those in Figure 4 have also been measured by Pennanen et al. (1997).

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

Calculated and measured ice surface temperature evolutions for this same winter day are shown in Figure 5. These measurements were performed with an infrared temperature sensor. Again the transient response is very well predicted by the simulation. The systematic small difference between predicted and measured values is again attributed to the assumed initial conditions used for the simulation. The time evolution of this variable can be easily explained based on the prevailing conditions in the rink. Between midnight and 3 a.m. the lights are off while the refrigeration system is operating; therefore, the temperature of the ice surface decreases. Then the compressors are turned off from 3 a.m. to 5:30 a.m. and the ice surface temperature increases. After 5:30 a.m., the compressors and some lights are turned on and the ice temperature decreases steadily until 9:30 a.m. The peak temperatures occurring at more or less regular intervals between 10 a.m. and 10 p.m. are due to the resurfacing operation. During this operation, the refrigeration system is operated at maximum capacity to meet the peak load conditions. The difference between the measured and calculated peak temperatures is due to the fact that the measurement does not always coincide with the application of the hot water used for resurfacing purposes.

Further transient measurements and comparisons of calculated and measured values will become available during the second year (2005-06) of the project. However, since the comparisons in Figures 3, 4, and 5 (as well as other comparisons not presented here for lack of space) are quite encouraging, we have started using the model and the numerical code to calculate heat fluxes into the ice and the brine. Some preliminary results are presented in the next section of the paper.

[FIGURE 5 OMITTED]

RESULTS AND DISCUSSION

The ASHRAE request for proposals specified that the two-dimensional transient numerical model simulating ice sheet heat transfer should be used to determine how the cooling load is affected by climatic conditions as well as by changes in design and operating characteristics. After consultation with the members of ASHRAE TC 10.2, the three site/season combinations previously specified were selected to analyze the effects of the climate. The chosen design characteristics include two subfloor insulation thicknesses (0 and 7.6 cm [3 in.]), subfloor heating (yes or no), and three ceiling emissivities (0.9, 0.5, and 0.1). The operating characteristics include three ice thicknesses (2.5, 5.1, and 7.6 cm [1, 2, and 3 in.]), two operating schedules (see appendix), and three average brine temperatures (-7[degrees]C, -9[degrees]C, and -11[degrees]C [19.4[degrees]F, 15.8[degrees]F, and 12.2[degrees]F). The following paragraphs present calculated results for several different combinations of these parameters.

Base Case

The base case corresponds to the configuration of Figure 1 and the meteorological conditions of Figure 2 (spring in Pittsburgh). For this case there is no subfloor heating, the ice and subfloor insulation thicknesses are, respectively, 5.1 cm (2 in.) and 7.6 cm (3 in.), the average brine temperature in the floor pipes is -9[degrees]C (15.8[degrees]F), the ceiling emissivity is 0.9, and the weekend operating schedule (see appendix) applies.

Figure 6 shows the evolution of the temperature at different positions in the calculation domain over two 24-hour periods with identical meteorological conditions (see Figure 2). The first few hours must be disregarded because of the influence of the arbitrary initial conditions. Steady periodic conditions prevail for the second 24-hour period and, therefore, all subsequent discussion refers to the time interval between t = 24 hours and t = 48 hours. During that period, the temperature of the ground-insulation interface is constant. This is due to the large thermal inertia of the ground and the fact that the ground temperature at a depth of 2 m (6.5 ft) was assumed to be constant. Since its predicted value is just below the freezing point, it is probably advisable to install ground heating in this case to avoid freezing and heaving.

[FIGURE 6 OMITTED]

On the other hand, the temperature of the ice surface decreases slightly between midnight and approximately 8 a.m. since occupation, lights, and resurfacing are nonexistent during that period (c.f. appendix) while the brine temperature is constant. After 8 a.m., however, when activities pick up, the temperature of the ice surface varies considerably. Its peaks, which reach the freezing point, coincide with the resurfacing operations (c.f. appendix), which are assumed to last ten minutes each. After each of the morning resurfacing operations, the temperature of the ice surface returns to a more or less constant value of approximately -4[degrees]C (24.8[degrees]F). During the afternoon, however, when the frequency of the resurfacing operations increases, the temperature of the ice surface after each resurfacing tends to increase. Thus, at 9 p.m. (t = 45 hours), this temperature is approximately -2.5[degrees]C (27.5[degrees]F). This is due to the cumulative effect of the resurfacing operations and the approximately constant cooling capacity of the brine. The effect of resurfacing also influences the average peaks. The dew point of the air near the ceiling remains constant, and it is important to note that it is considerably lower than the ceiling temperature. Thus, there is no danger of condensation on the ceiling for the conditions under consideration. Finally, the ceiling temperature follows fairly closely the variation of the outside temperature (c.f. Figure 2).

Figure 7 shows the corresponding heating load profiles over the same period. Curve #1 shows that the heat flux from the ground to the brine is time independent. This is consistent with the results of the previous figure regarding the temperatures of the brine and the interface between the ground and the insulation. Curve #2 corresponds to the heat flux into the ice due to condensation. It is constant since the ice temperature and the state of the ventilation air change little during a typical day. Curve #3 shows the heat flux into the ice due to convection from the air in the ice rink. This quantity is nearly constant most of the time, except when resurfacing takes place. During these time intervals, the ice temperature increases (c.f. Figure 6); therefore, the convective heat flux decreases. Curve #4 depicts the radiation flux into the ice. It is lower during the night since the lights are off and the ceiling temperature is relatively low (see Figure 6). It becomes more important during the day when the effect of these two sources increases. The short-duration, steep decreases of this flux occurring between 8 a.m. (t = 32) and 11 p.m. (t = 47) coincide with resurfacing operations and are due to the corresponding increase of the ice temperature (c.f. Figure 6). Curve #5 shows the total heat flux into the ice from above (convection + condensation + radiation + resurfacing). The pronounced peaks coincide with the resurfacing operations (c.f. appendix). The maximum heat flux into the ice occurs around 2 p.m. (t = 38) and is approximately equal to 300 W/[m.sup.2] (95.13 Btu/h x [ft.sup.2]). Finally, curve #6 shows the total heat flux into the brine from above. For the entire day, the integral of this quantity is the same as the corresponding integral for curve #5. However, the inertia of the ice sheet and the concrete floor are such that the peaks of curve #6 are considerably lower than those of curve #5. For the same reason, they occur slightly later and last longer than those on the ice surface, which is directly exposed to the hot water applied during resurfacing. Thus, the maximum heat flux into the brine from above (this quantity does not include the heat flux from the ground) occurs at about 10 p.m. (t = 46) and is approximately equal to 150 W/[m.sup.2] (47.56 Btu/h x [ft.sup.2]), i.e., half the corresponding value for the ice.

[FIGURE 7 OMITTED]

Table 3 shows the daily loads into the ice (integral of the quantities in Figure 7 from t = 24 to t = 48) and the corresponding average fluxes for the base case defined previously. For these conditions, radiation is the largest contributor to the total refrigeration load, followed by approximately equal contributions from convection and resurfacing operations. It should be noted that the daily load into the ice (integral of curve #6) differs by only 0.5% from the corresponding value for the ice (line 6 of Table 3). This excellent result is a good indicator of the precision of the numerical calculations. The total pump work corresponds to 0.248 kWh/[m.sup.2] (78.62 Btu/[ft.sup.2]), while the friction dissipation in the floor pipes represents 0.076 kWh/[m.sup.2] (24.09 Btu/[ft.sup.2]) on a daily basis. The former is approximately the same as the contribution of condensation, while the latter is very close to the energy transmitted to the brine from the ground (0.082 kWh/[m.sup.2] [25.99 Btu/[ft.sup.2]]). Depending on the fraction of the total pump work that is lost to the surroundings, the total refrigeration load transmitted to the chillers is between 2.41 and 2.58 kWh/[m.sup.2] (763.97 and 817.86 Btu/[ft.sup.2]). The corresponding average fluxes are, respectively, 101 and 108 W/[m.sup.2] (32.02 and 34.25 Btu/h x [ft.sup.2]).

Parametric Studies

Figures 8 through 11 show the effects of different parameters on the daily loads into the ice (by convection, radiation, condensation, and resurfacing operations) and into the brine (from the ground and due to dissipation in the floor pipes) as well as on their sum. In each case, one clearly identified parameter is varied while the others all retain the values specified in the description of the base case.

Figure 8 shows that the differences between weekends and weekdays are essentially due to the lower number of resurfacing operations in the latter case. The corresponding load is significantly lower (0.352 kWh/[m.sup.2] [111.58 Btu/[ft.sup.2]]) and compensates for the small increases of the convection and radiation loads, which are due to a slight decrease of the average ice temperature. The latter is a direct result of the reduction in the daily number of resurfacing operations (9 on weekdays, 12 on weekends). The total is 2.41kWh/[m.sup.2] (763.97 Btu/[ft.sup.2]) on weekends and 2.34 kWh/[m.sup.2] (741.78 Btu/[ft.sup.2]) on weekdays. It is important to note that the load into the brine from the ground is totally unaffected by the operation schedule.

Figure 9 illustrates the effects of the climate by comparing the base case (spring in Pittsburgh) with winter conditions in Edmonton and summer conditions in Pittsburgh. Three of the four daily loads into the ice are lowest in winter and highest in summer, due to the corresponding changes of the outdoor dry-bulb temperature and absolute humidity. The influence of the climate on ground gains is very small, while daily loads from resurfacing and dissipation are independent of climatic conditions. The total is 1.98 kWh/[m.sup.2] (627.66 Btu/[ft.sup.2]) for winter in Edmonton, 2.41 kWh/[m.sup.2] (763.97 Btu/[ft.sup.2]) for spring in Pittsburgh, and 2.89 kWh/[m.sup.2] (916.13 Btu/[ft.sup.2]) for summer in Pittsburgh.

Figure 10 shows the effects of ceiling emissivity. The daily radiation load into the ice decreases considerably as the ceiling emissivity decreases. On the other hand, the corresponding load due to convection increases slightly due to a small decrease of the temperature of the ice surface. However, the effect on the total is very positive: this value is 2.41, 2.13, and 1.80 kWh/[m.sup.2] (763.97, 675.21, and 570.60 Btu/[ft.sup.2]), respectively, for [[epsilon].sub.c] = 0.9, 0.5, and 0.1. It is important to note that the load into the brine from the ground is totally unaffected by the ceiling emissivity.

Finally, Figure 11 shows that, as expected, the convective, radiative, and ground loads increase significantly as the average brine temperature decreases. The condensation load varies very slightly, while resurfacing and dissipation loads are not dependent on this parameter. The total is 2.25, 2.41, and 2.53 kWh/[m.sup.2] (713.25, 763.97, and 802.01 Btu/[ft.sup.2]), respectively, for average brine temperatures of -7[degrees]C, -9[degrees]C, and -11[degrees]C (19.4[degrees]F, 15.8[degrees]F, and 12.2[degrees]F).

[FIGURE 8 OMITTED]

[FIGURE 10 OMITTED]

CONCLUSION

Turbulent flow with simultaneous mass diffusion in a ventilated ice rink coupled to conduction through its envelope and radiation between the inside surfaces of the latter has been modeled using the standard k-[epsilon] model and a modified version of Gebhart's method.

The model and numerical code used for the solution of the coupled partial differential equations were validated by comparing their predictions with measured values in a municipal ice rink situated in Montreal. They were then used to evaluate the daily profiles of heat fluxes into the ice and the brine for a representative ice rink and steady periodic meteorological conditions. The results for the base case (spring conditions in Pittsburgh) show that the maximum heat flux into the ice occurs during resurfacing and is approximately equal to 300 W/[m.sup.2] (95.13 Btu/h x [ft.sup.2]). The corresponding value for the brine is approximately 150 W/[m.sup.2] (47.56 Btu/h x [ft.sup.2]). The daily loads into the ice were also calculated: for the base case, approximately 43% of the total is due to radiation, while convection and resurfacing each contribute approximately 19.6% of the total.

[FIGURE 9 OMITTED]

[FIGURE 11 OMITTED]

The effects of the operation schedule, climate, ceiling emissivity, and average brine temperature on the daily refrigeration loads were also evaluated and presented.

For the conditions under investigation, the results show that a decrease in the ceiling emissivity from 0.9 to 0.1 reduces the radiation energy portion of the total refrigeration load from 43% to 21.7%. This results in a reduction of 25% of the total ice load.

ACKNOWLEDGMENTS

The authors acknowledge the financial support received from ASHRAE for the research project and the technical support of ASHRAE TC 10.2, Automatic Ice Making Plants/ Skating Rinks.

NOMENCLATURE 

A = area, [m.sup.2] ([ft.sup.2])

C = absolute humidity, kg water/kg dry air (lb/lb)

[C.sub.p] = specific heat, kJ/kg x K (Btu/lb x [degrees]F)

[h.sub.ig] = specific enthalpy of sublimation, J/kg (Btu/lb)

[H.sub.sys] = hydraulic head, m (ft)

k = turbulent kinetic energy, [m.sup.2]/[s.sup.2] ([ft.sup.2]/[h.sup.2])

m = mass flow rate, kg/s (lb/h)

M = mass, kg (lb)

[q.sub.cd] = conductive heat flow rate, W (Btu/h)

[q.sub.cv] = convective heat flow rate, W (Btu/h)

q = heat flow rate due to condensation, W (Btu/h)

[q.sub.r] = radiant heat flow rate, W (Btu/h)

[q.sub.s] = heat flow rate due to resurfacing, W (Btu/h)

[q.sub.p] = heat flow rate due to pumping, W (Btu/h)

Q = volume flow, L/s (cfm)

[S.sub.[phi]] = source term in Equation 1

T = temperature, [degrees]C ([degrees]F)

t = time, s (h)

[bar.V]= velocity vector, m/s (ft/h)

V = volume, [m.sup.3] ([ft.sup.3])

X, Y = Cartesian coordinates, m (ft)

[beta] = thermal expansion coefficient, [K.sup.-1] ([R.sup.-1])

[beta]* = coefficient of mass expansion

[[GAMMA].sub.[phi]] = diffusion coefficient in Equation 1

[DELTA]T = temperature difference, K (R)

[epsilon] = dissipation rate of turbulent kinetic energy, [m.sup.2]/[s.sup.3] ([ft.sup.2]/[h.sup.3])

[[epsilon].sub.c] = ceiling emissivity

[lambda] = thermal conductivity, W/m x K (Btu/h x ft x [degrees]F)

[rho] = density, kg/[m.sup.3] (lb/[ft.sup.3])

[phi] = field variable in Equation 1

Indices

i = inside

o = outside

b = brine

f = flood water

gr = ground

p = pump

REFERENCES

ASHRAE. 2002. 2002 ASHRAE Handbook--Refrigeration, chapter 34. Atlanta: American Society of Heating, Refrigerating and Air-Conditioning Engineers, Inc.

ASHRAE. 2005. 2005 ASHRAE Handbook--Fundamentals. Atlanta: American Society of Heating, Refrigerating, and Air-Conditioning Engineers, Inc.

Bellache, O., M. Ouzzane, and N. Galanis. 2005a. Numerical prediction of ventilation patterns and thermal processes in ice rinks. Building and Environment 40:417-26.

Bellache, O., M. Ouzzane, and N. Galanis. 2005b. Coupled conduction, convection, radiation heat transfer with simultaneous mass transfer in ice rinks. Numerical Heat Transfer, Part A 48:219-38.

Chen, Q. 1995. Comparison of different k-[epsilon] models for indoor air flow computations. Numerical Heat Transfer, Part B 28:353-69.

Chen, Q., and J. Srebric. 2002. A procedure for verification, validation, and reporting of indoor environment CFD analyses. HVAC & R Research 8(2):201-16.

Chen, Q., and J. Van der Kooi. 1988. Accuracy--A program for combined problems of energy analysis, indoor airflow and air quality. ASHRAE Transactions 94(2):196-214.

Davidson, L. 1989. Ventilation by displacement in a three-dimensional room: A numerical study. Building and Environment 24(4):363-72.

DOE. 2005. ENERGYPLUS. Washington, DC: US Department of Energy. www.eere.energygov/buildings.

Launder, B.E., and D.B. Spalding. 1974. The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Energy 3:269-89.

Lavoie, M., R. Sunye, and D. Giguere. 2000. Potentiel d'economies d'energie en refrigeration dans les arenas du Quebec. Report prepared by the CANMET Energy Technology Center, Varennes, Canada.

Li, Y. 1992. Simulation of flow and heat transfer in ventilated rooms. Transactions of the Dept. of Mechanics/ Applied Comp. Fluid Dynamics, Royal Institute of Technology, Stockholm, Sweden.

Li, Y., L. Fuchs, and M. Sandberg. 1993. Numerical prediction of airflow and heat radiation interaction in a room with displacement ventilation. Energy and Buildings 20(1):27-43.

MacCracken, M.M., and J.P. Scott. 2003. Work statement: Develop and verify methods for determining ice sheet cooling loads. ASHRAE TC 10.2. Atlanta: American Society of Heating, Refrigerating and Air-Conditioning Engineers, Inc.

Nielsen, P.V., A. Restivo, and J.H. Whitelow. 1979. Buoyancy affected flows in ventilated rooms. Numerical Heat Transfer 2:115-27.

Ouzzane, M., R. Sunye, R. Zmeureanu, D. Giguere, J.P. Scott, and O. Bellache. 2006. Cooling load and environmental measurements in a Canadian indoor ice rink (RP-1289). ASHRAE Annual Meeting, Quebec City, Quebec, June 24-28.

Pennanen, A.S., R.O. Salonen, S. Alm, M.J. Jantunen, and P. Pasanen. 1997. Characterization of air quality problems in five Finnish indoor ice arenas. J. of the Air & Waste Management Association 47(10):1079-86.

Yang, C., P. Demokritou, and Q. Chen. 2000. Ventilation and air quality in indoor ice skating arenas. ASHRAE Transactions 160(2):339-46.

Yuan, X., Q. Chen, A. Moser, and P. Suter. 1992. Numerical simulation of air flows in gymnasia. Indoor Environment 1:224-33.

APPENDIX

Operating Schedule

             Lighting  Lighting  Resurfacing
Time (hour)  Weekday   Weekend   Weekday

 1            10%       10%        0%
 2            10%       10%        0%
 3            10%       10%        0%
 4            10%       10%        0%
 5            10%       10%        0%
 6            10%       10%        0%
 7            10%      100%        0%
 8            10%      100%        0%
 9            10%      100%        0%
10            10%      100%        0%
11            10%      100%        0%
12            10%      100%      100%
13            10%      100%        0%
14           100%      100%      100%
15           100%      100%        0%
16           100%      100%      100%
17           100%      100%        0%
18           100%      100%      100%
19           100%      100%      100%
20           100%      100%      100%
21           100%      100%      100%
22           100%      100%      100%
23           100%      100%        0%
24           100%      100%      100%

             Resurfacing  Occupancy  Occupancy
Time (hour)  Weekend      Weekday    Weekend

 1             0%          0%         0%
 2             0%          0%         0%
 3             0%          0%         0%
 4             0%          0%         0%
 5             0%          0%         0%
 6             0%          0%         0%
 7             0%          0%         0%
 8           100%         10%        20%
 9           100%         10%        20%
10             0%         10%        20%
11           100%         10%        20%
12           100%         10%        20%
13             0%         10%        20%
14           100%         20%        20%
15             0%         20%        20%
16           100%         20%        20%
17           100%         20%        50%
18           100%         20%        50%
19           100%         20%        50%
20           100%         20%        50%
21           100%         20%        20%
22             0%         20%        20%
23           100%         20%        20%
24             0%         20%        20%

The total power input due to lighting is 28 kW (95.53 MBtu/h) and the maximum number of spectators is fixed at 200. The heat generated by the spectators is taken equal to 100 W (341.2 Btu/h) (ASHRAE 2005) and is assumed to be uniformly distributed over the stands. The resurfacing heat load is calculated according to chapter 34 of the 2002 ASHRAE Handbook--Refrigeration (2002) and is assumed to be uniformly distributed over the ice surface and over a time interval of ten minutes, resulting in a heat flux of 236 W/[m.sup.2] (74.8 Btu/h x [ft.sup.2]).

O. Bellache, PhD

N. Galanis, PhD

M. Ouzzane, PhD

R. Sunye, PhD

Member ASHRAE

D. Giguere

O. Bellache and M. Ouzzane are research scientists, R. Sunye is a project manager, and D. Giguere is a technical expert at the CANMET Energy Technology Centre-Varennes, Quebec, Canada. N. Galanis is a professor at the Department of Mechanical Engineering, Universite de Sherbrooke, Sherbrooke, Quebec, Canada.

Table 1. Values of [phi], [[GAMMA].sub.[phi]], and [S.sub.[phi]]
in Equation 1

Equation           [phi]      [[GAMMA].sub.[phi]]

Mass               1          0
X-momentum         U          [mu] + [[mu].sub.t]
Y-momentum         V          [mu] + [[mu].sub.t]
Energy             T          [mu]/Pr + [[mu].sub.t]/[Pr.sub.t]
Absolute humidity  C          ([mu] + [[mu].sub.t])/[[sigma].sub.c]
Turbulent kinetic  k          [[mu].sub.t]/[[sigma].sub.k]
  energy
Dissipation rate   [epsilon]  [[mu].sub.t]/[[sigma].sub.[epsilon]]

Equation           [S.sub.[phi]]

Mass               0
X-momentum         -[partial derivative]P/[partial derivative]x
Y-momentum         -[partial derivative]P/[partial derivative]y -
                     [rho]g[beta](T - [T.sub.0]) - [rho]g[beta]*(C -
                     [C.sub.0])
Energy             [S.sub.T]
Absolute humidity  [S.sub.c]
Turbulent kinetic  G - [rho][epsilon] + [G.sub.B]
energy
Dissipation rate   [[epsilon]([c.sub.[epsilon]1]G -
                     [c.sub.[epsilon]2][rho][epsilon])/k] +
                     [c.sub.[epsilon]3][G.sub.B]([epsilon]/k)

[[mu].sub.t] = [rho][c.sub.[mu]][c.sub.d][k.sup.2]/[epsilon];
  [Pr.sub.t] = 1
[G.sub.B] = -g[beta]([[mu].sub.t]/[Pr.sub.t])([partial derivative]T/
[partial derivative][x.sub.i]); G =
[[mu].sub.t]([partial derivative][U.sub.i]/
[partial derivative][x.sub.j] + [partial derivative][U.sub.j]/
[partial derivative][x.sub.i]) x ([partial derivative][U.sub.i]/
[partial derivative][x.sub.j])
[c.sub.[epsilon]1] = 1.44; [c.sub.[epsilon]2] = 1.92;
[c.sub.[epsilon]3] = 1.44; [c.sub.[mu]] = 0.09; [c.sub.d] = 0.1643;
[[sigma].sub.k] = 1; [[sigma].sub.[epsilon]] = 1.3;
[[sigma].sub.c] = 1

Table 2. Effect of Grid Points and Iterations on Calculated Results

                Mass Balance                    Energy Balance
                [summation]([dot.m.sub.ou] -    [summation]([E.sub.ou] -
                [dot.m.sub.in])/[dot.m.sub.in]  [E.sub.in])/[E.sub.in]

Grid 100 x 50   0                               0.02
  [10.sup.4]
    iterations
Grid 200 x 100  0                               0.01
  [10.sup.4]
    iterations
Grid 100 x 50   0                               0.02
  [20.sup.4]
    iterations
Grid 200 x 100  0                               0.01
  [20.sup.4]
  iterations

                Radiation Load,       Condensation Load,
                W/[m.sup.2]           W/[m.sup.2]
                (Btu/h x [ft.sup.2])  (Btu/h x [ft.sup.2])

Grid 100 x 50    47.7                 11.8
  [10.sup.4]    (15.12)               (3.74)
    iterations
Grid 200 x 100   47.8                 11.8
  [10.sup.4]    (15.15)               (3.74)
    iterations
Grid 100 x 50    47.7                 11.8
  [20.sup.4]    (15.12)               (3.74)
    iterations
Grid 200 x 100   47.8                 11.8
  [20.sup.4]    (15.15)               (3.74)
    iterations

                Convection Load,      Ground Load,
                W/[m.sup.2]           W/[m.sup.2]
                (Btu/h x [ft.sup.2])  (Btu/h x [f.sup.2])

Grid 100 x 50   19.2                   2.7
  [10.sup.4]    (6.08)                (0.86)
    iterations
Grid 200 x 100  19.3                   2.8
  [10.sup.4]    (6.12)                (0.89)
    iterations
Grid 100 x 50   19.3                   2.7
  [20.sup.4]    (6.12)                (0.86)
    iterations
Grid 200 x 100  19.3                   2.8
  [20.sup.4]    (6.12)                (0.89)
    iterations

                Refrigeration
                Load,
                W/[m.sup.2]
                (Btu/h x [ft.sup.2])

Grid 100 x 50    81.4
  [10.sup.4]    (25.81)
    iterations
Grid 200 x 100   81.7
  [10.sup.4]    (25.91)
    iterations
Grid 100 x 50    81.5
  [20.sup.4]    (25.84)
    iterations
Grid 200 x 100   81.7
  [20.sup.4]    (25.91)
    iterations

Table 3. Ice Sheet Loads, Base Case (#1)

               Daily Refrigeration Load        W/[m.sup.2]
               kWh/[m.sup.2] (Btu/[ft.sup.2])  (Btu/h x [ft.sup.2])

Convection     0.475 (150.574)                 19.8 (6.278)
Radiation      1.045 (331.262)                 43.5 (13.793)
Condensation   0.260 (82.419)                  10.8 (3.424)
Resurfacing    0.470 (149.989)                 19.6 (6.215)
Heat gain      2.250 (713.245)                 93.8 (29.743)
above the ice