Evaluation of various turbulence models in predicting airflow and turbulence in enclosed environments by CFD: part 2--comparison with experimental data from literature.

INTRODUCTION

The companion paper (Zhai et al. 2007) reviewed the recent development and applications of computational fluid dynamics (CFD) approaches and turbulence models for predicting air motion in enclosed spaces. The review identified eight prevalent and/or recently proposed turbulence models for indoor airflow prediction. These models include the indoor zero-equation model (0-eq.) by Chen and Xu (1998), the RNG k-[epsilon] model by Yakhot and Orszag (1986), a low Reynolds number k-[epsilon] model (LRN-LS) by Launder and Sharma (1974), the SST k-[omega] model (SST) by Menter (1994), a modified v2f model (v2f-dav) by Davidson et al. (2003), a Reynolds-stress model (RSM-IP) by Gibson and Launder (1978), the large eddy simulation (LES) with a dynamic subgrid scale model (LES-Dyn) (Germano et al. 1991; Lilly 1992), and the detached eddy simulation (DES-SA) by Shur et al. (1999). This paper evaluates and compares the selected turbulence models for several indoor benchmark cases that represent the primary flow mechanism of air movement in enclosed environments.

All of the turbulence model equations mentioned above can be written in a general form as follows:

[rho][[[partial derivative][[bar].[phi]]/[[partial derivative]t]] + [[rho][bar].[u.sub.j]][[[partial derivative][bar].[phi]]/[[partial derivative][x.sub.j]]][[partial derivative]/[[partial derivative][x.sub.j]]][[[GAMMA].sub.[phi],eff][[partial derivative][bar.[phi]]]/[[partial derivative][x.sub.j]]] = [S.sub.[phi]] (1)

where [phi] represents variables, [[GAMMA].sub.[phi],eff] represents the effective diffusion coefficient, and [S.sub.[phi]] represents the source term of an equation. Table 1 briefly summarizes the mathematical expressions of the eight turbulence models selected. In Table 1, [u.sub.i] is the velocity component in i direction, T is the air temperature, k is the kinetic energy of turbulence, [epsilon] is the dissipation rate of turbulent kinetic energy, [omega] is the specific dissipation rate of turbulent kinetic energy, P is the air pressure, H is the air enthalpy, [[mu].sub.t] is the eddy viscosity,[G.sub.[phi]] is the turbulence production for, [phi] and is the rate of the strain. The other coefficients are case-specific, and only some of those introduced here are important.

[TABLE 1 OMITTED]

For the 0-eq. model, V is the velocity magnitude and l is the wall distance. The [G.sub.B] is the buoyancy production term for the RNG k-[epsilon] model. For the LRN-LS model,[f.sub.[mu]], [C*.sub.[epsilon]1], and, [C*.sub.[epsilon]2] are the three modified coefficients (i.e., damping functions) to the standard k-[epsilon] model, and D and E are two additional terms. These five major modifications to the LRN model are responsible for improving model performance near the wall. In the SST model, Y is the dissipation term in the k and [omega] equations, [F.sub.1] and [F.sub.2] are blending functions that control the switch between the transformed k-[epsilon] model and the standard k-[omega] model, and [D.sub.[omega]] is produced from the transformed k-[epsilon] model, so it vanishes in the k-[omega] model when the blending function [F.sub.1] equals one. In the v2f-dav model (Davidson et al. 2003),[[bar].v'.sup.2] is the fluctuating velocity normal to the nearest wall. The variable f is part of the [[bar].v'.sup.2] source term that accounts for nonlocal blocking of the wall normal stress. The variable f is implicitly expressed by an elliptical partial differential equation. So the scalar f in principle can be solved by the same partial-differential-equation solver as the other variables. Note that T in the v2f-dav model also represents the turbulence time scale. In the RSM model, the [[phi].sub.lm] is the pressure-strain term and requires further modeling. In the present study, a liner pressure-strain model by Gibson and Launder (1978) is used.

In the LES model, the overbar represents the filtering. Expressions [[tau].sub.ij.sup.S] and [h.sub.j.sup.S] represent the subgrid-scale (SGS) stress and heat flux. Lilly's (1992) SGS model adopts the Boussinesq hypothesis and derives methods to calculate the coefficient, [C.sub.S], in the eddy-viscosity expression automatically. The presented DES model (Shur et al. 1999) couples the LES model with a one-equation RANS model (Spalart and Allmaras 1992). This one-equation model directly solves a modified eddy viscosity rather than the turbulence kinetic energy as most one-equation models do. The d variable is wall distance and [f.sub.v1] and [f.sub.v2] are the damping functions. Due to space constraints for this paper, a more detailed description of these models is not possible. Since many of the models are available in some commercial software, one could also refer to a user manual (e.g., FLUENT [2005]) for detailed model descriptions.

NUMERICAL METHOD

This study used commercial CFD software, FLUENT version 6.2 (FLUENT 2005) to conduct all the numerical investigations discussed in the next section. Most of the models shown in Table 1 are available in FLUENT except for the modified v2f-dav model. We applied user-defined scalar (UDS) transport equations and coded user-defined functions (UDF) to describe the governing equations of the k, [epsilon], and [[bar].v'.sup.2], as well as the elliptical partial differential equation for f. The RANS models used the second-order upwind scheme for all of the variables except pressure. The discretization of pressure is based on a staggered scheme, PRESTO! (FLUENT 2005). The SIMPLE algorithm was adopted to couple the pressure and momentum equations. If the sum of absolute normalized residuals for all of the cells in flow domain became less than [10.sup.-6] for energy and [10.sup.-3] for other variables, the solution was considered converged. Grid dependence of each case was checked using two to four different grids to ensure that grid resolution would not have a notable impact on the results.

RESULTS AND ANALYSIS

This study evaluated the performance of the eight selected models by simulating the distributions of airflow, air temperature, and turbulence quantities in four different enclosed environments. The four cases under investigation are natural convection in a tall cavity, forced convection in a model room with partitions, mixed convection in a square cavity, and strong buoyancy flow in a model fire room. The first three cases are benchmark cases that represent the most basic flow features in an enclosed environment. The fourth case is more challenging and can be used to test model robustness in the present research. Figure 1 shows the geometric and airflow information of the four cases. Detailed comparison and result analyses are discussed in the following subsections.

[FIGURE 1 OMITTED]

Natural Convection in a Tall Cavity

Natural convection in an enclosed environment is attributed to the buoyancy effect caused by the existence of gravity and fluid density differential. Typical examples of natural convection in an enclosed environment include thermal plumes generated by heat sources, and airflow near a wall or a window generated by a temperature difference between the surface and the air, etc.

Betts and Bokhari (2000) conducted an experimental investigation of natural convection in a tall cavity, as shown in Figure 1a. The dimensions of the cavity were 2.18 m high x 0.076 m wide x 0.52 m deep. The cold and hot walls had uniform temperatures of 15.1[degrees]C and 34.7[degrees]C, respectively. The Rayleigh number based on the cavity width was 0.86 x [10.sup.6]. The large ratio of cavity height and width ensured that the airflow in the core region was fully turbulent despite a relatively low Rayleigh number. Their experimental data also showed that the airflow pattern was approximately two-dimensional in the vicinity of the center plane. The air temperature inside the cavity was measured by a thermocouple. The air velocity was measured by a single component laser-Doppler anemometry (LDA) system.

The numerical results presented were based on a 25 x 150 nonuniform two-dimensional grid for all RANS models except the LRN-LS and a 25 x 150 x 50 three-dimensional grid for LES and DES. The corresponding [y.sup.+] was about 0.3 for the first grid close to the walls, while the [y.sup.+] for LRN-LS grid was less than 0.1. Since the calculated [y.sup.+] was rather small, the enhanced wall treatment (FLUENT 2005) was adopted for RNG k-[epsilon] and RSM-IP models. The same treatment was also adopted in the other cases.

Figure 2 compares the simulated results with the measured data. The zero-equation model produced significant errors on the mean air velocity from y/H = 0.3 to 0.7, although the predicted air-temperature profiles seem acceptable. The LRN-LS model could not correctly predict the air-temperature profile near the top and bottom walls. The DES model did not perform well. For velocity prediction, it had similar accuracy compared to the indoor zero-equation model. The present DES adopted the S-A one-equation model near the walls, and the inaccuracy of DES results was mainly associated with the S-A model performance. While the results from the other models reasonably agreed with the experimental data for temperature and vertical velocity, the v2f-dav model exhibited the best agreement.

[FIGURE 2 OMITTED]

For the normal Reynolds stress, the v2f-dav and the LES results best agreed with the measurements. The other models only predicted a similar Reynolds-stress profile but failed to give the correct magnitudes. The normal component of Reynolds stress,[[bar].v'.sup.2], can be written as shown in Equation 2.

[[bar].v'.sup.2] = [[2/3][rho]k] - 2[[micro].sub.t][[[partial derivative].sub.v]/[[partial derivative].sub.y]] (2)

The product of the turbulent viscosity and the partial derivative of normal velocity is negligible compared to the turbulence kinetic energy. The underprediction of [[bar].v'.sup.2] for the other models is due to the underprediction of the turbulence kinetic energy. For all the selected RANS models except the zero-equation model, the k equation has a very similar form. The source of the k equation includes three terms: production, dissipation, and buoyancy production. The dissipation term is generally smaller than the production in this case, unless very close to wall. The buoyancy production term, [G.sub.B] = ([[mu].sub.t]/[[sigma].sub.T, t])g/T[[partial derivative]T/[partial derivative]y], is about one-to-two orders smaller than that of the turbulence production, which is in the following form:

[G.sub.k] = [[mu].sub.t]{2[([[partial derivative]u]/[[partial derivative]x]).sup.2] + [([[partial derivative]v]/[[partial derivative]y]).sup.2]] + [([[partial derivative]v]/[[partial derivative]x] + [[partial derivative]u]/[[partial derivative]y]).sup.2]} (3)

Clearly, the [[partial derivative]v/[partial derivative]x] is much larger than other partial derivatives in this particular case. Therefore, the inaccurate prediction [[partial derivative]v/[partial derivative]x] of influenced the prediction of the turbulence kinetic energy and led to the underprediction of the normal Reynolds stress. As shown in Figure 2b, the RNG k-[epsilon], LRN, and RSM-IP models all predicted a lower [[partial derivative]v/[partial derivative]x] in the center region of the cavity compared to experimental data. So these models may underpredict the overall production of kinetic energy, leading to the inaccuracy of the simulated turbulent stress.

Furthermore, for all of the eddy-viscosity models, the momentum equations have the same form, while the only difference is the expression of the eddy viscosity. Except for the v2f model, the other eddy-viscosity models used [square root of]k as the velocity scale in determining the eddy viscosity. The inaccuracy of predicted k, therefore, adversely affected the prediction of the mean velocity field and the eddy viscosity. In the present case, the v2f-dav model provided a more suitable eddy-viscosity expression and, thus, achieved a better accuracy.

Forced Convection in a Room with Partitions

Forced convection is often encountered in enclosed spaces with mechanical ventilation systems. Air jets coming out of diffusers and airflow from fans are typical examples of forced convection in rooms. This study used a forced-convection case with experimental data from Ito et al. (2000) to analyze the performance of the turbulence models. Figure 1b shows the sketch and dimensions of the cross-sectional view of the room. The air-supply diffuser was located at the upper-left corner with a height of 0.02 m in y direction. The exhaust outlet was at the upper-right corner with the same size as the inlet. Two partitions, each 0.5 m high, were located at the lower part of the room.

The room airflow was isothermal and the temperature was about 25[degrees]C. The mean and turbulent velocities were measured by a two-component laser-Doppler velocimeter (LDV). The measured data showed that the mean air velocity did not vary much in the spanwise direction and the airflow was approximately two-dimensional. The mean air supply velocity was 3.0 m/s with a turbulent intensity of 1.6%. The Reynolds number based on the inlet condition was about 4000. Although the room airflow was not laminar, the turbulence level within the domain was very low. Figure 1b also shows the measured mean velocity field. The main airstream was attached to the ceiling and traveled down to Zone 3. The main circulation was clockwise. The secondary counter-clockwise circulations were observed in Zones 1 and 2.

Figures 3 and 4 show the detailed comparison of the model predictions with the experimental data of velocity and turbulent quantities on the vertical and horizontal center lines labeled as red dash-dotted lines in Figure 1b. All the turbulence models could accurately predict the air-jet flow pattern near the ceiling. For the prediction of various circulations in the room, the model performance varies. The indoor zero-equation model predicted a reversed U velocity profile in Zone 2, as shown in Figure 4. Similarly, the SST k-[omega] model predicted a reversed V velocity profile in Zone 1, indicating a wrong air circulation pattern. Meanwhile, the SST k-[omega] model produced the least satisfactory agreement with the experimental data compared to the other models. In principle, the SST k-[omega] model should have similar performance as the k-[epsilon] model in regions far away from the walls. The blending function,[F.sub.1], in Table 1 was designed to automatically switch to the k-[epsilon] model formulation outside the boundary layers where [F.sub.1] should vanish. Figure 5 shows the computed value of [F.sub.1] in the domain. The [F.sub.1] value was large even far away from the walls. So the SST k-[omega] model essentially worked as a model rather than a model in this case. In fact, the SST k-[omega] model used limiting functions to ensure blending functions vanished outside the wall region (Menter 1994). Those limiting functions have been tested as valid for many turbulent flows. Nevertheless, the turbulence level in room air is generally very low, and those limiting functions may not be always valid. It is necessary to modify the blending functions in order to improve the accuracy of the SST k-[omega] model in such low-turbulence airflow. The DES model underpredicted the vertical velocity fluctuation and resolved Reynolds stress, as shown in Figure 3. Again, the errors may come from the RANS form of the DES model, although no clear evidence can be provided at present. Further investigation of this and other DES models may improve our understanding of the inherent problem. Besides the zero equation and the SST k-[omega] models, the other models have similar accuracy, although some small disparities exist.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

[FIGURE 5 OMITTED]

Mixed Convection in a Square Cavity

Mixed convection is the most common airflow form in an air-conditioned environment. Blay et al. (1992) studied a mixed-convection flow in a square cavity using both experimental and computational methods. Shown in Figure 1c, air was discharged from the inlet slot at the ceiling level and exhausted at the floor level on the opposite wall, and the floor was heated. The measured inlet conditions were [u.sub.in] = 0.57 m/s;[v.sub.in] = 0; [T.sub.in] = 15[degrees]C;[[epsilon].sub.in] = 0; and [k.sub.in] = 1.25 x [10.sup.-3] [m.sup.2]/[s.sup.2]. The wall temperature, [T.sub.w], was 15[degrees]C and the floor temperature,[T.sub.fl], was 35.5[degrees]C. The Reynolds number based on the inlet condition was 684.

All the RANS models were simulated two-dimensionally. With a grid resolution of 60x60 in the two-dimensional domain, the [y.sup.+] for the first grid was about 1. For LES and DES, 30 uniform grids were placed in the spanwise direction, while the grid resolution in cross section was the same as that for RANS. Figure 6 presents the numerical results.

[FIGURE 6 OMITTED]

In general, all of the numerical simulations agreed reasonably with the experimental data for the air-temperature profile. For the temperature prediction, the LES and the v2f-dav agreed better with the measured data than the other models. Most of the models can also calculate the turbulence kinetic energy fairly well except for the SST k-[omega] model. Overall, the SST model predicted the turbulence kinetic energy to be 50% lower than the measurement, while this result was similar to that of a standard k-[omega] model (the results not shown here). As discussed in the forced convection case, the SST model might not switch to the k-[epsilon] model in regions far from the walls when the flow turbulence level is relatively low. Special care must be taken to apply the SST model in such flow regimes, while some modifications on the model blending functions may be needed.

Strong Natural Convection in a Model Fire Room

The three cases studied above represent the typical flow mechanisms in enclosed environments. Some models performed reasonably well for both cases, while the others did not. Another case with extreme buoyancy conditions was employed to further test the robustness of those models in a more challenging scenario: a model fire room with strong buoyancy flow. This case was designed by Murakami et al. (1995), who measured detailed mean and turbulence quantities. The chamber dimensions were 1.8 m long x 1.2 m wide x 1.2 m high, as shown in Figure 1d. The total heat power input from the heat sources was 9.1 kW, with an average surface temperature higher than 500[degrees]C. The opening size was 0.4 m wide x 0.9 m high. The air flowed through the opening between the chamber and its outside enclosure. The outer enclosure was about 8000 [m.sup.3]. All of the walls of the chamber were well insulated. Velocity vectors were measured by a two-component LDV. The air and wall temperatures were measured by thermocouples.

For the numerical simulations, this study included a part of the outside enclosure into the computation domain. The inclusion can avoid the numerical instability by setting pressure boundary condition for the opening of the model room, so the opening can be treated as an interior one. In order to save computing time, the grids in the included outside enclosure were rather coarse. But the grids inside the chamber remained reasonably fine (70 [L] x 55 [W] x 60 [H] = 231,000 cells) with an averaged [y.sup.+] = 1 for the first grids near the walls. In numerical simulation, the outer enclosure size was set to 10 m (X) x 20 m (Y) x 5 m (Z). Including the grids for the outer enclosure, the total grid number was 350,000. Since the grid in the enclosure is too coarse for an LES, the LES did not include the outside enclosure part.

Inside the chamber, the temperature difference can be on the order of several hundred Kelvin. The large temperature gradient caused significant air-density change, while the flow was still incompressible (low Ma number). Thus, the density variation needs to be considered by using the ideal gas law rather than the Boussinesq approximation as follows:

[rho] = [[p.sub.op]/[RT]] (4)

Where [P.sub.op] is the operating pressure in the room (here, 1 atm at sea level),R is the gas constant, and T is the air temperature.

The comparison between the numerical and experimental results in four measurement locations is shown in Figure 7. The RSM model could not converge with the grid distribution and boundary conditions used. This is likely due to the coarse grid used in the outer space. Thus, its results are not presented in the figure. For the mean temperature and velocity, all of the predictions presented agreed reasonably well with the measured data. Note that the buoyancy production term was added to the v2f-dav model. Due to the strong buoyancy effect in the chamber, the buoyancy production term is not negligible in the source terms of the turbulent kinetic energy equation.

[FIGURE 7 OMITTED]

Generally, all of the models gave a very good prediction of the mean temperature and velocity. The models, however, display obviously different performance in predicting the velocity fluctuation. The RNG k-[epsilon] model underpredicted the fluctuating velocity in the lower part of the room. By examining Equation 2 for the expression of [[bar].v'.sup.2], the [rho]k term, in this case, is at least one order of magnitude larger than the other on the right-hand side of Equation 2. Similar to the natural convection case in the tall cavity, the underprediction of [[bar].v'.sup.2] is likely due to the underprediction of k. The v2f and the RNG k-[epsilon] models have the exact same form of the k equation except for the expression of eddy viscosity, while the v2f-dav model has better accuracy. Thus, the inaccurate [[bar].v'.sup.2] prediction is possibly associated with the eddy-viscosity formulation of the k-[epsilon] model as well as the near-wall treatment. In addition to the v2f-dav model, the SST and LRN models also have much better performance than the high Reynolds number k-[epsilon] model. The SST model could successfully switch between the k-[omega] and k-[epsilon] models, as the airflow in the model room is more turbulent compared to the previous two cases. The near-wall treatment of a k-[omega] formulation worked better than that of a high Reynolds number k-[epsilon] model in this case. The present low Reynolds modification (Launder and Sharma 1974) also worked reasonably well and improved the performance of high Reynolds number k-[epsilon] models.

The LES results agreed with the measurement well in the center of the room but were less accurate near the wall. The errors near the walls are likely due to the limited ability of the subgrid scale model used here (Lilly 1992). The SGS model is essentially an algebraic model for subgrid scale eddies. The algebraic nature of the SGS model affects the LES accuracy close to the wall, especially when the turbulence is not locally in equilibrium. Therefore, the LES could be less accurate than some of the eddy-viscosity models.

The accuracy of the zero-equation model in this case was remarkable. Its results are comparable to those more advanced models. In addition, the zero-equation model was fastest and most stable. In fact, the other RANS models tested used the flow and temperature field calculated by the zero-equation model as initial conditions. Otherwise, the convergence of some models would be very difficult or even impossible.

DISCUSSION

The model accuracy has been analyzed by comparing the predicted results with the experiment data. In addition to accuracy, the computing cost is another important aspect that relates to model performance.

All of the RANS and DES simulations for the first three cases were conducted on a personal computer with a Pentium IV 3.0 GHz CPU with 1 GB memory. The RANS simulations for the fourth case were conducted in a computer cluster of three dual-CPU nodes with 2GHz CPU speed and 3.7 GB memory for each node. All of the LES simulations (except the fourth case) were performed on a computer cluster with four nodes, and each node had an AMD Opteron (64-bit) 2.6GHz CPU with 1 GB memory. Generally, four factors influenced the computing time: (1) grid resolution, (2) discretization scheme, (3) degree of nonlinearity of the model, and (4) number of PDEs the model contains. By fixing the first two factors, the difference in computing time is mainly attributed to the turbulence model itself. Taking the two-dimensional mixed-convection case, for example, all simulations were based on the same grid resolution (60 x 60). The indoor zero-equation model required two minutes of computing time, the RNG k-[epsilon] model seven minutes, the SST k-[omega] model eight minutes, the v2f-dav model around 13 minutes, the LRN-LS 15 minutes, and the RSM-IP 35 minutes, respectively. The DES used a 14-day computing time to finish the calculation of 20,000 time steps of 0.01 s each. The LES took eight days to finish the transient calculation on the cluster.

Table 2 summarizes the relative computing time along with the model accuracy as discussed. In general, the RNG k-[epsilon] and SST k-[omega] models required roughly two-to-four times the computing time as the indoor zero-equation model, while the v2f-dav model and the LRN model required four-to-eight times as much as the 0-eq. model. The RSM model required the most computing time among all RANS models tested, but the convergence of the RSM model is not as satisfactory as others. Compared to RANS, the computing effort of LES is significantly longer. With the computer clusters, the LES could handle some indoor airflows with simple domains. However, it still can be prohibitively time consuming to use LES for very complicated enclosed environments. The DES required a computing time similar to that of LES, based on the same grid. By carefully designing the grid resolution, the DES could save more computing time from pure LES, although the overall computing time required is still very high.

Table 2. Summary of the Performance of the Turbulence Models
Tested by This Study

Cases        Compared Items       Turbulence Models

                              0-eq.      RNG         SST
                                     k-[epsilon]  k-[omega]

Natural     Mean temperature    B        A            A
convection

            Mean velocity       D        B            A

            Turbulence         n/a       C            C

Forced      Mean velocity       C        A            C
convection

            Turbulence         n/a       B            C

Mixed       Mean temperature    A        A            A
convection

            Mean velocity       A        B            B

            Turbulence         n/a       A            D

Strong      Mean temperature    A        A            A
buoyancy
flow

            Mean velocity       B        A            A

            Turbulence         n/a       C            A

Computing time (unit)           1              2-4

Cases        Compared Items      Turbulence Models

                              LRN-LS  V2f-dav  RSM-IP


Natural     Mean temperature    C       A        A
convection

            Mean velocity       B       A        B

            Turbulence          C       A        C

Forced      Mean velocity       A       A        B
convection

            Turbulence          B       B        B

Mixed       Mean temperature    A       A        B
convection

            Mean velocity       B       A        A

            Turbulence          B       A        A

Strong      Mean temperature    A       A       n/c
buoyancy
flow

            Mean velocity       A       A       n/c

            Turbulence          B       B       n/c

Computing time (unit)              4-8         10-20

Cases        Compared Items    Turbulence Models

                                DES         LES

Natural     Mean temperature     C           A
convection

            Mean velocity        D           B

            Turbulence           C           A

Forced      Mean velocity        C           A
convection

            Turbulence           C           B

Mixed       Mean temperature     B           A
convection

            Mean velocity        B           B

            Turbulence           B           B

Strong      Mean temperature    n/a          B
buoyancy
flow

            Mean velocity       n/a          A

            Turbulence          n/a          B

Computing time (unit)          [10.sup.2]-[10.sup.3]

A = good, B = acceptable, C = marginal, D = poor, n/a = not
applicable, and n/c = not converged.

It is necessary to quantify the model accuracy criteria A, B, C, and D in Table 2. Due to the complexity, a strictly quantified and objective description is difficult. In general, this study used the relative error between prediction and measurement at measured points as a major criterion. If this error is less than 10% or larger than 50% at most measured points, the model accuracy is rated as A or D, respectively. While ratings A and D quantify the extremes, the difference between B and C can be more subtle. Rating B is given to predictions with a relative error less than 20%-30% at most measured points. Rating C is given to the remaining predictions. Note that the relative error calculations for the temperature were based on the nominal temperature differential, which has different definition for each case. In the natural convection case, the nominal temperature difference is the wall temperature difference between the hot and cold walls. For the mixed convection case, it is defined as the difference between the inlet and the outlet temperatures. For the strong buoyancy flow case, it is the difference between the measured local temperature and the environment temperature.

CONCLUSIONS

This study evaluated the overall performance of eight prevalent and/or recently proposed models for simulating airflows in enclosed environments. Four benchmark indoor flow cases were tested that represent the common flow regimes in enclosed environments. In general, the LES provides the most detailed flow features, while the computing time is much higher than the RANS models and the accuracy may not always be the highest. The DES also requires significant computing time in typical indoor flows. For these low-Reynolds-number flows, the DES does not save computing time and the accuracy becomes poorer compared to the LES with the same grid. Further investigation is needed to make conclusive remarks for using DES for airflow simulations in enclosed environments.

Among the RANS models, the v2f-dav and RNG k-[epsilon] models show the best overall performance compared to the other models in terms of accuracy, computing efficiency, and robustness. Both models are recommended for indoor airflow simulations. Although the present investigation has evaluated the selected turbulence models for four common scenarios, other important flow scenarios (e.g., three-dimensional wall jets in ventilated rooms and wind-driven natural-convection flows) should be further investigated.

The SST k-[omega] model did improve the accuracy in a strong buoyant flow scenario without significantly increasing computing time, compared to the RNG k-[epsilon] model. However, the SST k-[omega] model exhibited problems for low turbulence flows. The LRN-LS model and the v2f-dav model require similar computing time, but the LRN-LS model did not perform as well as the other. The RSM model performed reasonably well in the two-dimensional flows but encountered convergence problems in the three-dimensional buoyant flow.

Compared to these advanced turbulence models, the indoor zero-equation model is less accurate. Nevertheless, the model also has its merits. It is simple and always shows good convergence speed. Its results can be used as good initial fields for more advanced models to achieve converged results.

While the turbulence models have different performance in different flow categories, each airflow category favors specific turbulence models. The present study summarizes the best-suited turbulence models for each flow category studied. The v2f-dav and the LES are best suited for the low Rayleigh number natural convection flow in predicting air velocity, temperature, and turbulence quantities. In the forced convection flow with low turbulence levels, the RNG k-[epsilon], LRN-LS, v2f-dav, and LES all performed very well. The v2f-dav, the RNG k-[epsilon], and the indoor zero-equation model show the best accuracy and are suitable for mixed convection flows with low turbulence levels. Although the SST k-[omega] model is less accurate than other models in low Reynolds number flows, it works the best for the high Rayleigh number buoyancy-driven flow. Meanwhile, the LRN-LS and the v2f-dav model are also suitable for the high Rayleigh number flow in the present study.

ACKNOWLEDGMENTS

The authors would like to express their gratitude to Dr. Kazuhide Ito of Tokyo Polytechnic University , who kindly provided the details of his experimental data of the natural convection case. Z. Zhang and Q. Chen would like to acknowledge the financial support of the U.S. Federal Aviation Administration (FAA) Office of Aerospace Medicine through the Air Transportation Center of Excellence for Airliner Cabin Environment Research under cooperative agreement 04-C-ACE-PU. Although the FAA has sponsored this project, it neither endorses nor rejects the findings of this research. The presentation of this information is in the interest of invoking technical-community comment on the results and conclusions of the research.

NOMENCLATURE 

Symbols

C= concentration of scalar variables; constants and coefficients in Table 1

[c.sub.p] = specific heat capacity at constant pressure

D = the LRN modification term in the k-equation of the standard k-[epsilon] model

E = the LRN modification term in the [epsilon]-equation of the standard k-[epsilon] model

F = weighting function in SST k-[omega] model

f = model or the elliptical relaxation in the v2f model

G = the production term in various transport equations

H = enthalpy

h = subgrid scale turbulent heat flux in LES

k = turbulent kinetic energy

Re = Reynolds number

S = equations; magnitude of strain rate tensor 

T = air temperature or turbulence time scale

t = time

u = instantaneous air velocity

[[bar].v'.sup.2] = root mean square of fluctuating velocity normal to wall

x = Cartesian coordinates 

Y = dissipation term in various transport equations

Greek Symbols

[DELTA] = filter width of filter function in LES in Table 3.1

[epsilon] = the dissipation rate of turbulence kinetic energy

[phi] = the symbol of general variable in Equation 1 and Table 1

[TAU] = diffusion coefficient in Equation 1

[kappa] = von Karman constant (0.4187)

[mu] = molecular viscosity of a fluid

[[mu].sub.[tau]] = turbulent viscosity of fluid flow

[rho] = air density

[[sigma].sub.[phi],[tau]] = turbulent Prandtl numbers of variable

[tau] = stress

[OMEGA] = magnitude of rotation tensor rate

[omega] = rate of dissipation per unit turbulent kinetic energy

Superscripts

s = subgrid scale variables in LES

' = fluctuating quantities

Subscripts

i,j = Cartesian coordinates

l,m = Cartesian coordinates

REFERENCES

Betts, P.L., and I.H. Bokhari. 2000. Experiments on turbulent natural convection in an enclosed tall cavity. International Journal of Heat and Fluid Flow 21:675-83.

Blay, D., S. Mergur, and C. Niculae. 1992. Confined turbulent mixed convection in the presence of a horizontal buoyant wall jet. Fundamentals of Mixed Convection, ASME HTD 213:65-72.

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

Davidson, L., P.V. Nielsen, and A. Sveningsson. 2003. Modification of the V2F model for computing the flow in a 3d wall jet. Turbulence Heat and Mass Transfer 4:577-84.

FLUENT. 2005, FLUENT 6.2 Documentation, Fluent Inc., Lebanon, NH.

Germano, M., P. Piomelli, P. Moin, and W.H. Cabot. 1991. A dynamic subgrid-scale eddy viscosity model. Physics of Fluids 3(7):1760-765.

Gibson, M.M., and B.E. Launder. 1978. Ground effects on pressure fluctuations in the atmospheric boundary layer. Journal of Fluid Mechanics 86:491-511.

Ito, K., S. Kato, and S. Murakami. 2000. Model experiment of flow and temperature field in room for validating numerical simulation analysis of newly proposed ventilation effectiveness. Journal of Architecture, Planning and Environmental Engineering 534:49-56. In Japanese.

Launder, B.E., and B.L. Sharma. 1974. Application of the energy dissipation model of turbulence to the calculation of flow near a spinning disk. Letters in Heat Mass Transfer (1):131-38.

Lien, F., and G. Kalitzin. 2001. Computations of transonic flow with the v2f turbulence model. International Journal of Heat and Fluid Flow 22:53-61.

Lilly, D.K. 1992. A proposed modification of the Germano subgrid-scale closure model. Physics of Fluids 4:633-35.

Menter, F.R. 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal 32:1598-1605.

Murakami, S., S. Kato, and R. Yoshie. 1995. Measurement of turbulence statistics in a model fire room by LDV. ASHRAE Transactions 101(2):287-301.

Shur, M., P.R. Spalart, M. Strelets, and A. Travin. 1999. Detached-eddy simulation of an airfoil at high angle of attack. Proceedings of the 4th International Symposium on Engineering Turbulence Modeling and Experiments, Corsica, France, pp. 669-78.

Spalart, P., and S. Allmaras. 1992. A one-equation turbulence model for aerodynamic flows. Technical Report AIAA-92-0439, American Institute of Aeronautics and Astronautics, Reston, VA.

Yakhot, V., and S.A. Orszag. 1986. Renormalization group analysis of turbulence. Journal of Scientific Computing (1):3-51.

Zhai, J.Z., Z. Zhang, W. Zhang, and Q. Chen. 2007. Evaluation of various turbulence models in predicting airflow and turbulence in enclosed environments by CFD: Part 1-Summary of prevalent turbulence models. HVAC&R Research 13(6):853-70.

Zhao Zhang

Student Member ASHRAE

Zhiqiang (John) Zhai, PhD

Member ASHRAE

Wei Zhang, PhD

Member ASHRAE

Qingyan (Yan) Chen, PhD

Fellow ASHRAE

Zhiqiang (John) Zhai is an assistant professor in the Department of Civil, Environmental, and Architectural Engineering, University of Colorado, Boulder. Zhao Zhang is a graduate research assistant, Wei Zhang is an affiliate, and Qingyan (Yan) Chen is a professor of mechanical engineering in the School of Mechanical Engineering, Purdue University, West Lafayette, IN.

Received April 16, 2007; accepted August 29, 2007