# Simplified thermal model of spaces cooled with combined positive displacement ventilation and chilled ceiling system.

An improved plume-multilayer model was developed and validated to
represent thermal transport in enclosures conditioned by radiant cooling
and displacement ventilation systems. A novel approach was developed to
estimate wall plumes for non-isothermal surfaces using the similarity solution derived for power law representation of temperature difference
between the room air and the wall. The nonuniform wall plume flow rates
predicted by the model agreed well with flow rates in an enclosure produced by computational fluid dynamics (CFD) simulations. Wall plumes
associated with a nonuniformly heated wall with varying temperature
difference between the wall and the air were found to be substantially
higher than the plume predicted by isothermal wall correlation.

The wall plume model is integrated with a multilayer space thermal model to predict the stratification height in the space, the vertical distribution of wall and air temperatures as a function of the chilled ceiling temperature and space air supply conditions. The plume-multilayer model results were compared with results of three-dimensional CFD simulations using commercially available software. Three test cases were considered for the simulations at cooling loads of 40 W/[m.sup.2] (12.7 Btu/h x [ft.sup.2]), 67 W/[m.sup.2] (21.2 Btu/h x [ft.sup.2]), and 100 W/[m.sup.2] (31.7 Btu/h x [ft.sup.2]) and supply airflow rates per unit area of 22.5 (1.2 [ft.sup.3]/min x [ft.sup.2]), 30 (1.6 [ft.sup.3]/min x [ft.sup.2]), and 37.5 [m.sup.3]/h x [m.sup.2] (2.1 [ft.sup.3]/min x [ft.sup.2]), respectively. The vertical wall and average air temperatures for each layer agreed well with the results of the plume-multilayer model, showing maximum errors in values of average air temperature of 0.13[degrees]C (0.23[degrees]F), 0.37[degrees]C (0.66[degrees]F), and 0.3[degrees]C (0.54[degrees]F) for the low, medium, and high load cases, respectively. The simplified model accurately predicted the stratification height at a maximum error of [+ or -]0.05 m (0.16 ft) in the three test cases. The stratification height is overestimated by 35% if wall plumes are neglected in the plume-multilayer model.

INTRODUCTION

The energy efficiency of HVAC systems in buildings has a major impact on building energy consumption. One of the HVAC systems that has potential for energy savings is the combined displacement ventilation and chilled ceiling (DV/CC) system (Novoselac and Srebric 2002). In displacement ventilation (DV), air is supplied into the room at floor level, at a low velocity, and at a temperature slightly lower than the desired room temperature. Thus, the cooler air will displace the warmer room air that rises by buoyancy, creating two zones, as reported by Jiang et al. (1992) and Yuan et al. (1998). The lower occupied zone contains fresh, cool air, and the upper zone contains the warm, contaminated recirculation air, which is exhausted near the ceiling. Maintaining the interface between these two layers above the occupied zone is one of the many unique challenges of DV design (Yuan et al. 1998).

Vertical temperature variation and thermal draft are more critical in the design of DV systems. According to ASHRAE (2005), thermal comfort considerations limit the supply airflow rate (recommended velocity ~0.15 m/s [29.53 ft/min], maximum velocity 0.25 m/s [49.21 ft/min]) to avoid cold drafts, and the supply air temperature (minimum temperature limit is 18[degrees]C [64.4[degrees]F]). It also recommends that the maximum vertical temperature gradient in the room should be lower than 2.5[degrees]C/m (1.4[degrees]F/ft) if the percent dissatisfied is to be kept less than 5%. Yuan et al. (1998) indicated that the maximum removed load from the occupied zone by displacement systems is approximately 40.0 W/[m.sup.2] (12.7 Btu/h x [ft.sup.2]). The Behne (1999) design diagram sets the cooling load limit of DV when combined with a chilled ceiling to 100.0 W/[m.sup.2] (31.7 Btu/h x [ft.sup.2]) of floor area.

One of the primary challenges in modeling enclosures cooled by combined DV/CC systems is their inherent complexity and detail-oriented nature, which makes simulations computationally expensive and, hence, hinders their integration in building energy simulation programs for optimization and control purposes. The modeling of DV has followed several approaches ranging from full simulations using computational fluid dynamics (CFD) tools (Chen 1995; Rees et al. 2001) to experimentally based semi-empirical models (Skaaret 1998) to nodal models (Rees 1998) to plume two-zone models (Mundt 1996; Carrilho da Graca 2003).

CFD modeling can give reasonable accuracy when simulations are made by expert users and when the suitable boundary conditions are imposed during the simulation. Chen (1995) compared the numerical results of five different models and recommended that the renormalized group (RNG) k-[epsilon] model be used for predicting indoor airflows. Good agreement was also found when comparing CFD simulations and test chamber measurements, as reported by Rees et al. (2001). CFD modeling is the most accurate, but it is prohibitively expensive when studying multi-room ventilation geometries or transient performance over a period of days or months.

A simpler approach for DV modeling is to use experimentally based semi-empirical correlations developed for particular geometries based on simple and clear approximations (Skaaret 1998). These correlations are suitable for use in the design stage of DV systems, but they give reliable results only when the geometries under consideration are similar to those for which the correlations were developed. A DV nodal model reported by Rees (1998), which represents the space by a network of nodes upon which heat and mass conservation equations are applied, is more general than the semi-empirical approach and is computationally less expensive than the CFD approach. However, the interaction between the different nodes depends on precalculated airflow rates that should be obtained experimentally or by CFD simulations for the particular geometries under consideration.

The plume two-zone model is derived from fundamental equations of thermal transport from plumes induced by heat sources in the space. This approach eliminates the need to perform CFD simulations to calculate plume flow rates. Plume two-zone models take a physical approach to the problem, where the space is divided into two horizontal air layers and the mass flow in and out of each layer due to heat sources is calculated from the plume equations reported by Mundt (1996). Previous work on plume-multilayer wall models has either neglected wall plumes or used isothermal wall flow models, as in the work of Carrilho da Graca (2003), Mundt (1996), and Linden et al. (1990). The applicability of the plume two-zone model was studied by Carrilho da Graca (2003) in cases where plume convection was stronger than wall-driven natural convection. In his work, Carrilho da Graca (2003) developed a model that predicted three temperatures: floor-level temperature, occupied zone temperature, and a mixed-layer outflow temperature. The plumes were assumed to be of equal strength and located at the room floor. The wall plume convection was assumed negligible.

Wall plumes in hot weather conditions may not be neglected due to the high solar gain of the wall, which creates larger temperature differences between room air and wall surface. In many places with hot climates, building common practice uses concrete-based walls of high U-factor (U = 2.4 to 3.4 W/[m.sup.2] x K [0.4 to 0.6 Btu/h x [ft.sup.2] x [degrees]F]) that contributes to higher internal wall surface temperatures. The use of isothermal wall correlations in plume models in spaces where DV is combined with chilled ceilings is also incorrect. With DV, the nonuniformity in wall temperature is significant by the nature of the cooling mechanism. There is an apparent need for a modified model that deals with nonuniform wall conditions and heat sources at a low computational cost and without having to perform three-dimensional simulations to get the associated transport coefficients. A computationally efficient model that captures the physics of the problem would easily be incorporated with the HVAC system optimization and control program for reducing energy consumption under steady and dynamic conditions and for determining the building system efficiency. The use of a simple and accurate model is essential in finding means to improve the system response and energy consumption in response to changes to load within the enclosure.

For this study, a wall plume model was developed to accurately predict flow rates adjacent to a wall. The model solved the wall plume equations for two-dimensional boundary layer self-similar laminar flow while using a power law distribution for the temperature difference between the wall and the ambient conditions. The wall plume model was tested for accuracy when used in infinite and closed space domains. The wall-plume model was integrated with source plumes and a multilayer space model. The plume-multilayer model was validated by three-dimensional simulations using Airpak (2002) commercial software, and comparison was made with isothermal wall space conditions. The contribution of the work is the ability of the developed improved model to predict at low computational cost the DV/CC thermal transport, the stratification height, and the vertical air and wall temperature gradients.

MODEL MATHEMATICAL FORMULATION

The main feature of the plume-multilayer model is its ability to identify the fundamental physics of the problem in a simplified, computationally efficient approach. The modeling of the plumes generated by convective fluxes is a major part of the proposed model. The mass flow rates and the average temperature of the plumes are evaluated using expressions that take into account the strength of the convective flux and the position of the source as reported by Mundt (1996). Since the chilled ceiling temperature is relatively lower than the temperature of the other surfaces, the radiant heat exchange cannot be neglected. Unlike other models that assume conditions of adiabatic walls or isothermal walls, this model takes into account the nonuniformity of the wall temperature and the dependence of this temperature distribution on the internal and external conditions. The temperatures of the walls at different heights will be determined by applying the energy conservation equations to the wall surfaces. The model can predict the circulation (stratification) height. The convective and the radiative heat loads on the ceiling panel can also be predicted.

However, the model has limitations based on some of the assumptions. The heat sources are sufficiently distant to prevent plume interaction. Point heat sources are located at their actual positions, and heat sources with dimensions are located at equivalent virtual point source positions (Goodfellow and Esko 2002). It is also assumed that the equilibrium height of each plume is not less than three-fourths of the room height, where the equilibrium height is the elevation at which the density gradients disappear and the plume spreads horizontally, as reported by Goodfellow and Esko (2002).

Wall Plume Model

The stratification height and the temperatures of air and room walls are affected by the wall plumes. For an isothermal wall, the volume flow rate of the wall plume per unit width is given by Jaluria (1980) and Etheridge and Sandberg (1996) as

[Q.sub.w] =0.00287[DELTA][T.sup.1/4][Z.sup.3/4], (1)

where Z is the distance from the vertical wall leading edge and [DELTA]T is the temperature difference between the isothermal wall and the fluid far from the wall. For DV with a chilled ceiling system, the temperature differences between the top and the bottom of the wall are significant and, thus, the assumption of an isothermal wall may be overestimating or underestimating the wall plume flow rates. A more rigorous approach using the two-dimensional laminar boundary layer similarity analysis of Gebhart et al. (1988) in a stable, stratified environment was followed to determine the wall plumes. The two-dimensional transport boundary region governing equations of thermally induced motion are:

Continuity equation [[partial derivative]u/[partial derivative]x] + [[partial derivative]v/[partial derivative]Z] = 0 (2a)

Momentum equation u[[partial derivative]u/[partial derivative]x] + v[[partial derivative]v/[partial derivative]Z] = [v][[[partial derivative].sup.2]u/[partial derivative][Z.sup.2]] + g[beta](T - [T.sub.a])(2b)

Energy equation u[[partial derivative]T/[partial derivative]x] + v[[partial derivative]T/[partial derivative]Z] = a[[[partial derivative].sup.2]T/[partial derivative][Z.sup.2]] (2c)

subject to the following boundary conditions:

u(x,0) = 0

u(x,[infinity]) = 0

v(x,0) = 0

T(x,0) = [T.sub.w]

T(x,[infinity]) = [T.sub.a]

where

[T.sub.w] - [T.sub.a] = d(Z) (3)

where u is the horizontal component of the velocity, v is the vertical component, [upsilon] is the kinametic viscosity, [alpha] is the thermal diffusivity, g is the gravitational constant, and [beta] is the thermal expansion coefficient. The air and wall temperature difference decreases with height in spaces conditioned by a combined DV/CC system. A power law distribution represents well the temperature difference between the wall and the air d(Z) and is given by

d(Z) = [T.sub.w] - [T.sub.a] = N[Z.sup.n] (4)

where N is a constant. When n = 0, then an isothermal wall case is present if there is no variation in the air temperatures. The corresponding transformed boundary layer equations of momentum and energy derived by Gebhart et al. (1988) in terms of the transformation parameter [eta] are given in terms of a dimensionless stream function f([eta]) and dimensionless temperature [phi]([eta]) = (T - [T.sub.a])/d(Z) as follows:

f'" + (n + 3)ff" - (2n + 2)[f'.sup.2] + [phi] = 0 (5a)

[phi]" +Pr[(n + 3)f[phi]' - 4nf'[phi]] = 0 (5b)

subject to the following boundary conditions:

[[partial derivative]f/[partial derivative][eta]](0) = 0, [[partial derivative]f/[partial derivative][eta]]([infinity]) = 0, f(0) = 0, [phi]([infinity]) = 0, [phi](0) = 1 (5c)

where [eta] is defined as [eta](Z, y) = b(Z)y, and Pr is the Prandtl number. The function b(Z) is defined by

b(Z) = [1/Z]4[square root of ([Gr.sub.Z]/4)], (5d)

where [Gr.sub.Z] is the Grashof number ([Gr.sub.Z] = g[Z.sup.3]/[v.sup.2] [beta]d(Z)). The volume flow rate per unit width at any distance Z is calculated as

[Q.sub.w](Z) = 4v([Gr.sub.Z]/4)[.sup.1/4] [f.sub.[infinity]]. (6)

The volume flow rate per unit width of the wall plume [Q.sub.w] is directly related to the "shape" of the temperature distribution along the vertical wall rather than a computed average wall temperature. Gebhart et al. (1988) reported a range of the exponent n for which the solution is physical to be between -0.6 and unity (-3/5 < n < 1). The parameter [f.sub.[infinity]] depends only on the exponent n and is acquired from the numerical solution of ordinary differential Equations 5a and 5b using the shooting method (Charpa and Canale 2002) with first-order finite-difference discretization of f', f", f'", [phi]', and [phi]". The numerical solution is obtained for values of the exponent n between -0.6 and 1.0 in steps of 0.05 and increments of the independent variable [eta] at [DELTA][eta] = 0.0001. Figure 1 shows a plot of the value of the calculated parameter [f.sub.[infinity]] as a function of n in the desired range for Pr = 0.7. The solution for the special cases of the isothermal (n = 0) and constant heat flux (n = 0.2) wall presented by Jaluria (1980) are reproduced by the current general solution and are also shown in Figure 1. The [f.sub.[infinity]] - n relation is expressed as a third-order polynomial equation, with a correlation parameter value of 0.9998, as follows:

[FIGURE 1 OMITTED]

[f.sub.[infinity]] = 0.1145[n.sup.3] + 0.2271[n.sup.2] - 0.2978n + 0.6024 (7)

Since the real physical temperature difference distribution decreases with height in rooms conditioned with a DV/CC system, the maximum value of n should be less than zero. The range of exponent n for wall plumes is taken as -0.6 < n < 0.

The major challenge in the wall plume model is to develop a methodology to estimate the values of n and N that give accurate predictions of wall flow rates and yet provide a good power law curve fit of the calculated temperature difference between the wall and the air in the space.

Determination of the Power Law Parameters

The equations through which the mass flow rates are calculated depend on the power n and the coefficient N of the wall-to-air temperature difference distribution. For negative powers (n < 0), high temperature differences occur at low elevations and the difference is infinite at Z = 0. The similarity solution will result, then, in a relatively large mass flow rate at a low level, overestimating the real wall flow rates in the space. To overcome this problem, a low-height offset ([Z.sub.o]) of the wall is set to a uniform temperature equal to the local temperature at [Z.sub.o]. In practical cases, it can be assumed that the wall temperature is uniform at low levels. The flow rate of the plume at [Z.sub.o] resulting from the power law solution is matched to the flow rate produced at low height by the isothermal solution estimated from Equation 1. A valid curve fit for the wall-to-air temperature difference is then derived to give the highest correlation factor ([R.sup.2]) between the data values produced from the plume-multilayer model. The curve fit must also give the same flow rate given by the isothermal equation at [Z.sub.o] to offset the effect of the high flow rate in the region near Z = 0.

Thus, the calculation of the coefficient N and the power n used in determining the plume of a non-isothermal wall is subjected to two constraints. The first constraint is that N and n must give a flow rate close to the value given by the isothermal equation for low elevation at a selected height [Z.sub.o]. The second constraint is that the error between the corresponding temperature difference curve fit using the least squares method and the actual temperature difference is minimal or, alternatively, that the correlation parameter [R.sup.2] is maximal.

The method starts by determining the coefficient N using the least squares method for each power in the range of n = -0.6 and n = 0 in steps of [DELTA]n = 0.01 to generate the curve fits at these values of N and n. Then the corresponding wall plume flow rates are calculated from Equations 5 and 6 at the desired offset elevation [Z.sub.o]. These flow rates are then compared with the flow rate calculated at [Z.sub.o] from the isothermal equation (1), assuming the wall, up to the height [Z.sub.o], is at uniform temperature equal to the value of the wall temperature at [Z.sub.o]. The values of n and N that give the nearest flow rate to that of the isothermal equation with the highest [R.sup.2] value will be selected. The value of n could be refined by using smaller steps of [DELTA]n. The effect of the selection of the offset height [Z.sub.o] on the model results were tested at values of [Z.sub.o] = 0.15 m (0.49 ft), 0.3 m (0.98 ft), and 0.6 m (1.97 ft).

Testing the Wall Plume Model

The nonuniform thermal conditions of the wall used in the wall-plume model testing are a power law temperature difference distribution for Z > [Z.sub.o] and an isothermal wall at [T.sub.w]([Z.sub.o]) for Z [less then or equal to] [Z.sub.o]. The validation of the model is done by testing the accuracy of the followed procedure to determine the power law parameters that predict accurate flow rates.

The wall plume model methodology of selection of power law parameters and the accuracy of predicted flow rates are tested in the following ways:

* By making comparisons of the model wall flow rate predictions to results of Fluent (2004) CFD simulations of two-dimensional boundary layer flow adjacent to a nonuniformly heated wall in an infinite domain. The similarity solution from which the wall plume model was formulated is based on the assumption of an infinite domain. As a first step, the values of the flow rates calculated from the wall plume model are compared with flow rates calculated from infinite domain CFD tests. The aim of the infinite domain test is to check the suggested methodology for choosing the power n and the coefficient N that will be used in the flow rate calculation. The infinite domain test permits applying a power law temperature difference distribution while keeping the infinite domain air temperature constant. This condition is not as simple to apply in a closed domain, when the room vertical temperature variation is not known a priori. The second aim is to compare the flow rates from the CFD results with flow rates given by the isothermal equation.

* By making comparisons of model wall flow rate predictions to results of Fluent (2004) CFD simulations of two-dimensional flow in an enclosure with zero source plumes and in the presence of a chilled ceiling to isolate the effect of wall plumes due to nonuniformity of wall temperature. The enclosure test will also check the wall plume model applicability when both room air and room wall temperature distributions vary with height while maintaining a power law distribution of the temperature difference.

Infinite Domain CFD Test. Although many tests were performed to validate the model, the results are presented for only a few typical temperature difference distributions. The temperature difference distributions were based on the published experimental data of Rees (1998) and span the possible realistic range of temperature differences in rooms cooled by a combined DV/CC system. In infinite domains, the temperature difference between the wall and air is controlled by setting the air temperature at a constant and assigning the wall temperature distribution that satisfies the power law. The infinite domain test cases of natural convection of air at a uniform temperature of 23[degrees]C (73.15[degrees]C) adjacent to a vertically heated wall were simulated for a wall height of 4 m (13.12 ft) with seven wall temperature distributions, as shown in Figure 2. The first five temperature distributions had the maximum temperature difference at the lowest elevation less than 1.5[degrees]C (2.7[degrees]F). The cases FC6 and FC7 present limiting cases where the maximum temperature difference is equal to 0.72[degrees]C (1.3[degrees]F) and 2.8[degrees]C (5[degrees]F), respectively. The FC6 case is directly obtained from the experimental results of Rees (1998). The FC7 case was constrained by the condition that the flow at 3.0 m (9.8 ft) height remained laminar with the Rayleigh number below [10.sup.9]. The temperature difference distribution of case FC7 represented the maximum possible temperature difference at the lowest level that could maintain a laminar flow below a height of 3.0 m (9.8 ft).

The computational domain used in the Fluent (2004) simulations is extended in the horizontal direction to 2 m (6.56 ft). The heated wall temperatures were set such that for [Z.sub.o] [less then or equal to] Z [less then or equal to] 3.3 m (10.8 ft), the temperature difference between the wall and the air has a power law distribution. The heated wall temperature at Z [greater then or equal to] 3.3 m (10.8 ft) is set equal to that of the environment at 23[degrees]C (73.15[degrees]F) and for Z [less then or equal to] [Z.sub.o] the wall is isothermal at [T.sub.w]([Z.sub.o]) (e.g., case FC1: [T.sub.w]([Z.sub.o])=[T.sub.w] [0.3 m] = 24.5[degrees]C [75.8[degrees]F] for Z [less then or equal to] 0.3 m [0.98 ft]). The other boundaries of the domain were set at a constant temperature of 23[degrees]C (73.15[degrees]F), while the lower and upper boundaries were set as pressure inlet and outlet boundaries, respectively. The number of cells in the computational domain was 12,535. The 110 nodes in the vertical direction were evenly spaced, while the 127 nodes in the horizontal direction had variable distribution with a very fine grid close to the heated wall. The minimum and maximum values of the cell area were 1.135696 x [10.sup.-6] [m.sup.2] (1.22 x [10.sup.-5] [ft.sup.2]) and 3.965131 x [10.sup.-6] [m.sup.2] (4.268 x [10.sup.-5] [ft.sup.2]), while the minimum and maximum values of the cell face length were 4.5549 x [10.sup.-5] m (4.903 x [10.sup.-4] [ft.sup.2]) and 1.1172 x [10.sup.-1] m (1.2025 [ft.sup.2]). Convergence was reached after 1,159 iterations. Grid independence tests were also performed at a finer grid of 220 x 254 and showed no significant difference in velocity values. The velocity profiles are determined at different elevations. Figure 3 shows the variation of the vertical velocity component in the x direction at an elevation equal to 3.0 m (9.8 ft). The vertical velocity values are integrated over the x horizontal distance to obtain the wall flow rates.

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

Flow rates produced by Fluent (2004) software are compared with values obtained from the derived wall plume equations used in the multilayer model. Figure 4 shows the mass flow rates of the wall plumes calculated using the power law similarity solution for different fits as compared to the values calculated from two-dimensional CFD simulations for test case FC1 with offset height [Z.sub.o] = 0.3 m (0.98 ft). For the case FC1 distribution, the best power law parameters n and N resulting from applying the selection procedure of the model were found to be -0.23 and 0.5792, while the [R.sup.2] value was equal to 0.701. The mass flow rate calculated from the wall plume model at n = -0.23 agrees well with those calculated from the simulations. The maximum percentage error for n = -0.23 is 2.9%. The percentage error from the isothermal (average temperature, n = 0) case varied between 16% and 10%, while the percentage error for the best fit case of similarity wall plume model at n = -0.6 ([R.sup.2] = 1) was 45%. Similar results and accuracy of flow rates as those found for case FC1 were obtained for the other cases. Note that the wall plume flow rates below the offset height [Z.sub.o] (Z [less then or equal to] [Z.sub.o]) must be calculated based on the isothermal equation rather than the wall plume model.

[FIGURE 4 OMITTED]

Figure 5a shows the local Rayleigh number as a function of elevation for the FC6, FC7, and FC1 cases at [d.sub.max] = 0.72[degrees]C (1.3[degrees]F), 2.8[degrees]C (5[degrees]F), and 1.5[degrees]C (2.7[degrees]F), respectively. The isothermal region elevation is taken at 0.3 m (0.98 ft). The flow rates are calculated for the different fits of each case and compared to the values obtained from the isothermal equation. Results from the CFD simulation of each case are compared to the flow rates calculated from different fits using the wall plume model. The variations of the flow rate with elevation using the fit that matches the isothermal flow rate at 0.3 m (0.98 ft) elevation and using CFD predictions for cases FC6 and FC7 are shown in Figures 5b and 5c, respectively. The maximum error in flow rates is less than 10%. This indicates that the model is insensitive to the maximum temperature difference present in the distribution.

The effect of the elevation at which the comparison between the flow rates from the isothermal equation and the model is made (0.3 m in the above cases) was considered by selecting three different values of [Z.sub.o] for case FC8 at 0.15 m (0.49 ft), 0.3 m (0.98 ft), and 0.6 m (1.97 ft). The FC8 case is the same as the FC1 case but with the first 0.6 m of the heated wall set at a temperature of 24[degrees]C (75[degrees]F). The wall plume model was used to calculate the wall plume flow rate as a function of elevation. The optimal power n was -0.19, -0.21, and -0.25, and the corresponding optimal coefficients N were 0.48, 0.48, and 0.49 for the cases of 0.15 m (0.49 ft), 0.3 m (0.98 ft), and 0.6 m (1.97 ft) offset heights, respectively. Figure 6 shows the variation of wall plume model flow rates at optimal fits at 0.15 m (0.49 ft), 0.3 m (0.98 ft), and 0.6 m (1.97 ft) offset heights of case FC8. The maximum percentage difference in the predicted flow rate of the different fits is about 6.75% and occurs at a height of 0.15 m (0.49 ft). The percentage difference decreases as the elevation increases. This implies that the sensitivity of the model to the isothermal height introduces minimal error in the value of mass flow rates above [Z.sub.o].

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

The infinite domain test results show that the suggested methodology to calculate the flow rates using open-domain similarity solutions of adjusted power law parameters can accurately predict flow rates of a nonuniformly heated vertical wall with power law distribution with an isothermal initial starting length. The model reproduces the two-dimensional CFD flow rate predictions in an open domain with the same heated wall temperature distribution.

Two-Dimensional Enclosure CFD Test. The applicability of the wall plume model obtained from the infinite domain solution were tested in closed domains (enclosures) of typical room heights at zero internal heat sources. Two-dimensional simulations of thermally driven flow in a rectangular enclosure with nonuniformly heated walls in the presence of a chilled ceiling boundary were performed using CFD software (Fluent 2004). Unlike infinite domain, the CFD solution will generate room air temperature variations. These temperatures will then be used as input to the wall plume model to compare resulting flow rates with the two-dimensional enclosure flow rates produced by CFD and with flow rates based on isothermal flow equations. Note that when space dimensions are large, the wall plume flow rates obtained from CFD simulation will approach the values obtained from an infinite domain.

A symmetrical 4 m (13.1 ft) wide and 3 m (9.840 ft) high two-dimensional enclosure was simulated. Due to symmetry, only half of the physical domain was considered in the computational domain and the symmetry boundary was represented by an adiabatic wall. The upper boundary (ceiling) and the lower boundary (floor) of the rectangle are assigned constant surface temperature conditions of 18[degrees]C (64.4[degrees]F) and 23[degrees]C (75[degrees]F), respectively. Cases FE1, FE2, FE3, and FE4 are considered in an enclosure corresponding to various temperature difference distributions between the vertical wall and air with computed room air temperatures of 22.9[degrees]C (73.2[degrees]F), 22.5[degrees]C (72.5[degrees]F), 21.89[degrees]C (71.4[degrees]F), and 21.9[degrees]C (71.4[degrees]F), respectively, and corresponding maximum temperature differences of 4.9[degrees]C (8.8[degrees]F), 4.0[degrees]C (7.2[degrees]F), 2.6[degrees]C (4.7[degrees]F), and 2.3[degrees]C (4.1[degrees]F), respectively, as shown in Figure 7a. In Figure 7b, a variation of the vertical velocity component is plotted in the horizontal direction at an elevation of 2.7 m (8.86 ft) for case FE1. It is clear how the velocity profile in the boundary region is disturbed. When calculating the mass flow rate, the boundary layer region is considered to extend from the wall to the point where the velocity profile sharply changes slope.

[FIGURE 7 OMITTED]

Following the wall plume model methodology described to select the temperature difference distribution parameters for the best fit, the power n and the corresponding coefficient N are determined. For the four cases of FE1, FE2, FE3, and FE4, the optimal power law parameters (n, N) from the wall plume model were (-0.195, 2.526), (-0.26, 1.56), (-0.193, 1.64), and (-0.092, 1.67), respectively. The mass flow rates from the wall model were calculated using the obtained powers and corresponding coefficients and are compared with the mass flow rates obtained from CFD simulations. Figures 8a-8d show the variation of wall plume flow rates obtained from the model and from the two-dimensional enclosure simulations as a function of elevation for cases FE1 (a), FE2 (b), FE3 (c), and FE4 (d). The maximum percentage error between the predicted flow rates from the wall plume model and the two-dimensional enclosure simulations are less than 10% in all the four cases for elevations Z > [Z.sub.o]. The wall plume model predictions of flow rate are far more accurate than the isothermal wall predictions at the average temperature difference. The accuracy improves over the isothermal wall assumption with increased wall-to-air temperature difference, which is apparent when comparing case FE1 to case FE4 in Figures 8a and 8d.

The wall plume model has been shown to have an acceptable accuracy in predicting flow rates adjacent to nonuniformly heated walls in a two-dimensional enclosure with chilled ceiling. The next step is to implement the wall plume model in the space multilayer model to predict the DV/CC space vertical temperature gradients, stratification height, and chilled ceiling load.

Plume-Multilayer Model

The space is divided into four adjacent horizontal air layers with adjacent wall layers as shown schematically in Figures 9a and 9b. Each layer is represented by a uniform temperature T(i), where i corresponds to the air layer index (i = 1 to 4). Each layer interacts with side wall segments. The wall segment is represented by two nodal temperatures, the outer and the inner nodes. The thermal resistance diagram for the wall segment is shown in Figure 9c. The chilled ceiling temperature [T.sub.c] is constant and uniform. The fresh supply air enters from the side wall into the floor air layer at flow rate [M.sub.s] and temperature [T.sub.s] and the exhaust air leaves from the side wall of the upper air layer at [T.sub.(4)]. The wall temperature is designated by [T.sub.w](i,j), where j corresponds to the wall orientation (j = 1 to 4, where 1 = east, 2 = west, 3 = north, and 4 = south). The outside convection coefficients are [h.sub.ocv,f] and [h.sub.ocv,(i,j)] for the floor and walls, respectively.

The flow rates resulting from walls and heat sources are denoted by [M.sub.w(j,k)] and [M.sub.p(l,k)], respectively. The net circulated mass [M.sub.cir] at each boundary (interface surface between two adjacent air layers i and i + 1) will be calculated. The enclosure mass balance equation is given by

[M.sub.cir,k] = [M.sub.s] [r.summation over (l-1)] [M.sub.p,(l,k)] [4.summation over (j-1)] [M.sub.w,(l,k)], (8a)

where the wall plume mass flow rate of any wall j is given by

[M.sub.w,(j,k)] = [rho][Q.sub.w,(j,k)][Y.sub.(j)]. (8b)

In Equation 8, [rho] is the air density in kg/[m.sup.3] (lbm/[ft.sup.3]), [Y.sub.(j)] is the width of wall j, k refers to the interface level between the consecutive air layers (k = 1, 2, or 3 for a space divided into four layers), and r is the number of heat sources in the first air layer. If the circulation mass is positive, then the circulated mass is moving upward; if the circulated mass is negative, then it is moving downward. The determination of this circulated mass is very important because the stratification height is defined at the level where the circulated mass is zero.

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

The energy balance is then written for each horizontal air layer and for each wall segment. The first air layer energy balance is given by

[r.summation over (l=1)][[PHI].sub.(l,1)] + [4.summation over (i=1)] [h.sub.cv,(i,1)][A.sub.w,(i,1)] ([T.sub.w,(i,1)] - [T.sub.a,(1)]) + [h.sub.cv,f]([T.sub.f] - [T.sub.a,(1)]) + [M.sub.s][H.sub.s] = [r.summation over (j=1)] [M.sub.p,(j,1)][H.sub.p,(j,1)] + [4.summation over (i=1)] [M.sub.w,(i,1)][H.sub.w,(i,1)] + [M.sub.cir,1][H.sub.cir,1], (9)

where [PHI] is the heat source flux and H is the enthalpy in J/kg (Btu/lbm). For the left-hand side of Equation 9, the first term is the summation of the heat fluxes in the first layer, the second and third terms are heat convection from the walls and the floor, and the fourth term is the enthalpy of the supply air. For the right-hand side of Equation 9, the first and second terms are the enthalpies of the heat sources and wall plumes, and the third term is the enthalpy of the circulated mass at the boundary. The adopted correlations for laminar and turbulent convective heat transfer coefficients for vertical walls are, respectively:

[h.sub.cv] = 1.42 ([DELTA]T/Z)[.sup.0.25] (laminar convection with Z dimension in meters) (10a)

[h.sub.cv] = 1.31([DELTA]T)[.sup.0.33] (turbulent convection) (10b)

while the adopted correlation for the floor or ceiling is that of Min et al. (1956)given by

[h.sub.cv,f] = 2.13([DELTA]T)[.sup.0.31] (10c)

The value used by Rees (1998) was 5.9 W/[m.sup.2] x K (1.04 Btu/h x [ft.sup.2] x [degrees]F) and is close to the correlation of Min et al. (1956). The energy balance equations on the two intermediate layers and the upper layer can be similarly derived. It is also worthy to note that the number of horizontal layers can be increased from four to any number without any restrictions. This will make the model more powerful when dealing with more complicated geometries (presence of windows) and furniture settings.

The energy balance on the outer wall node gives

[h.sub.cv,o]([T.sub.out] - [T.sub.so]) + [q.sub.solar] = [U.sub.w]([T.sub.so(i,j)] [T.sub.si,(i,j)]), (11)

where [U.sub.w] is the wall conductance and [q.sub.solar] is the solar flux into the wall. The energy balance on the wall inner node yields

[U.sub.w]([T.sub.so(i,j)] - [T.sub.si,(i,j)]) +[q.sub.l] = [q.sub.cv] + [q.sub.rd], (12)

where [q.sub.l] is the radiative heat flux from light loads, [q.sub.cv] is the convective heat flux, and [q.sub.rd] is the radiative heat flux. The convective heat flux [q.sub.cv] on each part of the wall is given by

[q.sub.cv] = [h.sub.cv]([T.sub.si] - [T.sub.a]). (13)

The radiative heat flux [q.sub.rd,i] on a surface of index i is given by

[q.sub.rd,i] = [J.summation over (j=1)][epsilon][sigma]([T.sub.i.sup.4] - [T.sub.j.sup.4] [F.sub.i] [right arrow]j, (14)

where J is the number of surfaces in the space.

Radiation between the ceiling panel and other surfaces in the enclosure is an important parameter for heat exchange. The floor and the four walls of the enclosure are considered to participate in the radiant heat transfer as black bodies. However, the air within the enclosure is assumed to be a nonparticipating medium.

Modeling of heat source plumes in the presence of temperature gradients will follow the expressions given by Mundt (1996). The plume flow rate through the layer boundary depends on the plume source strength and vertical temperature gradient. In a region of uniform temperature, the volume flow rate [Q.sub.p] at an elevation Z above the ground level of the plume source with convective heat flux [PHI] is given by Mundt (1996) as

[Q.sub.p] = 0.005[[PHI].sup.1/3] (Z - [Z.sub.p])[.sup.5/3], (15)

where [Z.sub.p] is the source height from the ground level. Plumes are influenced by temperature stratification. The flow rate of a plume generated by a heat source of flux [PHI] in a space with temperature gradient [dT.sub.a]/dZ can be calculated by the following steps (Mundt 1996):

[Q.sub.p] = 0.00238[[PHI].sup.3/4]([dT.sub.a]/dZ)[.sup.-5/8] [B.sub.1] (16a)

[B.sub.1] = 0.004 + 0.039[A.sub.1] + 0.38[A.sub.1.sup.2] - 0.062[A.sub.1.sup.3] (16b)

[A.sub.1] = 2.86(Z - [Z.sub.p])([dT.sub.a]/dZ)[.sup.3/8] [[PHI].sup.-1/4] (16c)

In the presence of a temperature gradient [dT.sub.a]/dZ, the equilibrium height ([Z.sub.eq]), where the temperature difference between the plume and the ambient air disappears, is given by Mundt (1996) as

[Z.sub.eq] = 0.74[PHI][.sup.1/4]([dT.sub.a]/dz)[.sup.-3/8] (17)

The plume equilibrium heights for a point source of 100 W (341 Btu/h) are 2.34 m (7.7 ft) and 1.55 m (5 ft) for air temperature gradients of 1 K/m (0.55[degrees]F/ft) and 3 K/m (1.65[degrees]F/ft), respectively.

To solve the energy balance equations of the enclosure air layers, the values of the airflow rates passing from one layer to the other due to the generated plumes from the heat sources and the walls are calculated from the plume equation for the sources and using the wall plume model tested in this work, respectively.

The developed plume-multilayer mathematical model can now be translated into an efficient numerical model that solves for the vertical temperature distribution in a room conditioned by a combined DV/CC system in the presence of heat sources and nonuniform room wall conditions. The energy and mass balances represent a set of coupled nonlinear equations that are iteratively solved to obtain the layers and wall segment temperatures.

The numerical solution starts by assuming a temperature distribution for air layers and wall segments for a given temperature difference profile. The program proceeds by calculating the heat source plumes and the wall plumes. The set of algebraic equations resulting from the energy balance on the air layers and the walls were solved and new values of temperatures are calculated. The temperatures are updated and this procedure is repeated until convergence is reached when the change in values of temperature is less than 0.001[degrees]C (0.0018[degrees]F). A curve fit is then made for the wall-to-air temperature difference. If this curve fit matches the one assumed initially, then convergence is reached; otherwise, the new curve fit is used and the procedure is repeated until values of wall and air temperatures converge. The accuracy of the plume-multilayer model was tested by comparison to three-dimensional CFD simulations of an enclosure with heat sources using Airpak (2002) software interfaced with Fluent (2004) software.

PLUME-MULTILAYER MODEL TESTING

The validation of the complete modified plume-multilayer model is made by comparing the model outputs with values obtained from three-dimensional simulations of an enclosure conditioned by a DV/CC system using Airpak (2002) software, which was specifically developed for indoor airflow simulations. Chen and Xu (1998) tested the accuracy and reliability of a zero-equation turbulence model for DV-cooled rooms by experimentation and comparison to known CFD models. The agreement between the computed velocity and temperature vertical distributions was reported by Chen and Xu (1998) to be reasonably good. For that reason, the zero-equation turbulence model was adopted in the Airpak three-dimensional simulations of the test case comparison with the plume-multilayer model. The comparison between both models were based on the vertical temperature gradient, wall temperatures, cooling panel load, and stratification height. The simulations were done on a single geometry for three different heat loads (low, medium, and high).

The enclosure shown in Figure 10 has dimensions of 5 x 4 x 3 m (16.4 x 13.1 x 9.8 ft). Air is introduced into the room through a 0.6 x 1 m (1.96 x 3.28 ft) grille positioned at an elevation of 0.1 m (0.328 ft) above the ground and at a distance of 2 m (6.56 ft) from the edge of the south wall. The supply mass flow rate temperature is 22[degrees]C (71.5[degrees]F). The return air leaves the room through the exhaust grille, which is placed at the north wall, 2 m (6.56 ft) from its edge and 0.1 m (0.328 ft) below the ceiling level. The internal heat sources have a rectangular shape of 0.1 x 0.1 m (0.328 x 0.328 ft) in area and are placed at a height of 0.2 m (0.66 ft) above the floor. The heat flux of the internal heat sources is set to be convective only. The external heat load is conducted through the walls from outside conditions at fixed Tout and zero solar gain. The wall conductance is set to 1.53 W/K x [m.sup.2] (0.27 Btu/h x [ft.sup.2] x [degrees]F). The outside heat transfer coefficient is 21.1 W/K x [m.sup.2] (0.372 Btu/h x [ft.sup.2] x [degrees]F) for the walls and the floor. The supply mass flow rate, ceiling temperature, heat sources, and outside air temperature varied between low, medium, and high load cases. The different parameters used in the various cases are provided in Table 1.

The three cases were chosen based on the design chart of Behne (1999). The supply volume flow rates per unit area for the low, medium, and high load cases were 22.5 [m.sup.3]/h x [m.sup.2], 30 [m.sup.3]/h x [m.sup.2], and 37.5 [m.sup.3]/h x [m.sup.2] (1.23 [ft.sup.3]/min x [ft.sup.2], 1.64 [ft.sup.3]/min x [ft.sup.2], and 2.05 [ft.sup.3]/min x [ft.sup.2]), respectively. The low-load case is a lower limit case for the use of the combined DV/CC system. For the medium-load case, Behne (1999) recommends the use of DV or combining DV with a chilled ceiling. For the high load case, only the combined DV/CC system is advised by Behne (1999). The problem is set as steady with the radiation model turned on. The view factors are calculated and stored by the Airpak (2002) software. The number of nodes selected for the simulations is 11,950 (21 x 22 x 22) after checking solution grid independence at 76,015 (40 x 40 x 40) nodes and 231,141 nodes (60 x 60 x 60). Average heat flux load on the ceiling panel for the medium-load case was 51.5 W/[m.sup.2] (16.3 Btu/h) on the coarse grid compared to 50.47 W/[m.sup.2] (15.32 Btu/h) for the finest grid result. The convergence criteria were 0.001 m/s for the flow and [10.sup.-6] W for energy. The under-relaxation factors were the default values: 0.7 for pressure, 0.3 for momentum, 1.0 for temperature and viscosity, and 0.1 for body forces. Convergence was reached after 200-250 iterations.

[FIGURE 10 OMITTED]

The plume-multilayer model was applied on the three test cases dividing the room into four air layers of equal height and using the offset height [Z.sub.o] equal to the height of the first layer. The developed procedure for determining n and N is applied to each case to iterate for a converged solution starting from an initial linear temperature distribution between the floor and the ceiling temperatures. The wall plume flow rates, the average temperature of each air layer, and the associated average wall segment temperature were then obtained from the converged solution. The values of the parameters (N, n) derived from the model for the low, medium, and high load test cases were (0.48, -0.16), (0.82, -0.2), and (1.61, -0.14), respectively. The vertical temperature gradient in the room is one of the main parameters that affects thermal comfort inside the room, while the stratification height affects indoor air quality.

The Airpak (2002) software calculates the temperature at each node in the thermal space. For the vertical temperature gradient comparison to be valid, an average temperature for each layer must be calculated from the Airpak results. The average air layer temperature in Airpak is calculated by averaging the temperature of all the grid points in the layer volume. Figure 11 shows the average air temperature of each air layer of the four layers calculated from the plume-multilayer model and Airpak for the (a) low, (b) medium, and (c) high load cases. The maximum differences in the temperature values are 0.13[degrees]C (0.23[degrees]F), 0.37[degrees]C (0.66[degrees]F), and 0.3[degrees]C (0.54[degrees]F) for the low, medium, and high load cases, respectively.

The second comparison of the model values is based on the wall temperature distribution. The calculation of the temperatures of the wall parts at the different layers is done by averaging the wall temperature values over the wall width at each level using Airpak (2002) simulation results. The values of wall temperatures for the east and west walls are the same due to the symmetry of the room. Figure 12 shows the variation of the average vertical wall temperature associated with air layers for all the walls as predicted by the plume-multilayer model and Airpak for the low, medium, and high load cases for the east/west and north/south temperatures. The maximum error in the east/west wall temperature is about 0.41[degrees]C (0.72[degrees]F) for the medium load case. For the north and south walls, the maximum wall temperature differs by 0.5[degrees]C (0.9[degrees]F) and 0.8[degrees]C (1.44[degrees]F), respectively, due to the existence of the supply grille on the south wall and the return grille on the north wall. For the medium load case, a maximum temperature difference of about 0.83[degrees]C (1.5[degrees]F) exists between the model and Airpak results at the lower levels and a temperature difference of about 0.69[degrees]C (1.24F) exists at the higher levels. In general, wall plumes are driven by the temperature difference between the wall and the room air. This difference is inherently not uniform in DV systems, and that effect is compounded when the walls are not well insulated, which could result in a larger temperature difference between the wall and the room air.

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

The panel cooling load is calculated from Airpak (2002) software by multiplying the resulting panel heat flux by the ceiling area. The panel cooling loads computed by the plume-multilayer model for low, medium, and high loads were 574 W (1,958 Btu/h), 905.4 W (3,089 Btu/h), and 1,430 W (4,878 Btu/h), respectively, compared to three-dimensional simulation results of 530 W (1,808 Btu/h), 1,022 W (3,486 Btu/h), and 1506.4 W (5,139 Btu/h), respectively. The maximum error is less than 11% between model predictions and three-dimensional simulations. The intention for accurate estimation of the flow rate is to be able to predict the stratification height, which is an important parameter for thermal comfort and air quality considerations. The model prediction for the stratification height inside the room is based only on mass balance equations. The stratification height was predicted by the model to be at a level of 1.372 m (4.5 ft), 1.16 m (3.8 ft), and 1.15 m (3.77 ft) for the low, medium, and high load cases, respectively.

The three-dimensional flow simulation inside the enclosure represents a case of mixed convection, while the present model assumes no interaction between sources and the wall. The stratification height is distinguished in the three-dimensional flow enclosure results by the highest level at which no noticeable downward flow exists. The upward and downward mass flow rates were calculated in the space at different elevations by integrating the vertical velocity component distribution over the horizontal surface area. The variation of the upward and downward mass flow rates with elevation for the low, medium, and high load cases is shown in Figures 13a, 13b, and 13c, respectively. The up and down mass flow rates remain constant below the stratification height region and then slightly change slope in the upper region. The stratification heights predicted by the model are shown directly in Figure 12, which indicates that their values fall closely in the region when a slight change in upward mass flow rate takes place. The model predictions are within 0.05 m (0.16 ft) of the three-dimensional model.

[FIGURE 13 OMITTED]

The importance of including wall plumes in DV with chilled ceiling models is seen by examining the stratification height and the vertical temperature gradient when wall plumes are neglected or isothermal plumes are used. Figure 14 shows the variation of (a) the stratification height and (b) the vertical temperature gradient, with supply flow rate at a fixed room perimeter-to-area ratio while varying the ceiling temperature to meet the steady-state load of the medium case. For a cooling load of 67.3W/[m.sup.2] (21.31 Btu/h x [ft.sup.2]), a zero wall-plume value overestimates the stratification height by 35% within the space, while an isothermal plume overestimates the stratification height by 8% to 10%. The vertical temperature gradient is decreased by 0.07 K/m (0.04[degrees]F/ft) when wall plumes are in the model compared to zero wall plume.

The model accurately predicted the stratification height when compared to the three-dimensional CFD model. The added complexity of calculating the wall plumes over isothermal correlation using the power law similarity model improved predictions of wall flow rates and stratification height to become comparable to three-dimensional full CFD simulations at a much lower computational cost.

CONCLUSIONS

A new, improved simple model was developed to represent room heat transfer with a combined DV/CC system. The model is based on energy conservation equations in air layers while deriving airflow rates exchanged between adjacent layers from plume equations of Mundt (1996) and from a new wall plume model developed in this study. The heat source induced flow rate relied on plume equations of point and area heat sources of Mundt (1996). The wall plume flow rates are based on a similarity solution of two-dimensional boundary layer flow with power law temperature difference distribution between the heated wall and the environment air. The determination of the power law parameters of the temperature difference distribution was subject to rigorous testing on the selected exponent and the temperature difference coefficient that predicts the wall plumes most closely to the real physical situation occurring in the space.

The plume-multilayer model using the developed wall plume model was validated in comparison with detailed three-dimensional simulation of the enclosure, and results of the average air and wall temperature variations with elevation, the stratification height, and the panel cooling load agreed well with the plume-multilayer model predictions. The model has very good accuracy at a low computational cost and with sufficient flexibility in the geometry, source strengths, and supply flow rates of real systems. The flexibility stems from the physical consideration of the mechanism of convection that takes place in the enclosures conditioned by DV and a chilled ceiling. Future work will extend the model to include turbulent plume regimes and to use the model in transient application when changes in the load or external environment take place.

ACKNOWLEDGMENTS

The authors acknowledge the support of the University Research Board of American University of Beirut research grant DDF-113040-688324, the Lebanese National Council for Scientific Research grant LCR-113040-002208, and ASHRAE grant 113040-025043.

[FIGURE 14 OMITTED]

NOMENCLATURE

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

CC = chilled ceiling

CFD = computational fluid dynamics

d = temperature difference between the wall and the air, [degrees]C ([degrees]F)

[d.sub.max] = maximum temperature difference between the wall and the air, [degrees]C ([degrees]F)

DV = displacement ventilation

f = dimensionless stream function

[F.sub.i[right arrow]j]= view factor from surface i to surface j

g = gravitational acceleration

[Gr.sub.Z] = Grashof number, [Gr.sub.Z] = g[Z.sup.3]/[v.sup.2] [beta]d(Z)

H = enthalpy, J/kg (Btu/1bm)

h = heat transfer coefficient, W/[m.sup.2] x K (Btu/h x [ft.sup.2] x [degrees]F)

J = number of surfaces in the space

L = room length, m

M = mass flow rate, kg/s (lbm/h)

n = power of the curve fit

N = power law coefficient

Pr = Prandtl number

q = heat flux, W/[m.sup.2] (Btu/h x [ft.sup.2])

[Q.sub.p] = plume volumetric flow rate, [m.sup.3]/s ([ft.sup.3]/min)

[Q.sub.w] = wall plume volumetric flow rate per unit width, [m.sup.3]/m x s ([ft.sup.3]/ft x min)

r = number of heat sources in the first air layer

Re = Reynolds number

[R.sup.2] = data fit correlation parameter

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

t = time, s

U = wall conductance, W/[m.sup.2] x K (Btu/h x [ft.sup.2] x [degrees]F)

u = horizontal component of the velocity, m/s (ft/min)

v = vertical component of the velocity, m/s (ft/min)

W = width of the room, m (ft)

[Y.sub.(j)] = width of wall j, m (ft)

Z = elevation above ground level, m (ft)

[Z.sub.o] = offset elevation, m (ft)

[Z.sub.p] = source height from floor level, m (ft)

[Z.sub.eq] = plume equilibrium height (ft)

Greek Symbols

[alpha] = thermal diffusivity, [m.sup.2]/s ([ft.sup.2/s]

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

[DELTA]T = difference between average or isothermal wall temperature and far fluid temperature [degrees]C ([degrees]F)

[PHI] = convective heat flux, W (Btu/h)

[phi] = dimensionless temperature, (T - [T.sub.a])/d(Z)

[eta] = transformation parameter, [eta] = y/z4[square root of ([Gr.sub.z]/4)]

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

[upsilon] = kinematic viscosity, [m.sup.2]/s ([ft.sup.2]/s)

Subscripts

a = air

c = ceiling

cir = circulation

cv = convection

f = floor

l = light

o = outer

p = plume

rd = radiation

s = supply

solar = solar radiation

si = inner surface wall temperature

so = outer surface wall temperature

w = wall

REFERENCES

Airpak. 2002. Airpak 2.1, Computational Fluid Dynamics Software for Air Conditioning Applications, Fluent, Inc., Lebanon, NH. www.fluent.com.

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

Behne, M. 1999. Indoor air quality in rooms with cooled ceilings: Mixing ventilation or rather displacement ventilation? Energy and Buildings 30(2):155-66.

Carrilho da Graca, G.C.C. 2003. Simplified models for heat transfer in rooms. PhD thesis, University of California, San Diego.

Chapra, S.C., and R.P. Canale. 2002. Numerical Methods for Engineers, with Programming and Software Applications, 4th ed. Boston: W.C. Brown/McGraw-Hill.

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

Chen, Q., and W. Xu. 1998. A zero-equation turbulence model for indoor airflow simulation. Energy and Buildings 28(2):137-44.

Etheridge, D., and M. Sandberg. 1996. Building Ventilation: Theory and Measurement. Chichester: John Wiley and Sons.

Fluent. 2004. Computational Fluid Dynamics Software, Fluent 6.2.16, Fluent, Inc., Lebanon, NH. www.fluent.com.

Gebhart, B., Y. Jaluria, R. Mahajan, and B. Sammakia. 1988. Buoyancy Induced Flows and Transport. New York: Hemisphere Publishing Corporation.

Goodfellow, H.D., and T. Esko, eds. 2002. Industrial Ventilation Design Guidebook. San Diego, CA: Academic Press.

Jaluria, Y. 1980. Natural Convection, Heat and Mass Transfer. New York: Pergamon Press.

Jiang, Z., Q. Chen, and A. Moser. 1992. Indoor airflow with cooling panel and radiative/convective heat source. ASHRAE Transactions 98(1):33-42.

Linden, P.F., G.F. Lane-Serff, and D.A. Smeed. 1990. Emptying filing boxes: The fluid mechanics of natural ventilation. Journal of Fluid Mechanics 212:300-35.

Min, T.C., L.F. Schutrum, G.V. Parmelee, and J.D. Vouris. 1956. Natural convection and radiation in a panel heated room. ASHRAE Transactions 62:337-58.

Mundt, E. 1996. The Performance of Displacement Ventilation Systems: Experimental and Theoretical Studies. Stockholm, Sweden: Royal Institute of Technology.

Novoselac, A., and J. Srebric. 2002. A critical review of the performance and design of combined cooled ceiling and displacement ventilation systems. Energy and Buildings 34:497-509.

Rees, S.J. 1998. Modeling of displacement ventilation and chilled ceiling systems using nodal models. PhD thesis, Loughborough University, London, UK.

Rees, S.J., J.J. McGuirk, and P. Haves. 2001. Numerical investigation of transient buoyant flow in a room with a displacement ventilation and chilled ceiling system. International Journal of Heat and Mass Transfer 44(16):3067-80.

Skaaret, E. 1998. A semi-empirical flow model for low-velocity air supply in displacement ventilation. Proceedings of Roomvent 98, Stockholm, Sweden 1:85-92.

Yuan, X., Q. Chen, and L. Glicksman. 1998. A critical review of displacement ventilation. ASHRAE Transactions 104(1):78-90.

Mohamad Ayoub

Nesreen Ghaddar, PhD

Member ASHRAE

Kamel Ghali, PhD

Received February 8, 2006; accepted April 28, 2006

Mohamad Ayoub is a graduate research assistant and Nesreen Ghaddar is a professor and Qatar Chair for Energy Studies in the Department of Mechanical Engineering, American University of Beirut, Beirut, Lebanon. Kamel Ghali is an associate professor in the Department of Mechanical Engineering, Beirut Arab University, Beirut, Lebanon.

The wall plume model is integrated with a multilayer space thermal model to predict the stratification height in the space, the vertical distribution of wall and air temperatures as a function of the chilled ceiling temperature and space air supply conditions. The plume-multilayer model results were compared with results of three-dimensional CFD simulations using commercially available software. Three test cases were considered for the simulations at cooling loads of 40 W/[m.sup.2] (12.7 Btu/h x [ft.sup.2]), 67 W/[m.sup.2] (21.2 Btu/h x [ft.sup.2]), and 100 W/[m.sup.2] (31.7 Btu/h x [ft.sup.2]) and supply airflow rates per unit area of 22.5 (1.2 [ft.sup.3]/min x [ft.sup.2]), 30 (1.6 [ft.sup.3]/min x [ft.sup.2]), and 37.5 [m.sup.3]/h x [m.sup.2] (2.1 [ft.sup.3]/min x [ft.sup.2]), respectively. The vertical wall and average air temperatures for each layer agreed well with the results of the plume-multilayer model, showing maximum errors in values of average air temperature of 0.13[degrees]C (0.23[degrees]F), 0.37[degrees]C (0.66[degrees]F), and 0.3[degrees]C (0.54[degrees]F) for the low, medium, and high load cases, respectively. The simplified model accurately predicted the stratification height at a maximum error of [+ or -]0.05 m (0.16 ft) in the three test cases. The stratification height is overestimated by 35% if wall plumes are neglected in the plume-multilayer model.

INTRODUCTION

The energy efficiency of HVAC systems in buildings has a major impact on building energy consumption. One of the HVAC systems that has potential for energy savings is the combined displacement ventilation and chilled ceiling (DV/CC) system (Novoselac and Srebric 2002). In displacement ventilation (DV), air is supplied into the room at floor level, at a low velocity, and at a temperature slightly lower than the desired room temperature. Thus, the cooler air will displace the warmer room air that rises by buoyancy, creating two zones, as reported by Jiang et al. (1992) and Yuan et al. (1998). The lower occupied zone contains fresh, cool air, and the upper zone contains the warm, contaminated recirculation air, which is exhausted near the ceiling. Maintaining the interface between these two layers above the occupied zone is one of the many unique challenges of DV design (Yuan et al. 1998).

Vertical temperature variation and thermal draft are more critical in the design of DV systems. According to ASHRAE (2005), thermal comfort considerations limit the supply airflow rate (recommended velocity ~0.15 m/s [29.53 ft/min], maximum velocity 0.25 m/s [49.21 ft/min]) to avoid cold drafts, and the supply air temperature (minimum temperature limit is 18[degrees]C [64.4[degrees]F]). It also recommends that the maximum vertical temperature gradient in the room should be lower than 2.5[degrees]C/m (1.4[degrees]F/ft) if the percent dissatisfied is to be kept less than 5%. Yuan et al. (1998) indicated that the maximum removed load from the occupied zone by displacement systems is approximately 40.0 W/[m.sup.2] (12.7 Btu/h x [ft.sup.2]). The Behne (1999) design diagram sets the cooling load limit of DV when combined with a chilled ceiling to 100.0 W/[m.sup.2] (31.7 Btu/h x [ft.sup.2]) of floor area.

One of the primary challenges in modeling enclosures cooled by combined DV/CC systems is their inherent complexity and detail-oriented nature, which makes simulations computationally expensive and, hence, hinders their integration in building energy simulation programs for optimization and control purposes. The modeling of DV has followed several approaches ranging from full simulations using computational fluid dynamics (CFD) tools (Chen 1995; Rees et al. 2001) to experimentally based semi-empirical models (Skaaret 1998) to nodal models (Rees 1998) to plume two-zone models (Mundt 1996; Carrilho da Graca 2003).

CFD modeling can give reasonable accuracy when simulations are made by expert users and when the suitable boundary conditions are imposed during the simulation. Chen (1995) compared the numerical results of five different models and recommended that the renormalized group (RNG) k-[epsilon] model be used for predicting indoor airflows. Good agreement was also found when comparing CFD simulations and test chamber measurements, as reported by Rees et al. (2001). CFD modeling is the most accurate, but it is prohibitively expensive when studying multi-room ventilation geometries or transient performance over a period of days or months.

A simpler approach for DV modeling is to use experimentally based semi-empirical correlations developed for particular geometries based on simple and clear approximations (Skaaret 1998). These correlations are suitable for use in the design stage of DV systems, but they give reliable results only when the geometries under consideration are similar to those for which the correlations were developed. A DV nodal model reported by Rees (1998), which represents the space by a network of nodes upon which heat and mass conservation equations are applied, is more general than the semi-empirical approach and is computationally less expensive than the CFD approach. However, the interaction between the different nodes depends on precalculated airflow rates that should be obtained experimentally or by CFD simulations for the particular geometries under consideration.

The plume two-zone model is derived from fundamental equations of thermal transport from plumes induced by heat sources in the space. This approach eliminates the need to perform CFD simulations to calculate plume flow rates. Plume two-zone models take a physical approach to the problem, where the space is divided into two horizontal air layers and the mass flow in and out of each layer due to heat sources is calculated from the plume equations reported by Mundt (1996). Previous work on plume-multilayer wall models has either neglected wall plumes or used isothermal wall flow models, as in the work of Carrilho da Graca (2003), Mundt (1996), and Linden et al. (1990). The applicability of the plume two-zone model was studied by Carrilho da Graca (2003) in cases where plume convection was stronger than wall-driven natural convection. In his work, Carrilho da Graca (2003) developed a model that predicted three temperatures: floor-level temperature, occupied zone temperature, and a mixed-layer outflow temperature. The plumes were assumed to be of equal strength and located at the room floor. The wall plume convection was assumed negligible.

Wall plumes in hot weather conditions may not be neglected due to the high solar gain of the wall, which creates larger temperature differences between room air and wall surface. In many places with hot climates, building common practice uses concrete-based walls of high U-factor (U = 2.4 to 3.4 W/[m.sup.2] x K [0.4 to 0.6 Btu/h x [ft.sup.2] x [degrees]F]) that contributes to higher internal wall surface temperatures. The use of isothermal wall correlations in plume models in spaces where DV is combined with chilled ceilings is also incorrect. With DV, the nonuniformity in wall temperature is significant by the nature of the cooling mechanism. There is an apparent need for a modified model that deals with nonuniform wall conditions and heat sources at a low computational cost and without having to perform three-dimensional simulations to get the associated transport coefficients. A computationally efficient model that captures the physics of the problem would easily be incorporated with the HVAC system optimization and control program for reducing energy consumption under steady and dynamic conditions and for determining the building system efficiency. The use of a simple and accurate model is essential in finding means to improve the system response and energy consumption in response to changes to load within the enclosure.

For this study, a wall plume model was developed to accurately predict flow rates adjacent to a wall. The model solved the wall plume equations for two-dimensional boundary layer self-similar laminar flow while using a power law distribution for the temperature difference between the wall and the ambient conditions. The wall plume model was tested for accuracy when used in infinite and closed space domains. The wall-plume model was integrated with source plumes and a multilayer space model. The plume-multilayer model was validated by three-dimensional simulations using Airpak (2002) commercial software, and comparison was made with isothermal wall space conditions. The contribution of the work is the ability of the developed improved model to predict at low computational cost the DV/CC thermal transport, the stratification height, and the vertical air and wall temperature gradients.

MODEL MATHEMATICAL FORMULATION

The main feature of the plume-multilayer model is its ability to identify the fundamental physics of the problem in a simplified, computationally efficient approach. The modeling of the plumes generated by convective fluxes is a major part of the proposed model. The mass flow rates and the average temperature of the plumes are evaluated using expressions that take into account the strength of the convective flux and the position of the source as reported by Mundt (1996). Since the chilled ceiling temperature is relatively lower than the temperature of the other surfaces, the radiant heat exchange cannot be neglected. Unlike other models that assume conditions of adiabatic walls or isothermal walls, this model takes into account the nonuniformity of the wall temperature and the dependence of this temperature distribution on the internal and external conditions. The temperatures of the walls at different heights will be determined by applying the energy conservation equations to the wall surfaces. The model can predict the circulation (stratification) height. The convective and the radiative heat loads on the ceiling panel can also be predicted.

However, the model has limitations based on some of the assumptions. The heat sources are sufficiently distant to prevent plume interaction. Point heat sources are located at their actual positions, and heat sources with dimensions are located at equivalent virtual point source positions (Goodfellow and Esko 2002). It is also assumed that the equilibrium height of each plume is not less than three-fourths of the room height, where the equilibrium height is the elevation at which the density gradients disappear and the plume spreads horizontally, as reported by Goodfellow and Esko (2002).

Wall Plume Model

The stratification height and the temperatures of air and room walls are affected by the wall plumes. For an isothermal wall, the volume flow rate of the wall plume per unit width is given by Jaluria (1980) and Etheridge and Sandberg (1996) as

[Q.sub.w] =0.00287[DELTA][T.sup.1/4][Z.sup.3/4], (1)

where Z is the distance from the vertical wall leading edge and [DELTA]T is the temperature difference between the isothermal wall and the fluid far from the wall. For DV with a chilled ceiling system, the temperature differences between the top and the bottom of the wall are significant and, thus, the assumption of an isothermal wall may be overestimating or underestimating the wall plume flow rates. A more rigorous approach using the two-dimensional laminar boundary layer similarity analysis of Gebhart et al. (1988) in a stable, stratified environment was followed to determine the wall plumes. The two-dimensional transport boundary region governing equations of thermally induced motion are:

Continuity equation [[partial derivative]u/[partial derivative]x] + [[partial derivative]v/[partial derivative]Z] = 0 (2a)

Momentum equation u[[partial derivative]u/[partial derivative]x] + v[[partial derivative]v/[partial derivative]Z] = [v][[[partial derivative].sup.2]u/[partial derivative][Z.sup.2]] + g[beta](T - [T.sub.a])(2b)

Energy equation u[[partial derivative]T/[partial derivative]x] + v[[partial derivative]T/[partial derivative]Z] = a[[[partial derivative].sup.2]T/[partial derivative][Z.sup.2]] (2c)

subject to the following boundary conditions:

u(x,0) = 0

u(x,[infinity]) = 0

v(x,0) = 0

T(x,0) = [T.sub.w]

T(x,[infinity]) = [T.sub.a]

where

[T.sub.w] - [T.sub.a] = d(Z) (3)

where u is the horizontal component of the velocity, v is the vertical component, [upsilon] is the kinametic viscosity, [alpha] is the thermal diffusivity, g is the gravitational constant, and [beta] is the thermal expansion coefficient. The air and wall temperature difference decreases with height in spaces conditioned by a combined DV/CC system. A power law distribution represents well the temperature difference between the wall and the air d(Z) and is given by

d(Z) = [T.sub.w] - [T.sub.a] = N[Z.sup.n] (4)

where N is a constant. When n = 0, then an isothermal wall case is present if there is no variation in the air temperatures. The corresponding transformed boundary layer equations of momentum and energy derived by Gebhart et al. (1988) in terms of the transformation parameter [eta] are given in terms of a dimensionless stream function f([eta]) and dimensionless temperature [phi]([eta]) = (T - [T.sub.a])/d(Z) as follows:

f'" + (n + 3)ff" - (2n + 2)[f'.sup.2] + [phi] = 0 (5a)

[phi]" +Pr[(n + 3)f[phi]' - 4nf'[phi]] = 0 (5b)

subject to the following boundary conditions:

[[partial derivative]f/[partial derivative][eta]](0) = 0, [[partial derivative]f/[partial derivative][eta]]([infinity]) = 0, f(0) = 0, [phi]([infinity]) = 0, [phi](0) = 1 (5c)

where [eta] is defined as [eta](Z, y) = b(Z)y, and Pr is the Prandtl number. The function b(Z) is defined by

b(Z) = [1/Z]4[square root of ([Gr.sub.Z]/4)], (5d)

where [Gr.sub.Z] is the Grashof number ([Gr.sub.Z] = g[Z.sup.3]/[v.sup.2] [beta]d(Z)). The volume flow rate per unit width at any distance Z is calculated as

[Q.sub.w](Z) = 4v([Gr.sub.Z]/4)[.sup.1/4] [f.sub.[infinity]]. (6)

The volume flow rate per unit width of the wall plume [Q.sub.w] is directly related to the "shape" of the temperature distribution along the vertical wall rather than a computed average wall temperature. Gebhart et al. (1988) reported a range of the exponent n for which the solution is physical to be between -0.6 and unity (-3/5 < n < 1). The parameter [f.sub.[infinity]] depends only on the exponent n and is acquired from the numerical solution of ordinary differential Equations 5a and 5b using the shooting method (Charpa and Canale 2002) with first-order finite-difference discretization of f', f", f'", [phi]', and [phi]". The numerical solution is obtained for values of the exponent n between -0.6 and 1.0 in steps of 0.05 and increments of the independent variable [eta] at [DELTA][eta] = 0.0001. Figure 1 shows a plot of the value of the calculated parameter [f.sub.[infinity]] as a function of n in the desired range for Pr = 0.7. The solution for the special cases of the isothermal (n = 0) and constant heat flux (n = 0.2) wall presented by Jaluria (1980) are reproduced by the current general solution and are also shown in Figure 1. The [f.sub.[infinity]] - n relation is expressed as a third-order polynomial equation, with a correlation parameter value of 0.9998, as follows:

[FIGURE 1 OMITTED]

[f.sub.[infinity]] = 0.1145[n.sup.3] + 0.2271[n.sup.2] - 0.2978n + 0.6024 (7)

Since the real physical temperature difference distribution decreases with height in rooms conditioned with a DV/CC system, the maximum value of n should be less than zero. The range of exponent n for wall plumes is taken as -0.6 < n < 0.

The major challenge in the wall plume model is to develop a methodology to estimate the values of n and N that give accurate predictions of wall flow rates and yet provide a good power law curve fit of the calculated temperature difference between the wall and the air in the space.

Determination of the Power Law Parameters

The equations through which the mass flow rates are calculated depend on the power n and the coefficient N of the wall-to-air temperature difference distribution. For negative powers (n < 0), high temperature differences occur at low elevations and the difference is infinite at Z = 0. The similarity solution will result, then, in a relatively large mass flow rate at a low level, overestimating the real wall flow rates in the space. To overcome this problem, a low-height offset ([Z.sub.o]) of the wall is set to a uniform temperature equal to the local temperature at [Z.sub.o]. In practical cases, it can be assumed that the wall temperature is uniform at low levels. The flow rate of the plume at [Z.sub.o] resulting from the power law solution is matched to the flow rate produced at low height by the isothermal solution estimated from Equation 1. A valid curve fit for the wall-to-air temperature difference is then derived to give the highest correlation factor ([R.sup.2]) between the data values produced from the plume-multilayer model. The curve fit must also give the same flow rate given by the isothermal equation at [Z.sub.o] to offset the effect of the high flow rate in the region near Z = 0.

Thus, the calculation of the coefficient N and the power n used in determining the plume of a non-isothermal wall is subjected to two constraints. The first constraint is that N and n must give a flow rate close to the value given by the isothermal equation for low elevation at a selected height [Z.sub.o]. The second constraint is that the error between the corresponding temperature difference curve fit using the least squares method and the actual temperature difference is minimal or, alternatively, that the correlation parameter [R.sup.2] is maximal.

The method starts by determining the coefficient N using the least squares method for each power in the range of n = -0.6 and n = 0 in steps of [DELTA]n = 0.01 to generate the curve fits at these values of N and n. Then the corresponding wall plume flow rates are calculated from Equations 5 and 6 at the desired offset elevation [Z.sub.o]. These flow rates are then compared with the flow rate calculated at [Z.sub.o] from the isothermal equation (1), assuming the wall, up to the height [Z.sub.o], is at uniform temperature equal to the value of the wall temperature at [Z.sub.o]. The values of n and N that give the nearest flow rate to that of the isothermal equation with the highest [R.sup.2] value will be selected. The value of n could be refined by using smaller steps of [DELTA]n. The effect of the selection of the offset height [Z.sub.o] on the model results were tested at values of [Z.sub.o] = 0.15 m (0.49 ft), 0.3 m (0.98 ft), and 0.6 m (1.97 ft).

Testing the Wall Plume Model

The nonuniform thermal conditions of the wall used in the wall-plume model testing are a power law temperature difference distribution for Z > [Z.sub.o] and an isothermal wall at [T.sub.w]([Z.sub.o]) for Z [less then or equal to] [Z.sub.o]. The validation of the model is done by testing the accuracy of the followed procedure to determine the power law parameters that predict accurate flow rates.

The wall plume model methodology of selection of power law parameters and the accuracy of predicted flow rates are tested in the following ways:

* By making comparisons of the model wall flow rate predictions to results of Fluent (2004) CFD simulations of two-dimensional boundary layer flow adjacent to a nonuniformly heated wall in an infinite domain. The similarity solution from which the wall plume model was formulated is based on the assumption of an infinite domain. As a first step, the values of the flow rates calculated from the wall plume model are compared with flow rates calculated from infinite domain CFD tests. The aim of the infinite domain test is to check the suggested methodology for choosing the power n and the coefficient N that will be used in the flow rate calculation. The infinite domain test permits applying a power law temperature difference distribution while keeping the infinite domain air temperature constant. This condition is not as simple to apply in a closed domain, when the room vertical temperature variation is not known a priori. The second aim is to compare the flow rates from the CFD results with flow rates given by the isothermal equation.

* By making comparisons of model wall flow rate predictions to results of Fluent (2004) CFD simulations of two-dimensional flow in an enclosure with zero source plumes and in the presence of a chilled ceiling to isolate the effect of wall plumes due to nonuniformity of wall temperature. The enclosure test will also check the wall plume model applicability when both room air and room wall temperature distributions vary with height while maintaining a power law distribution of the temperature difference.

Infinite Domain CFD Test. Although many tests were performed to validate the model, the results are presented for only a few typical temperature difference distributions. The temperature difference distributions were based on the published experimental data of Rees (1998) and span the possible realistic range of temperature differences in rooms cooled by a combined DV/CC system. In infinite domains, the temperature difference between the wall and air is controlled by setting the air temperature at a constant and assigning the wall temperature distribution that satisfies the power law. The infinite domain test cases of natural convection of air at a uniform temperature of 23[degrees]C (73.15[degrees]C) adjacent to a vertically heated wall were simulated for a wall height of 4 m (13.12 ft) with seven wall temperature distributions, as shown in Figure 2. The first five temperature distributions had the maximum temperature difference at the lowest elevation less than 1.5[degrees]C (2.7[degrees]F). The cases FC6 and FC7 present limiting cases where the maximum temperature difference is equal to 0.72[degrees]C (1.3[degrees]F) and 2.8[degrees]C (5[degrees]F), respectively. The FC6 case is directly obtained from the experimental results of Rees (1998). The FC7 case was constrained by the condition that the flow at 3.0 m (9.8 ft) height remained laminar with the Rayleigh number below [10.sup.9]. The temperature difference distribution of case FC7 represented the maximum possible temperature difference at the lowest level that could maintain a laminar flow below a height of 3.0 m (9.8 ft).

The computational domain used in the Fluent (2004) simulations is extended in the horizontal direction to 2 m (6.56 ft). The heated wall temperatures were set such that for [Z.sub.o] [less then or equal to] Z [less then or equal to] 3.3 m (10.8 ft), the temperature difference between the wall and the air has a power law distribution. The heated wall temperature at Z [greater then or equal to] 3.3 m (10.8 ft) is set equal to that of the environment at 23[degrees]C (73.15[degrees]F) and for Z [less then or equal to] [Z.sub.o] the wall is isothermal at [T.sub.w]([Z.sub.o]) (e.g., case FC1: [T.sub.w]([Z.sub.o])=[T.sub.w] [0.3 m] = 24.5[degrees]C [75.8[degrees]F] for Z [less then or equal to] 0.3 m [0.98 ft]). The other boundaries of the domain were set at a constant temperature of 23[degrees]C (73.15[degrees]F), while the lower and upper boundaries were set as pressure inlet and outlet boundaries, respectively. The number of cells in the computational domain was 12,535. The 110 nodes in the vertical direction were evenly spaced, while the 127 nodes in the horizontal direction had variable distribution with a very fine grid close to the heated wall. The minimum and maximum values of the cell area were 1.135696 x [10.sup.-6] [m.sup.2] (1.22 x [10.sup.-5] [ft.sup.2]) and 3.965131 x [10.sup.-6] [m.sup.2] (4.268 x [10.sup.-5] [ft.sup.2]), while the minimum and maximum values of the cell face length were 4.5549 x [10.sup.-5] m (4.903 x [10.sup.-4] [ft.sup.2]) and 1.1172 x [10.sup.-1] m (1.2025 [ft.sup.2]). Convergence was reached after 1,159 iterations. Grid independence tests were also performed at a finer grid of 220 x 254 and showed no significant difference in velocity values. The velocity profiles are determined at different elevations. Figure 3 shows the variation of the vertical velocity component in the x direction at an elevation equal to 3.0 m (9.8 ft). The vertical velocity values are integrated over the x horizontal distance to obtain the wall flow rates.

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

Flow rates produced by Fluent (2004) software are compared with values obtained from the derived wall plume equations used in the multilayer model. Figure 4 shows the mass flow rates of the wall plumes calculated using the power law similarity solution for different fits as compared to the values calculated from two-dimensional CFD simulations for test case FC1 with offset height [Z.sub.o] = 0.3 m (0.98 ft). For the case FC1 distribution, the best power law parameters n and N resulting from applying the selection procedure of the model were found to be -0.23 and 0.5792, while the [R.sup.2] value was equal to 0.701. The mass flow rate calculated from the wall plume model at n = -0.23 agrees well with those calculated from the simulations. The maximum percentage error for n = -0.23 is 2.9%. The percentage error from the isothermal (average temperature, n = 0) case varied between 16% and 10%, while the percentage error for the best fit case of similarity wall plume model at n = -0.6 ([R.sup.2] = 1) was 45%. Similar results and accuracy of flow rates as those found for case FC1 were obtained for the other cases. Note that the wall plume flow rates below the offset height [Z.sub.o] (Z [less then or equal to] [Z.sub.o]) must be calculated based on the isothermal equation rather than the wall plume model.

[FIGURE 4 OMITTED]

Figure 5a shows the local Rayleigh number as a function of elevation for the FC6, FC7, and FC1 cases at [d.sub.max] = 0.72[degrees]C (1.3[degrees]F), 2.8[degrees]C (5[degrees]F), and 1.5[degrees]C (2.7[degrees]F), respectively. The isothermal region elevation is taken at 0.3 m (0.98 ft). The flow rates are calculated for the different fits of each case and compared to the values obtained from the isothermal equation. Results from the CFD simulation of each case are compared to the flow rates calculated from different fits using the wall plume model. The variations of the flow rate with elevation using the fit that matches the isothermal flow rate at 0.3 m (0.98 ft) elevation and using CFD predictions for cases FC6 and FC7 are shown in Figures 5b and 5c, respectively. The maximum error in flow rates is less than 10%. This indicates that the model is insensitive to the maximum temperature difference present in the distribution.

The effect of the elevation at which the comparison between the flow rates from the isothermal equation and the model is made (0.3 m in the above cases) was considered by selecting three different values of [Z.sub.o] for case FC8 at 0.15 m (0.49 ft), 0.3 m (0.98 ft), and 0.6 m (1.97 ft). The FC8 case is the same as the FC1 case but with the first 0.6 m of the heated wall set at a temperature of 24[degrees]C (75[degrees]F). The wall plume model was used to calculate the wall plume flow rate as a function of elevation. The optimal power n was -0.19, -0.21, and -0.25, and the corresponding optimal coefficients N were 0.48, 0.48, and 0.49 for the cases of 0.15 m (0.49 ft), 0.3 m (0.98 ft), and 0.6 m (1.97 ft) offset heights, respectively. Figure 6 shows the variation of wall plume model flow rates at optimal fits at 0.15 m (0.49 ft), 0.3 m (0.98 ft), and 0.6 m (1.97 ft) offset heights of case FC8. The maximum percentage difference in the predicted flow rate of the different fits is about 6.75% and occurs at a height of 0.15 m (0.49 ft). The percentage difference decreases as the elevation increases. This implies that the sensitivity of the model to the isothermal height introduces minimal error in the value of mass flow rates above [Z.sub.o].

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

The infinite domain test results show that the suggested methodology to calculate the flow rates using open-domain similarity solutions of adjusted power law parameters can accurately predict flow rates of a nonuniformly heated vertical wall with power law distribution with an isothermal initial starting length. The model reproduces the two-dimensional CFD flow rate predictions in an open domain with the same heated wall temperature distribution.

Two-Dimensional Enclosure CFD Test. The applicability of the wall plume model obtained from the infinite domain solution were tested in closed domains (enclosures) of typical room heights at zero internal heat sources. Two-dimensional simulations of thermally driven flow in a rectangular enclosure with nonuniformly heated walls in the presence of a chilled ceiling boundary were performed using CFD software (Fluent 2004). Unlike infinite domain, the CFD solution will generate room air temperature variations. These temperatures will then be used as input to the wall plume model to compare resulting flow rates with the two-dimensional enclosure flow rates produced by CFD and with flow rates based on isothermal flow equations. Note that when space dimensions are large, the wall plume flow rates obtained from CFD simulation will approach the values obtained from an infinite domain.

A symmetrical 4 m (13.1 ft) wide and 3 m (9.840 ft) high two-dimensional enclosure was simulated. Due to symmetry, only half of the physical domain was considered in the computational domain and the symmetry boundary was represented by an adiabatic wall. The upper boundary (ceiling) and the lower boundary (floor) of the rectangle are assigned constant surface temperature conditions of 18[degrees]C (64.4[degrees]F) and 23[degrees]C (75[degrees]F), respectively. Cases FE1, FE2, FE3, and FE4 are considered in an enclosure corresponding to various temperature difference distributions between the vertical wall and air with computed room air temperatures of 22.9[degrees]C (73.2[degrees]F), 22.5[degrees]C (72.5[degrees]F), 21.89[degrees]C (71.4[degrees]F), and 21.9[degrees]C (71.4[degrees]F), respectively, and corresponding maximum temperature differences of 4.9[degrees]C (8.8[degrees]F), 4.0[degrees]C (7.2[degrees]F), 2.6[degrees]C (4.7[degrees]F), and 2.3[degrees]C (4.1[degrees]F), respectively, as shown in Figure 7a. In Figure 7b, a variation of the vertical velocity component is plotted in the horizontal direction at an elevation of 2.7 m (8.86 ft) for case FE1. It is clear how the velocity profile in the boundary region is disturbed. When calculating the mass flow rate, the boundary layer region is considered to extend from the wall to the point where the velocity profile sharply changes slope.

[FIGURE 7 OMITTED]

Following the wall plume model methodology described to select the temperature difference distribution parameters for the best fit, the power n and the corresponding coefficient N are determined. For the four cases of FE1, FE2, FE3, and FE4, the optimal power law parameters (n, N) from the wall plume model were (-0.195, 2.526), (-0.26, 1.56), (-0.193, 1.64), and (-0.092, 1.67), respectively. The mass flow rates from the wall model were calculated using the obtained powers and corresponding coefficients and are compared with the mass flow rates obtained from CFD simulations. Figures 8a-8d show the variation of wall plume flow rates obtained from the model and from the two-dimensional enclosure simulations as a function of elevation for cases FE1 (a), FE2 (b), FE3 (c), and FE4 (d). The maximum percentage error between the predicted flow rates from the wall plume model and the two-dimensional enclosure simulations are less than 10% in all the four cases for elevations Z > [Z.sub.o]. The wall plume model predictions of flow rate are far more accurate than the isothermal wall predictions at the average temperature difference. The accuracy improves over the isothermal wall assumption with increased wall-to-air temperature difference, which is apparent when comparing case FE1 to case FE4 in Figures 8a and 8d.

The wall plume model has been shown to have an acceptable accuracy in predicting flow rates adjacent to nonuniformly heated walls in a two-dimensional enclosure with chilled ceiling. The next step is to implement the wall plume model in the space multilayer model to predict the DV/CC space vertical temperature gradients, stratification height, and chilled ceiling load.

Plume-Multilayer Model

The space is divided into four adjacent horizontal air layers with adjacent wall layers as shown schematically in Figures 9a and 9b. Each layer is represented by a uniform temperature T(i), where i corresponds to the air layer index (i = 1 to 4). Each layer interacts with side wall segments. The wall segment is represented by two nodal temperatures, the outer and the inner nodes. The thermal resistance diagram for the wall segment is shown in Figure 9c. The chilled ceiling temperature [T.sub.c] is constant and uniform. The fresh supply air enters from the side wall into the floor air layer at flow rate [M.sub.s] and temperature [T.sub.s] and the exhaust air leaves from the side wall of the upper air layer at [T.sub.(4)]. The wall temperature is designated by [T.sub.w](i,j), where j corresponds to the wall orientation (j = 1 to 4, where 1 = east, 2 = west, 3 = north, and 4 = south). The outside convection coefficients are [h.sub.ocv,f] and [h.sub.ocv,(i,j)] for the floor and walls, respectively.

The flow rates resulting from walls and heat sources are denoted by [M.sub.w(j,k)] and [M.sub.p(l,k)], respectively. The net circulated mass [M.sub.cir] at each boundary (interface surface between two adjacent air layers i and i + 1) will be calculated. The enclosure mass balance equation is given by

[M.sub.cir,k] = [M.sub.s] [r.summation over (l-1)] [M.sub.p,(l,k)] [4.summation over (j-1)] [M.sub.w,(l,k)], (8a)

where the wall plume mass flow rate of any wall j is given by

[M.sub.w,(j,k)] = [rho][Q.sub.w,(j,k)][Y.sub.(j)]. (8b)

In Equation 8, [rho] is the air density in kg/[m.sup.3] (lbm/[ft.sup.3]), [Y.sub.(j)] is the width of wall j, k refers to the interface level between the consecutive air layers (k = 1, 2, or 3 for a space divided into four layers), and r is the number of heat sources in the first air layer. If the circulation mass is positive, then the circulated mass is moving upward; if the circulated mass is negative, then it is moving downward. The determination of this circulated mass is very important because the stratification height is defined at the level where the circulated mass is zero.

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

The energy balance is then written for each horizontal air layer and for each wall segment. The first air layer energy balance is given by

[r.summation over (l=1)][[PHI].sub.(l,1)] + [4.summation over (i=1)] [h.sub.cv,(i,1)][A.sub.w,(i,1)] ([T.sub.w,(i,1)] - [T.sub.a,(1)]) + [h.sub.cv,f]([T.sub.f] - [T.sub.a,(1)]) + [M.sub.s][H.sub.s] = [r.summation over (j=1)] [M.sub.p,(j,1)][H.sub.p,(j,1)] + [4.summation over (i=1)] [M.sub.w,(i,1)][H.sub.w,(i,1)] + [M.sub.cir,1][H.sub.cir,1], (9)

where [PHI] is the heat source flux and H is the enthalpy in J/kg (Btu/lbm). For the left-hand side of Equation 9, the first term is the summation of the heat fluxes in the first layer, the second and third terms are heat convection from the walls and the floor, and the fourth term is the enthalpy of the supply air. For the right-hand side of Equation 9, the first and second terms are the enthalpies of the heat sources and wall plumes, and the third term is the enthalpy of the circulated mass at the boundary. The adopted correlations for laminar and turbulent convective heat transfer coefficients for vertical walls are, respectively:

[h.sub.cv] = 1.42 ([DELTA]T/Z)[.sup.0.25] (laminar convection with Z dimension in meters) (10a)

[h.sub.cv] = 1.31([DELTA]T)[.sup.0.33] (turbulent convection) (10b)

while the adopted correlation for the floor or ceiling is that of Min et al. (1956)given by

[h.sub.cv,f] = 2.13([DELTA]T)[.sup.0.31] (10c)

The value used by Rees (1998) was 5.9 W/[m.sup.2] x K (1.04 Btu/h x [ft.sup.2] x [degrees]F) and is close to the correlation of Min et al. (1956). The energy balance equations on the two intermediate layers and the upper layer can be similarly derived. It is also worthy to note that the number of horizontal layers can be increased from four to any number without any restrictions. This will make the model more powerful when dealing with more complicated geometries (presence of windows) and furniture settings.

The energy balance on the outer wall node gives

[h.sub.cv,o]([T.sub.out] - [T.sub.so]) + [q.sub.solar] = [U.sub.w]([T.sub.so(i,j)] [T.sub.si,(i,j)]), (11)

where [U.sub.w] is the wall conductance and [q.sub.solar] is the solar flux into the wall. The energy balance on the wall inner node yields

[U.sub.w]([T.sub.so(i,j)] - [T.sub.si,(i,j)]) +[q.sub.l] = [q.sub.cv] + [q.sub.rd], (12)

where [q.sub.l] is the radiative heat flux from light loads, [q.sub.cv] is the convective heat flux, and [q.sub.rd] is the radiative heat flux. The convective heat flux [q.sub.cv] on each part of the wall is given by

[q.sub.cv] = [h.sub.cv]([T.sub.si] - [T.sub.a]). (13)

The radiative heat flux [q.sub.rd,i] on a surface of index i is given by

[q.sub.rd,i] = [J.summation over (j=1)][epsilon][sigma]([T.sub.i.sup.4] - [T.sub.j.sup.4] [F.sub.i] [right arrow]j, (14)

where J is the number of surfaces in the space.

Radiation between the ceiling panel and other surfaces in the enclosure is an important parameter for heat exchange. The floor and the four walls of the enclosure are considered to participate in the radiant heat transfer as black bodies. However, the air within the enclosure is assumed to be a nonparticipating medium.

Modeling of heat source plumes in the presence of temperature gradients will follow the expressions given by Mundt (1996). The plume flow rate through the layer boundary depends on the plume source strength and vertical temperature gradient. In a region of uniform temperature, the volume flow rate [Q.sub.p] at an elevation Z above the ground level of the plume source with convective heat flux [PHI] is given by Mundt (1996) as

[Q.sub.p] = 0.005[[PHI].sup.1/3] (Z - [Z.sub.p])[.sup.5/3], (15)

where [Z.sub.p] is the source height from the ground level. Plumes are influenced by temperature stratification. The flow rate of a plume generated by a heat source of flux [PHI] in a space with temperature gradient [dT.sub.a]/dZ can be calculated by the following steps (Mundt 1996):

[Q.sub.p] = 0.00238[[PHI].sup.3/4]([dT.sub.a]/dZ)[.sup.-5/8] [B.sub.1] (16a)

[B.sub.1] = 0.004 + 0.039[A.sub.1] + 0.38[A.sub.1.sup.2] - 0.062[A.sub.1.sup.3] (16b)

[A.sub.1] = 2.86(Z - [Z.sub.p])([dT.sub.a]/dZ)[.sup.3/8] [[PHI].sup.-1/4] (16c)

In the presence of a temperature gradient [dT.sub.a]/dZ, the equilibrium height ([Z.sub.eq]), where the temperature difference between the plume and the ambient air disappears, is given by Mundt (1996) as

[Z.sub.eq] = 0.74[PHI][.sup.1/4]([dT.sub.a]/dz)[.sup.-3/8] (17)

The plume equilibrium heights for a point source of 100 W (341 Btu/h) are 2.34 m (7.7 ft) and 1.55 m (5 ft) for air temperature gradients of 1 K/m (0.55[degrees]F/ft) and 3 K/m (1.65[degrees]F/ft), respectively.

To solve the energy balance equations of the enclosure air layers, the values of the airflow rates passing from one layer to the other due to the generated plumes from the heat sources and the walls are calculated from the plume equation for the sources and using the wall plume model tested in this work, respectively.

The developed plume-multilayer mathematical model can now be translated into an efficient numerical model that solves for the vertical temperature distribution in a room conditioned by a combined DV/CC system in the presence of heat sources and nonuniform room wall conditions. The energy and mass balances represent a set of coupled nonlinear equations that are iteratively solved to obtain the layers and wall segment temperatures.

The numerical solution starts by assuming a temperature distribution for air layers and wall segments for a given temperature difference profile. The program proceeds by calculating the heat source plumes and the wall plumes. The set of algebraic equations resulting from the energy balance on the air layers and the walls were solved and new values of temperatures are calculated. The temperatures are updated and this procedure is repeated until convergence is reached when the change in values of temperature is less than 0.001[degrees]C (0.0018[degrees]F). A curve fit is then made for the wall-to-air temperature difference. If this curve fit matches the one assumed initially, then convergence is reached; otherwise, the new curve fit is used and the procedure is repeated until values of wall and air temperatures converge. The accuracy of the plume-multilayer model was tested by comparison to three-dimensional CFD simulations of an enclosure with heat sources using Airpak (2002) software interfaced with Fluent (2004) software.

PLUME-MULTILAYER MODEL TESTING

The validation of the complete modified plume-multilayer model is made by comparing the model outputs with values obtained from three-dimensional simulations of an enclosure conditioned by a DV/CC system using Airpak (2002) software, which was specifically developed for indoor airflow simulations. Chen and Xu (1998) tested the accuracy and reliability of a zero-equation turbulence model for DV-cooled rooms by experimentation and comparison to known CFD models. The agreement between the computed velocity and temperature vertical distributions was reported by Chen and Xu (1998) to be reasonably good. For that reason, the zero-equation turbulence model was adopted in the Airpak three-dimensional simulations of the test case comparison with the plume-multilayer model. The comparison between both models were based on the vertical temperature gradient, wall temperatures, cooling panel load, and stratification height. The simulations were done on a single geometry for three different heat loads (low, medium, and high).

The enclosure shown in Figure 10 has dimensions of 5 x 4 x 3 m (16.4 x 13.1 x 9.8 ft). Air is introduced into the room through a 0.6 x 1 m (1.96 x 3.28 ft) grille positioned at an elevation of 0.1 m (0.328 ft) above the ground and at a distance of 2 m (6.56 ft) from the edge of the south wall. The supply mass flow rate temperature is 22[degrees]C (71.5[degrees]F). The return air leaves the room through the exhaust grille, which is placed at the north wall, 2 m (6.56 ft) from its edge and 0.1 m (0.328 ft) below the ceiling level. The internal heat sources have a rectangular shape of 0.1 x 0.1 m (0.328 x 0.328 ft) in area and are placed at a height of 0.2 m (0.66 ft) above the floor. The heat flux of the internal heat sources is set to be convective only. The external heat load is conducted through the walls from outside conditions at fixed Tout and zero solar gain. The wall conductance is set to 1.53 W/K x [m.sup.2] (0.27 Btu/h x [ft.sup.2] x [degrees]F). The outside heat transfer coefficient is 21.1 W/K x [m.sup.2] (0.372 Btu/h x [ft.sup.2] x [degrees]F) for the walls and the floor. The supply mass flow rate, ceiling temperature, heat sources, and outside air temperature varied between low, medium, and high load cases. The different parameters used in the various cases are provided in Table 1.

The three cases were chosen based on the design chart of Behne (1999). The supply volume flow rates per unit area for the low, medium, and high load cases were 22.5 [m.sup.3]/h x [m.sup.2], 30 [m.sup.3]/h x [m.sup.2], and 37.5 [m.sup.3]/h x [m.sup.2] (1.23 [ft.sup.3]/min x [ft.sup.2], 1.64 [ft.sup.3]/min x [ft.sup.2], and 2.05 [ft.sup.3]/min x [ft.sup.2]), respectively. The low-load case is a lower limit case for the use of the combined DV/CC system. For the medium-load case, Behne (1999) recommends the use of DV or combining DV with a chilled ceiling. For the high load case, only the combined DV/CC system is advised by Behne (1999). The problem is set as steady with the radiation model turned on. The view factors are calculated and stored by the Airpak (2002) software. The number of nodes selected for the simulations is 11,950 (21 x 22 x 22) after checking solution grid independence at 76,015 (40 x 40 x 40) nodes and 231,141 nodes (60 x 60 x 60). Average heat flux load on the ceiling panel for the medium-load case was 51.5 W/[m.sup.2] (16.3 Btu/h) on the coarse grid compared to 50.47 W/[m.sup.2] (15.32 Btu/h) for the finest grid result. The convergence criteria were 0.001 m/s for the flow and [10.sup.-6] W for energy. The under-relaxation factors were the default values: 0.7 for pressure, 0.3 for momentum, 1.0 for temperature and viscosity, and 0.1 for body forces. Convergence was reached after 200-250 iterations.

[FIGURE 10 OMITTED]

The plume-multilayer model was applied on the three test cases dividing the room into four air layers of equal height and using the offset height [Z.sub.o] equal to the height of the first layer. The developed procedure for determining n and N is applied to each case to iterate for a converged solution starting from an initial linear temperature distribution between the floor and the ceiling temperatures. The wall plume flow rates, the average temperature of each air layer, and the associated average wall segment temperature were then obtained from the converged solution. The values of the parameters (N, n) derived from the model for the low, medium, and high load test cases were (0.48, -0.16), (0.82, -0.2), and (1.61, -0.14), respectively. The vertical temperature gradient in the room is one of the main parameters that affects thermal comfort inside the room, while the stratification height affects indoor air quality.

The Airpak (2002) software calculates the temperature at each node in the thermal space. For the vertical temperature gradient comparison to be valid, an average temperature for each layer must be calculated from the Airpak results. The average air layer temperature in Airpak is calculated by averaging the temperature of all the grid points in the layer volume. Figure 11 shows the average air temperature of each air layer of the four layers calculated from the plume-multilayer model and Airpak for the (a) low, (b) medium, and (c) high load cases. The maximum differences in the temperature values are 0.13[degrees]C (0.23[degrees]F), 0.37[degrees]C (0.66[degrees]F), and 0.3[degrees]C (0.54[degrees]F) for the low, medium, and high load cases, respectively.

The second comparison of the model values is based on the wall temperature distribution. The calculation of the temperatures of the wall parts at the different layers is done by averaging the wall temperature values over the wall width at each level using Airpak (2002) simulation results. The values of wall temperatures for the east and west walls are the same due to the symmetry of the room. Figure 12 shows the variation of the average vertical wall temperature associated with air layers for all the walls as predicted by the plume-multilayer model and Airpak for the low, medium, and high load cases for the east/west and north/south temperatures. The maximum error in the east/west wall temperature is about 0.41[degrees]C (0.72[degrees]F) for the medium load case. For the north and south walls, the maximum wall temperature differs by 0.5[degrees]C (0.9[degrees]F) and 0.8[degrees]C (1.44[degrees]F), respectively, due to the existence of the supply grille on the south wall and the return grille on the north wall. For the medium load case, a maximum temperature difference of about 0.83[degrees]C (1.5[degrees]F) exists between the model and Airpak results at the lower levels and a temperature difference of about 0.69[degrees]C (1.24F) exists at the higher levels. In general, wall plumes are driven by the temperature difference between the wall and the room air. This difference is inherently not uniform in DV systems, and that effect is compounded when the walls are not well insulated, which could result in a larger temperature difference between the wall and the room air.

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

The panel cooling load is calculated from Airpak (2002) software by multiplying the resulting panel heat flux by the ceiling area. The panel cooling loads computed by the plume-multilayer model for low, medium, and high loads were 574 W (1,958 Btu/h), 905.4 W (3,089 Btu/h), and 1,430 W (4,878 Btu/h), respectively, compared to three-dimensional simulation results of 530 W (1,808 Btu/h), 1,022 W (3,486 Btu/h), and 1506.4 W (5,139 Btu/h), respectively. The maximum error is less than 11% between model predictions and three-dimensional simulations. The intention for accurate estimation of the flow rate is to be able to predict the stratification height, which is an important parameter for thermal comfort and air quality considerations. The model prediction for the stratification height inside the room is based only on mass balance equations. The stratification height was predicted by the model to be at a level of 1.372 m (4.5 ft), 1.16 m (3.8 ft), and 1.15 m (3.77 ft) for the low, medium, and high load cases, respectively.

The three-dimensional flow simulation inside the enclosure represents a case of mixed convection, while the present model assumes no interaction between sources and the wall. The stratification height is distinguished in the three-dimensional flow enclosure results by the highest level at which no noticeable downward flow exists. The upward and downward mass flow rates were calculated in the space at different elevations by integrating the vertical velocity component distribution over the horizontal surface area. The variation of the upward and downward mass flow rates with elevation for the low, medium, and high load cases is shown in Figures 13a, 13b, and 13c, respectively. The up and down mass flow rates remain constant below the stratification height region and then slightly change slope in the upper region. The stratification heights predicted by the model are shown directly in Figure 12, which indicates that their values fall closely in the region when a slight change in upward mass flow rate takes place. The model predictions are within 0.05 m (0.16 ft) of the three-dimensional model.

[FIGURE 13 OMITTED]

The importance of including wall plumes in DV with chilled ceiling models is seen by examining the stratification height and the vertical temperature gradient when wall plumes are neglected or isothermal plumes are used. Figure 14 shows the variation of (a) the stratification height and (b) the vertical temperature gradient, with supply flow rate at a fixed room perimeter-to-area ratio while varying the ceiling temperature to meet the steady-state load of the medium case. For a cooling load of 67.3W/[m.sup.2] (21.31 Btu/h x [ft.sup.2]), a zero wall-plume value overestimates the stratification height by 35% within the space, while an isothermal plume overestimates the stratification height by 8% to 10%. The vertical temperature gradient is decreased by 0.07 K/m (0.04[degrees]F/ft) when wall plumes are in the model compared to zero wall plume.

The model accurately predicted the stratification height when compared to the three-dimensional CFD model. The added complexity of calculating the wall plumes over isothermal correlation using the power law similarity model improved predictions of wall flow rates and stratification height to become comparable to three-dimensional full CFD simulations at a much lower computational cost.

CONCLUSIONS

A new, improved simple model was developed to represent room heat transfer with a combined DV/CC system. The model is based on energy conservation equations in air layers while deriving airflow rates exchanged between adjacent layers from plume equations of Mundt (1996) and from a new wall plume model developed in this study. The heat source induced flow rate relied on plume equations of point and area heat sources of Mundt (1996). The wall plume flow rates are based on a similarity solution of two-dimensional boundary layer flow with power law temperature difference distribution between the heated wall and the environment air. The determination of the power law parameters of the temperature difference distribution was subject to rigorous testing on the selected exponent and the temperature difference coefficient that predicts the wall plumes most closely to the real physical situation occurring in the space.

The plume-multilayer model using the developed wall plume model was validated in comparison with detailed three-dimensional simulation of the enclosure, and results of the average air and wall temperature variations with elevation, the stratification height, and the panel cooling load agreed well with the plume-multilayer model predictions. The model has very good accuracy at a low computational cost and with sufficient flexibility in the geometry, source strengths, and supply flow rates of real systems. The flexibility stems from the physical consideration of the mechanism of convection that takes place in the enclosures conditioned by DV and a chilled ceiling. Future work will extend the model to include turbulent plume regimes and to use the model in transient application when changes in the load or external environment take place.

ACKNOWLEDGMENTS

The authors acknowledge the support of the University Research Board of American University of Beirut research grant DDF-113040-688324, the Lebanese National Council for Scientific Research grant LCR-113040-002208, and ASHRAE grant 113040-025043.

[FIGURE 14 OMITTED]

NOMENCLATURE

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

CC = chilled ceiling

CFD = computational fluid dynamics

d = temperature difference between the wall and the air, [degrees]C ([degrees]F)

[d.sub.max] = maximum temperature difference between the wall and the air, [degrees]C ([degrees]F)

DV = displacement ventilation

f = dimensionless stream function

[F.sub.i[right arrow]j]= view factor from surface i to surface j

g = gravitational acceleration

[Gr.sub.Z] = Grashof number, [Gr.sub.Z] = g[Z.sup.3]/[v.sup.2] [beta]d(Z)

H = enthalpy, J/kg (Btu/1bm)

h = heat transfer coefficient, W/[m.sup.2] x K (Btu/h x [ft.sup.2] x [degrees]F)

J = number of surfaces in the space

L = room length, m

M = mass flow rate, kg/s (lbm/h)

n = power of the curve fit

N = power law coefficient

Pr = Prandtl number

q = heat flux, W/[m.sup.2] (Btu/h x [ft.sup.2])

[Q.sub.p] = plume volumetric flow rate, [m.sup.3]/s ([ft.sup.3]/min)

[Q.sub.w] = wall plume volumetric flow rate per unit width, [m.sup.3]/m x s ([ft.sup.3]/ft x min)

r = number of heat sources in the first air layer

Re = Reynolds number

[R.sup.2] = data fit correlation parameter

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

t = time, s

U = wall conductance, W/[m.sup.2] x K (Btu/h x [ft.sup.2] x [degrees]F)

u = horizontal component of the velocity, m/s (ft/min)

v = vertical component of the velocity, m/s (ft/min)

W = width of the room, m (ft)

[Y.sub.(j)] = width of wall j, m (ft)

Z = elevation above ground level, m (ft)

[Z.sub.o] = offset elevation, m (ft)

[Z.sub.p] = source height from floor level, m (ft)

[Z.sub.eq] = plume equilibrium height (ft)

Greek Symbols

[alpha] = thermal diffusivity, [m.sup.2]/s ([ft.sup.2/s]

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

[DELTA]T = difference between average or isothermal wall temperature and far fluid temperature [degrees]C ([degrees]F)

[PHI] = convective heat flux, W (Btu/h)

[phi] = dimensionless temperature, (T - [T.sub.a])/d(Z)

[eta] = transformation parameter, [eta] = y/z4[square root of ([Gr.sub.z]/4)]

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

[upsilon] = kinematic viscosity, [m.sup.2]/s ([ft.sup.2]/s)

Subscripts

a = air

c = ceiling

cir = circulation

cv = convection

f = floor

l = light

o = outer

p = plume

rd = radiation

s = supply

solar = solar radiation

si = inner surface wall temperature

so = outer surface wall temperature

w = wall

REFERENCES

Airpak. 2002. Airpak 2.1, Computational Fluid Dynamics Software for Air Conditioning Applications, Fluent, Inc., Lebanon, NH. www.fluent.com.

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

Behne, M. 1999. Indoor air quality in rooms with cooled ceilings: Mixing ventilation or rather displacement ventilation? Energy and Buildings 30(2):155-66.

Carrilho da Graca, G.C.C. 2003. Simplified models for heat transfer in rooms. PhD thesis, University of California, San Diego.

Chapra, S.C., and R.P. Canale. 2002. Numerical Methods for Engineers, with Programming and Software Applications, 4th ed. Boston: W.C. Brown/McGraw-Hill.

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

Chen, Q., and W. Xu. 1998. A zero-equation turbulence model for indoor airflow simulation. Energy and Buildings 28(2):137-44.

Etheridge, D., and M. Sandberg. 1996. Building Ventilation: Theory and Measurement. Chichester: John Wiley and Sons.

Fluent. 2004. Computational Fluid Dynamics Software, Fluent 6.2.16, Fluent, Inc., Lebanon, NH. www.fluent.com.

Gebhart, B., Y. Jaluria, R. Mahajan, and B. Sammakia. 1988. Buoyancy Induced Flows and Transport. New York: Hemisphere Publishing Corporation.

Goodfellow, H.D., and T. Esko, eds. 2002. Industrial Ventilation Design Guidebook. San Diego, CA: Academic Press.

Jaluria, Y. 1980. Natural Convection, Heat and Mass Transfer. New York: Pergamon Press.

Jiang, Z., Q. Chen, and A. Moser. 1992. Indoor airflow with cooling panel and radiative/convective heat source. ASHRAE Transactions 98(1):33-42.

Linden, P.F., G.F. Lane-Serff, and D.A. Smeed. 1990. Emptying filing boxes: The fluid mechanics of natural ventilation. Journal of Fluid Mechanics 212:300-35.

Min, T.C., L.F. Schutrum, G.V. Parmelee, and J.D. Vouris. 1956. Natural convection and radiation in a panel heated room. ASHRAE Transactions 62:337-58.

Mundt, E. 1996. The Performance of Displacement Ventilation Systems: Experimental and Theoretical Studies. Stockholm, Sweden: Royal Institute of Technology.

Novoselac, A., and J. Srebric. 2002. A critical review of the performance and design of combined cooled ceiling and displacement ventilation systems. Energy and Buildings 34:497-509.

Rees, S.J. 1998. Modeling of displacement ventilation and chilled ceiling systems using nodal models. PhD thesis, Loughborough University, London, UK.

Rees, S.J., J.J. McGuirk, and P. Haves. 2001. Numerical investigation of transient buoyant flow in a room with a displacement ventilation and chilled ceiling system. International Journal of Heat and Mass Transfer 44(16):3067-80.

Skaaret, E. 1998. A semi-empirical flow model for low-velocity air supply in displacement ventilation. Proceedings of Roomvent 98, Stockholm, Sweden 1:85-92.

Yuan, X., Q. Chen, and L. Glicksman. 1998. A critical review of displacement ventilation. ASHRAE Transactions 104(1):78-90.

Mohamad Ayoub

Nesreen Ghaddar, PhD

Member ASHRAE

Kamel Ghali, PhD

Received February 8, 2006; accepted April 28, 2006

Mohamad Ayoub is a graduate research assistant and Nesreen Ghaddar is a professor and Qatar Chair for Energy Studies in the Department of Mechanical Engineering, American University of Beirut, Beirut, Lebanon. Kamel Ghali is an associate professor in the Department of Mechanical Engineering, Beirut Arab University, Beirut, Lebanon.

Table 1. Input Parameters of the Test Cases Case Low Load Medium Load Mass flow rate, kg/s (lbm/h) 0.15 (1,188) 0.2 (1,587) Ceiling temperature, 19.0 (66.2) 18.0 (64.4) [degrees]C ([degrees]F) Internal heat sources, W (Btu/h) 2 x 100 4 x100 (2 x 341) (4 x 341) Outside air temperature 28.0 (82.4) 34.0 (93.2) [T.sub.out], [degrees]C ([degrees]F) Cooling load, W/[m.sup.2] 39.5 (12.5) 67.37 (21.35) (Btu/h x [ft.sup.2]) Case High Load Mass flow rate, kg/s (lbm/h) 0.25 (1,984) Ceiling temperature, 17.0 (62.6) [degrees]C ([degrees]F) Internal heat sources, W (Btu/h) 4 x 250 (4 x 852.5) Outside air temperature 34.0 (93.2) [T.sub.out], [degrees]C ([degrees]F) Cooling load, W/[m.sup.2] 98.21 (31.06) (Btu/h x [ft.sup.2])

Printer friendly Cite/link Email Feedback | |

Author: | Ayoub, Mohamad; Ghaddar, Nesreen; Ghali, Kamel |
---|---|

Publication: | HVAC & R Research |

Geographic Code: | 1USA |

Date: | Oct 1, 2006 |

Words: | 10689 |

Previous Article: | Simulated performance of natural and hybrid ventilation systems in an office building. |

Next Article: | Silicon microchannel cooling for high power chips. |

Topics: |