- Email: [email protected]
Numerical simulation of turbulent airflow in a ventilated room: Inlet turbulence parameters and solution multiplicity
Accepted Manuscript Title: Numerical Simulation of Turbulent Airflow in a Ventilated Room: Inlet Turbulence Parameters and Solution Multiplicity Author: Erhan Pulat Hıfzı Arda Ersan PII: DOI: Reference:
S0378-7788(15)00094-8 http://dx.doi.org/doi:10.1016/j.enbuild.2015.01.067 ENB 5673
To appear in:
ENB
Received date: Revised date: Accepted date:
12-8-2014 12-1-2015 30-1-2015
Please cite this article as: E. Pulat, H.A. Ersan, Numerical Simulation of Turbulent Airflow in a Ventilated Room: Inlet Turbulence Parameters and Solution Multiplicity, Energy and Buildings (2015), http://dx.doi.org/10.1016/j.enbuild.2015.01.067 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Numerical Simulation of Turbulent Airflow in a Ventilated Room: Inlet Turbulence Parameters and Solution Multiplicity Uludag University, Faculty of Engineering, Mechanical Engineering Department, Gorukle Campus, TR-16059 Bursa, Turkey 2 Uludag University, School of Natural and Applied Sciences, Mechanical Engineering Department, Gorukle Campus, TR-16059 Bursa, Turkey * Corresponding Author: Tel.: +90 224 4428176, Fax: +90 224 4428021, E-Mail: [email protected]
ip t
1
Erhan PULAT1,* and Hıfzı Arda ERSAN2
cr
Abstract
ed
M
an
us
In this study, airflow and temperature distributions in the well-known International Energy Agency (IEA) Annex 20 room are predicted numerically to investigate the effects of the inlet turbulence intensity and the length scale on the flow characteristics, while considering the possibility of solution multiplicity related to chaos theory. The flow is considered to be turbulent, steady, incompressible, and two-dimensional. Computations are performed using the standard k-ε, RNG k-ε, standard k-ω, and SST k-ω turbulence models with scalable and automatic wall functions, and the results are compared with numerical and experimental results from the literature. The validated turbulence model is then used to investigate the effects of the turbulence intensity and the length scale. At a low inlet turbulence intensity value (Tu=0.01), the length scale variation has no influence on the flow pattern. However, the length scale affects the flow pattern at a high inlet turbulence intensity value (Tu=0.4). At constant low and medium length scale values, an increase in the inlet turbulence intensity from 0.01 to 0.4 affects the flow pattern. However, an increase in the turbulence intensity has no influence on the flow pattern at a constant high length scale value.
Ac ce
1. Introduction
pt
Keywords: Room Ventilation, Eddy Viscosity Turbulence Models, Turbulence Intensity and Length Scale, Solution Multiplicity, CFD
On average, people spend 90% of their life indoors [1]. Good airflow is important in ventilated spaces to provide good air quality and thermal comfort for the occupants, to maintain specified thermal and flow conditions for industrial processes and to ensure the efficient use of energy. The importance of ventilation related energy use is increasing and may represent up to 50% of the total energy use of a building, particularly for certain typologies such as office buildings [1]. Computational fluid dynamics (CFD) is one of the primary methods used to assess indoor airflow [2, 3, 4], among other methods such as the Building Energy Simulation (BES) approach and zonal models. In recent years, CFD has taken a prominent role in the simulation of indoor environment airflow problems [5, 6]. However, indoor airflow is a complex system, and new modeling approaches such as time series modeling can be found in the literature [7]. Pioneering work using the CFD technique to simulate room air motion was performed by Nielsen [8]. The flow inside a building environment is generally turbulent [2], and the IEA Annex 20 room, described in Nielsen [9], is one of the experimental setups for studying the suitability of room turbulent air distribution models [10, 11]. Aside from these studies, which were mainly focused on the comparison of turbulence models, and although the effects of turbulent boundary conditions in ventilated enclosures were identified in Reynolds-averaged 1 Page 1 of 23
Ac ce
pt
ed
M
an
us
cr
ip t
Navier-Stokes (RANS) simulations by Cao and Meyers [12], turbulent conditions have been largely neglected in simulations in the indoor ventilation literature [13, 14]. As observed by Abdilganie et al. [15] and Cao and Meyers [12], further investigation is required. Solution multiplicity refers to the existence of more than one solution for the overall flow patterns in a room under the same boundary conditions. Different initial conditions can lead to different solutions [16]. This phenomenon is in essence related to the chaotic behavior in nonlinear dynamic systems that can be encountered in all fields of science, including engineering sciences [17]. Zheng et al. [18] studied the nonlinear equations for ventilation flows in an enclosure with a single opening. Theoretical and computational nonlinear dynamic analysis of natural ventilation in a two-zone building was performed by Yang et al. [19, 20]. Van Schijndel [21] performed numerical research of the chaotic behavior of the airflow in an ordinary ventilated room. Chen [22] and Lemaire [23] reported two different airflow pattern results from simulations of the IEA Annex 20 room. The numerical study of Xin et al. [24] was intended to reveal some of the possible sequences of bifurcations leading to chaotic mixed convection in confined spaces. Recently, Dillon et al. [25] investigated chaotic natural convection in a tall rectangular cavity with non-isothermal walls. In addition to the above studies, Yee et al. [26] paid particular attention to isolating the different nonlinear behaviors and spurious dynamics caused by some of the numerical uncertainties observed in practical CFD computations. The multiplicity of airflow in naturally ventilated buildings [27, 28], smoke flows [29], and mixed convection in 3D ventilated cavities [30] has been confirmed by experiments (Müllejans, 1966 in [16], [30, 31, 32]). As seen in the above CFD and experimental studies, different initial conditions in CFD studies can lead to different solutions in both steady and transient calculations, and one stable solution can shift to another when there is sufficient perturbation. For example, in natural ventilation, the airflow pattern and the solution multiplicity are dependent on air change parameters [27, 28, 31]. In the case of a ventilated room, chaotic computational behavior is observed when very small changes are made in the supply air temperature [21, 32]. In mixed convection in confined spaces, multiplicity is dependent on two main system parameters, the Re and Ar numbers [22, 23, 24]. Although turbulent inlet parameters were studied from the point of view of flow and heat transfer in some of the studies mentioned above [12, 13 14, 15], solution multiplicity was not considered. In addition to providing a comparison of the popular RANS-based eddy viscosity turbulence models incorporating the scalable wall function and automatic near-wall treatment approaches, the main objective of this study is to investigate the effects of the inlet turbulence intensity and the length scale on the flow characteristics of the IEA Annex 20 room while considering the possibility of solution multiplicity. After a brief explanation of turbulence models and near-wall treatment methods is provided, a validated RNG k-ε turbulence model with a scalable wall function is used for this objective. 2. Materials and methods 2.1. Geometry The geometry of the IEA Annex 20 room described by Nielsen [9] is shown in Figure 1. The dimensions of the two-dimensional room are H=3 m, L=9 m, h=0.168 m and t=0.48 m, and it is equipped with an inlet and outlet to avoid numerical problems, as in Voigt’s study [10]. In a recent study, Olmedo and Nielsen [33] checked the two dimensionality of the test case, and another aim of their study is to analyze the results obtained considering the case as steady or unsteady. It has been concluded that the IEA Annex20 case solved with time
2 Page 2 of 23
dependent equations and k-epsilon model gives a solution similar to the solution obtained with steady state three dimensional and the two dimensional equations. 2.2. Turbulence models and near-wall treatment methods
Ac ce
pt
ed
M
an
us
cr
ip t
Computations were performed using the standard k-ε [34], RNG k-ε [35], standard k-ω [36] and SST k-ω [37] turbulence models. These two-equation models with eddy viscosity are widely used, as they offer a good compromise between numerical effort and computational accuracy. The advantage of the standard k-ε model from a numerical point of view is that it is robust and reliable. Although the physics is treated in a simplistic manner, the model works surprisingly well in many types of flows [38]. In a recent study, Klinzing and Sparrow [39] compared five turbulence models commonly offered in commercially available CFD software, including FLUENT, ANSYS CFX, and STAR-CD. They found that the most suitable turbulence model for external flows was the standard k-ε model with enhanced nearwall treatment. Chen [40] previously evaluated the performance of five k-ε models for the simulation of forced-convection, natural-convection, and mixed-convection room airflow. The results from the RNG k-ε model were slightly better than the standard k-ε model, and the other models, including the low Reynolds number k-ε, the two-layer k-ε, and the two-scale kε models, were not stable. Therefore, Chen recommended the RNG k-ε model for indoor airflow simulations. The k-ω model was compared with the standard and low Reynolds number k-ε models to model recirculating flows in a ventilated space by Peng et al. [41]. To improve the prediction accuracy of the standard k-ω model, they proposed some modifications. The SST k-ω model is a modification of the k-ω model, but it predicts a recirculation zone larger than that measured in PIV experiments in a water scale model of the Annex 20 room [10]. The k-ε and k-ω models seek to solve a modified set of transport equations by introducing mean and fluctuating components ( U i U i u i ). The continuity equation is not altered, and the momentum and scalar transport equations contain turbulent flux terms in addition to the molecular diffusion fluxes. These are the Reynolds stresses, u i u j . These terms arise from the nonlinear convective term in the unaveraged equations. The resulting steady-state RANS equations, assuming air is Newtonian and incompressible, and dropping the bar for averaged quantities, are as follows: Continuity: U j 0 x j Momentum: U iU x j
j
(1)
p ij u i u xi x j
j
S
(2)
M
Energy: T ( U j h tot ) ( u j h ) U i ( ij u i u j ) S E x j x j x j
(3)
where is the molecular stress tensor. S M is a source term added for buoyancy calculations:
3 Page 3 of 23
S M ,buoy ( ref ) g
(4)
The density difference is evaluated using the Boussinesq model because the density variation is driven only by small temperature variations in the room ventilation. The buoyancy source term is approximated as
ref ref (T Tref )
(5)
1 T
(6)
p
cr
ip t
where is the thermal expansivity
us
The mean total enthalpy is given by
an
1 htot h U iU i k 2
(7)
The total enthalpy contains a contribution from the turbulent kinetic energy, k, given by 1 2 ui 2
M
k
(8)
ed
Eddy viscosity models use the gradient diffusion hypothesis to relate the Reynolds stresses to the mean velocity gradients and the eddy (turbulent) viscosity (the Boussinesq approach) [41, 42]. 2 ij k t U k 3 x k
pt
U i U j ui u j t x xi j
(9)
Ac ce
where t is the eddy viscosity that is linked to the turbulent kinetic energy and its dissipation rate via the relation
t C
k2
(10)
The main advantage of this approach from the observation that the representation of Reynolds stress term is of exactly the same form as that of the diffusion terms in the original equations. They can be combined if an effective viscosity is defined as the sum of the laminar viscosity and the turbulent viscosity,
eff t
(11)
In two-equation models, the turbulence velocity scale is computed from the turbulent kinetic energy, which is obtained from the solution of its transport equation. The turbulence length scale is estimated from two properties of the turbulence field, usually the turbulent kinetic 4 Page 4 of 23
energy and its dissipation rate. The dissipation rate of the turbulent kinetic energy is obtained from its transport equation. Standard versions of four models used have been explained briefly as follows. Standard k-ε model
ip t
The k-ε model introduces two new variables into the system of equations. In this model, C is a constant, 0.09. The values of k and ε come directly from the differential transport equations for the turbulent kinetic energy and the turbulence dissipation rate [34].
cr
RNG k-ε model
us
The RNG k-ε model is based on a renormalization group analysis of the Navier-Stokes equations. The transport equations for turbulence generation and dissipation are the same as those for the standard k-ε model, but the model constants differ, with the constant C 1 replaced by the function C 1RNG [35].
an
Standard k-ω model
k
(12)
ed
t
M
One of the advantages of the k-ω formulation is the near-wall treatment for low Reynolds number computations. The k-ω models assume that the turbulence viscosity is linked to the turbulent kinetic energy and turbulent frequency via the relation
pt
The k-ω model solves two transport equations, one for the turbulent kinetic energy, k, and one for the turbulent frequency, ω [36]. The Shear Stress Transport (SST) k-ω model
Ac ce
The SST k-ω model accounts for the transport of the turbulent shear stress and gives highly accurate predictions of the onset and amount of flow separation under adverse pressure gradients [37]. The standard k-ω model does not account for the transport of the turbulent shear stress. This results in an overprediction of the eddy viscosity. The proper transport behavior can be obtained by incorporating a limiter in the formulation of the eddy viscosity as follows:
t a1 k max(a1 , SF2 )
(13)
F2 is a blending function that restricts the limiter to the wall boundary layer, because the underlying assumptions are not correct for free shear flows. S is an invariant measure of the strain rate. The blending function is critical to the success of the method. Its formulation is based on the distance to the nearest surface and on the flow variables:
5 Page 5 of 23
F2 tanh arg 22
with
2 k 500 arg 2 max ' , 2 y y
(14)
Near-wall treatment methods
Ac ce
pt
ed
M
an
us
cr
ip t
When computing wall-bounded turbulent flows it is necessary to alter the turbulence model in the vicinity of the wall. Two basic approaches are generally used: the so-called wall function approach and the low Reynolds number model [43]. Near-wall treatment is crucial in any turbulence model, and the most popular approach is the wall function approach, due to its simplicity. Basically, the method requires that the first grid point near the wall be located in the fully turbulent region. In the low Reynolds number approach, the basic idea is to integrate the equations up to the wall. This approach requires the clustering of grid points near the wall, which increases the total number of grid points required for an accurate solution. In flows with steep pressure gradients or swirl, the wall function representation of the near-wall region is not sufficiently universal to allow accurate predictions [43]. However, according to Chow and Gao [44], the standard k-ε model provides steady numerical results without fluctuations and might be a suitable method for representing the separated and reattached flows if the wall functions are applied properly. Additionally, Mohammadi and Medic [45] have shown that the careful application of numerical methods and adequate numerical dissipation lead to a substantial improvement in the results of the classical k-ε model. A useful discussion of wall function treatments, including the mixed convection effect, is provided by Craft et al. [46]. An extension of this method is used in ANSYS CFX as follows. The main problem in the wall function approach is the danger of excessive mesh refinement, so that the physical assumptions are especially problematic in flows at lower Reynolds numbers in the near-wall mesh. To avoid this problem, the scalable wall function (SWF) formulation in ANSYS CFX is used in the k-ε model. The basic idea behind this approach is to limit the y+ value used in the logarithmic formulation by a lower bound of y͂+ = max(y+, 11.06). The computed y͂+ is not allowed to fall below this limit. However, in flows at low Reynolds numbers, this can cause an error in the displacement thickness, and it is therefore desirable to offer a formulation that will automatically switch from a wall function to a low Reynolds number near-wall formulation as the mesh is refined. The main idea behind this formulation, called the automatic near-wall treatment (ANWT), is to blend the wall value for ω between the logarithmic and the near-wall formulation. This method has been used in k-ω models [47]. A blending function that depends on y+ can therefore be defined that can be used to derive a formulation describing the behavior of ω in linear and logarithmic near-wall regions as follows [42, 47, 48]: 2 1 y vis y log2 y
(15)
A similar formulation is used to describe the relation between the velocity at the first grid point and the wall shear stress:
u
u 4 uVis
4
log 4
(16)
In both the wall function and the automatic near-wall treatment approaches, in the treatment of the energy equation near the wall, an algebraic formulation is required to link the temperature and the heat flux. The formulation of Kader [49] is used in ANSYS-CFX: 6 Page 6 of 23
Pr y e 2.12 ln 1 y Pr e
1
(17)
2
with
(18) (19)
ip t
1 where (Pr) 3.85 Pr 3 1.3 2.12 ln(Pr) 4 0.01(Pr y ) , 1 5 Pr 3 y
and the non-dimensional temperature is defined as
2.3. Boundary conditions and computational procedures
cr
qw Tw T , T c p u T
(20)
us
ghT0 u0
pt
Ar
ed
M
an
The modeled inlet conditions were a uniform velocity of 0.455 m/s, an inlet turbulence intensity of 4 % and a temperature of 20 ºC, consistent with Nielsen’s test case [9]. The turbulence length scale used was LS=0.0156, which is very close to Nielsen’s test case value of 0.0168. The outlet condition was a constant pressure of zero. The velocity was set to zero on the room walls and the inlet and outlet walls to impose a no-slip condition. Based on the conditions in the supply opening, the inlet Reynolds number was 5000. Air was treated as an ideal gas, and its density and viscosity were assumed to be 1.116 kg/m3 and 1.725 Pa·s [50]. The left and right wall and the ceiling were assumed to be adiabatic. A constant heat flux of 63.08 W/m2 was added along the floor in the non-isothermal case, similar to Lemaire’s study [23]. The effect of buoyancy is commonly expressed in the ventilation literature by the Archimedes number, Ar: (21)
Ac ce
where ΔT0 is the temperature difference between the inlet and outlet, and h and u0 are the height and velocity of the inlet, respectively. Archimedes number was 0.173 in the present study. The characteristics of turbulent flow are also commonly expressed by the turbulence intensity, the turbulence length scales, and the turbulent kinetic energy and its dissipation rate. The turbulence intensity, Tu, is the standard deviation divided by the mean velocity ( u ): Tu
u' x100% u
(22)
where the velocity fluctuations, u ' , are in the main direction of the flow. The turbulence length scales are the integral scale and the microscale. It is assumed that the turbulent motion consists of the superposition of eddies of various sizes. The integral scale identifies the average size of the largest eddies, whereas the microscale is a measure of the smallest eddies, which are mainly responsible for dissipation [51]. The ANSYS CFX R13 code uses the finite volume method, and pressure and velocity are coupled using the SIMPLEC algorithm with Rie-Chow interpolation. Buoyancy is modeled 7 Page 7 of 23
using the Boussinesq approach. A mesh independency study was performed by comparing three different mapped mesh structures (coarse, medium, and fine), corresponding to 12168, 19506 and 29054 nodes, using the standard k-ε model. The mesh structure with 29054 nodes was chosen by comparing the u/uo velocity profiles for the different numbers of nodes, as seen in Figure 2. Further information about mesh structures can be found in [52]. 3. Results and discussion
Ac ce
pt
ed
M
an
us
cr
ip t
In Figure 3, the velocity vectors and isotherms obtained under the specified boundary conditions from the popular RANS-based eddy viscosity turbulence models using the relatively new near-wall modeling approaches, namely the scalable wall function (SWF) and the automatic near-wall treatment (ANWT), are compared with the results of Lemaire’s [23] numerical study, because that study includes the non-isothermal case for the IEA Annex 20 room. As shown in Figure 3, the flow patterns obtained from the standard and RNG k-ε models are consistent with Lemaire’s flow pattern, that is, the main vortex is in the clockwise direction. In the case of the standard k-ω and SST k-ω models, however, the flow pattern changes completely: the cold jet falls down, causing a reversal of the main vortex. Due to this different flow pattern, the isotherms in the k-ω models are also not consistent with Lemaire’s isotherms. Additionally, the isotherms are predicted more accurately by the RNG k-ε model than the standard k-ε model when compared with Lemaire’s results, as shown in Figure 3a, 3b, and 3c. The streamlines in the RNG k-ε model are also qualitatively consistent with Nielsen’s [8] experimental results, as seen in Figure 4. For this reason, the effects of the inlet turbulence intensity and the length scale on the flow pattern and temperature distribution have been investigated with the RNG k-ε model. Additionally, as mentioned previously, this model was recommended for indoor airflow simulations by Chen [40]. Although one model might perform well for one flow but poorly for another, the performance of the RNG k-ε model was relatively stable [53]. Yousaf et al. [54] used ANSYS CFX to apply the SST k-ω and RNG kε models, using the same near-wall modeling approach as the present study, to more complex indoor airflows dominated by buoyancy effects generated by human occupancy and equipment. They reported differences of 6% and 8% in the temperature values near the inlet in the SST and RNG models, respectively; at the other sensor points, the RNG model predictions differed by less than 2.6%. The SST model appears to overestimate the temperature field. To investigate the effects of the turbulence intensity (Tu) and length scale (LS) on the flow pattern and temperature distribution, two Tu values (Tu = 0.01 and 0.4) and three LS values (LS = 0.00156, 0.0156 and 0.156) are considered, as shown in Figure 5. The seven cases shown in Table 1 are investigated using the RNG k-ε model. Case 4 is the comparison condition, representing the test case boundary conditions described in section 2.3. If the turbulence intensity is greater than 0.05, it is commonly assumed that the flow is in the high free-stream turbulence condition [55]. Although Tu value of 0.4 can be treated as very high for ventilation flows, high Tu values are encountered in air supply devices such as diffusers. For example, Cehlin and Moshfegh [56] was reported that turbulence intensity in the openings was around 0.1-0.3% while the intensity outside the diffuser in the near zone was up to 0.6. According to Cao et al. [57], in a ventilated room, due to wall effects, the magnitude of turbulence intensity can grow to a value greater than 0.2 from 0.02. In addition, the levels of inlet turbulence intensities investigated by Cao and Meyers [12] were 0.05 and 0.3, by Loomans [58] were 0.1 and 0.35, and by Joubert et al. [13] were 0.04 and 0.374. So high Tu values can be treated as reasonable, and as mentioned by van Hooff et al. [3], the flow pattern encountered in room ventilation can be challenging due to the presence of transitional flow, turbulence anisotropy, and adverse pressure gradient. Really, as mentioned by Cehlin and 8 Page 8 of 23
Ac ce
pt
ed
M
an
us
cr
ip t
Moshfegh [56] prediction of room airflows is a challenging task because of the complex nature of turbulence. Increase of turbulent intensity at the inlet also leads to an increased turbulent intensity in the rest of domain as shown by Cao and Meyers [12]. As known turbulence increases heat transfer rate and high Tu values may cause high heat transfer coefficients in ventilation studies together with room heating. From built environment point of view, understanding airborne contaminant transport and distribution in buildings is essential for addressing issues such as home and office air quality, worker exposure in manufacturing workrooms, product protection in clean rooms, accident response planning in nuclear facilities, avoiding fatalities in confined spaces and infection control in medical settings [59], and as mentioned by Cao and Meyers [12], when pollutant dispersion is also considered, the inlet turbulent quantities both turbulent intensity and length scales have an important effect on pollutant dispersion. As shown in Figure 5 (I), there is no influence of the length scale on the flow pattern or the temperature distribution at the low inlet turbulence intensity value (Tu=0.01). Although the length scale varies by a factor of a hundred (from 0.00156 to 0.156), the flow pattern and temperature distribution are not affected by this increase. However, at the high turbulence intensity value (Tu=0.4), at the low and medium length scale values (LS=0.00156 and 0.0156), the flow pattern changes completely, producing flow patterns similar to those obtained by the standard and SST k-ω models for the test case, as shown in Figure 3d and e. The cold jet falls down, causing a reversal of the main vortex. At the high length scale value (LS=0.156), the clockwise flow pattern is obtained again, and due to this change in the flow pattern, the temperature distribution is also clockwise. A similar flow pattern change was reported by Chen [22] and Lemaire [23] as a result of a small change in the Archimedes number (Ar). Their transition Archimedes numbers are different due to their different inlet sizes and Reynolds numbers [22]. This phenomenon was also reported by Heiselberg et al. [31] and Ezzouhri et al. [30] for different room dimensions, and was ascribed to a hysteresis effect. Ezzouhri et al. [30] concluded that the numerical bifurcation between the two flow structures was highly dependent on the turbulence level at the inlet. In the present study, at fixed small and medium length scale values, as the inlet turbulence intensity is increased from 0.01 to 0.4, the clockwise flow pattern is converted into a counterclockwise flow. At the high fixed inlet turbulence intensity value, if the length scale is decreased from 0.156 to 0.0156 or 0.00156, the clockwise flow pattern is also converted into a counterclockwise flow, as shown in Figure 5. These results are consistent with the findings of Hancock and Bradshaw [60], who showed for the first time that the free stream turbulence length scale was also an important parameter, in addition to the free stream turbulence intensity [55], and showed that the free stream turbulence effect decreased as the length scale increased. Really, if the length scale is high (LS=0.156) as seen from Figure 5c, flow pattern is not affected from increase in the free stream turbulence intensity from 0.01 to 0.4. A comparison of the present results with Lemaire’s [23] results for the same geometry using the standard k-ε model and Ezzouhri et al.’s [30] results for a different geometry using the LES approach is given in Figure 6. Additionally, Ezzouhri et al. [30] observed that the dynamic LES model not only correctly predicts the mean flow characteristics of the flow but is also able to correctly reproduce the flow bifurcation and the hysteresis effect. However, in the present study, as well as in Chen’s study [22] using the Lam-Bremhorst low Reynolds number k-ε model, Lemaire’s study [23] using the standard k-ε model, and Nielsen et al.’s study [61] using the standard k-ε model, it has been shown that RANS simulations have capability comparable to LES simulations. Similarly, in the recent study by Durrani et al. [62], the RNG k-ε model was used in addition to the LES model to represent turbulence in an unsteady RANS (URANS) simulation. They concluded that the URANS method was unable to capture the details of the flow structures
9 Page 9 of 23
pt
4. Conclusions
ed
M
an
us
cr
ip t
predicted by LES, but by its nature, the LES approach required approximately 500% more computation time than URANS. A better understanding of solution multiplicity in ventilation flows is very important because researchers and practitioners may encounter this phenomenon. For example, in the study of O’Halloran et al. [63], two turbulence models, the RNG k-ε model and the Reynolds stress model (RSM), were evaluated in simulations and compared with experimental particle image velocimeter (PIV) data. In their study, although qualitatively the RNG k-ε model and RSM results compared well with each other and with the PIV results, one difference was apparent in the RNG k-ε model result. They expressed this difference as follows: “In the recirculation area near the inlet of the jet, from approximately x=0 to 150 mm, the direction of flow is opposite to that of the PIV results and the RSM model results. The flow reversal appears to result from a horseshoe vortex that wraps around the inlet flow. The numerical reason for this difference could not be determined, neither by the researchers nor by the software manufacturer, and further investigation is needed to understand the cause of the differences.” This difference might be related to solution multiplicity. As shown in Figure 3, different turbulence models produce different solutions under the same boundary conditions. Additionally, as shown in Figure 5, a single turbulence model can also produce different solutions under different boundary conditions. However, as reported by Nielsen et al. [61] and Li and Nielsen [16], the existence of two solutions under the same boundary conditions is not only found experimentally, as in the results of Müllejans for the same Ar and Re numbers [Müllejans, 1966 in [16]]; two solutions can be obtained numerically from a single turbulence model for the same Ar and Re numbers [54, 61]. Consequently, in any validation process in the evaluation of a suitable turbulence model, one must consider whether the flow is expected to have multiple solutions. For example, in this study, the flow pattern predictions from the standard and SST k-ω models producing counterclockwise main vortices represent the second solution in each case, and the predictions of these models must be compared using the measurements from this second solution, if available, before making a final decision about their suitability.
Ac ce
Airflow and temperature distributions in the well-known IEA Annex 20 room have been predicted numerically to investigate the effects of the inlet turbulence intensity and the length scale on the flow characteristics while considering solution multiplicity. For this purpose, the RNG k-ε model is preferred because of its slightly better prediction of the temperature distribution, although similar flow patterns are obtained from the standard k-ε model. The inlet turbulence parameters have a great influence on the flow pattern, from the viewpoint of solution multiplicity, at similar Reynolds and Archimedes numbers, and the free stream turbulence effect decreases as the length scale increases. But to determine a threshold value of Tu that the flow pattern changes from cw direction to ccw direction, it is required consideration of a few additional Tu values. Such computations is planned as a future study. In addition, it has been shown that present steady RANS simulations have capability comparable to LES simulations required more computation time than Unsteady RANS in literature. Consequently, in any indoor airflow process or validation study, one must consider whether the flow is expected to have multiple solutions.
Nomenclature β
thermal expansivity [1/K]
cp
constant pressure specific heat [J/kgK] 10 Page 10 of 23
dissipation rate of turbulent kinetic energy [m2/s3]
g
gravitational acceleration [m/s2]
h
enthalpy [J/kg]
k
turbulent kinetic energy [m2/s2]
λ
thermal conductivity [W/mK]
μ
viscosity [kg/ms]
ν
kinematic viscosity [m2/s]
p
pressure [Pa]
qw
wall heat flux [W/m2]
ρ
density [kg/m3]
S
source term
T
temperature [K]
ω
turbulent frequency (ε / k) [1/s]
uτ
dimensionless friction velocity
U
velocity [m/s]
y+
dimensionless distance from the wall
cr us an M
ed
pt
References
ip t
ε
Working Group 16, Ventilation, good indoor air quality and rational use of energy, European Commission Joint Research Centre, Institute For Health and Consumer Protection, Physical and Chemical Exposure Unit, Report No 23, EUR 20741 EN, 2003.
[2]
Q. Chen, Computational Fluid Dynamics for HVAC: Successes and failures, ASHRAE Transactions 103-1 (1997) 178-187. T. van Hoof, B. Blocken, G.J.F. van Heijst, On the suitability of steady RANS CFD for forced mixing ventilation at transitional slot Reynolds numbers, Indoor Air 23 (2013) 236-249. E. Djunaedy, J.L.M. Hensen, M.G.L.C. Loomans, Development of a guideline for selecting a simulation tool for airflow prediction, in: Proceedings of the 8th IBPSA conference BuildingSimulation’2003, Eindhoven, 2003, pp. 267-274. R.A. Pitarma, J.E. Ramos, M.E. Ferreira, M.G. Carvalho, Computational fluid dynamics an advanced active tool in environmental management and education, Manag. of Environ. Qual.: An Int. J. 15 (2004) 102-110. P.V. Nielsen, Computational fluid dynamics and room air movement, Indoor Air 14-Suppl 7 (2004) 134-143. M. Maciejewska, A. Szczurek, G. Sikora, A.Wylomanska, Diffusive and subdiffusive dynamics of indoor microclimate: A time series modeling, Phys. Rev. E 86 (2012) 031128-1-11. P.V. Nielsen, Flow in air conditioned rooms, PhD thesis from the Technical University of Denmark, 1974.
[3] [4] [5] [6] [7] [8]
Ac ce
[1]
11 Page 11 of 23
[15] [16] [17] [18] [19] [20] [21] [22] [23]
[24] [25] [26] [27] [28] [29] [30]
ip t
cr
us
[14]
an
[13]
M
[12]
ed
[11]
pt
[10]
P.V. Nielsen, Specification of a two-dimensional test case – International Energy Agency, energy conservation in buildings and community systems, Annex20: Air flow pattern within buildings, Aalborg University Technical Report, 1990. L.K. Voigt, J.N. Sorensen, J.M. Pedersen, M. Brons, Review of four turbulence models using topology, in: Proceedings of the 8th IBPSA conference BuildingSimulation’2003, Eindhoven, 2003, pp. 1325-1331. M. Omri, N. Galanis, Evaluation of confined natural and forced convection predictions by different turbulence models, Int. J. of Numer. Methods for Heat and Fluid Flow 19 (2009) 5-24. S-J. Cao, J. Meyers, Influence of turbulent boundary conditions on RANS simulations of pollutant dispersion in mechanically ventilated enclosures with transitional slot Reynolds number, Build. and Environ. 59 (2013) 397-407. P. Joubert, A. Sandu, C. Beghein, F. Allard, Numerical study of the influence of inlet boundary conditions on the air movement in a ventilated enclosure, in: Proceedings of the 5th International Conference on Air Distribution in Rooms Roomvent’96, Yokohama,1996, pp. 235-242. S.A. Al-Sanea, M.F. Zedan, M.B. Al-Harbi, Effect of supply Reynolds number and room aspect ratio on flow and ceiling heat-transfer coefficient for mixing ventilation, Int. J. of Therm. Sciences 54 (2012) 176-187. A.M. Abdilghanie, L.R. Collins, D.A. Caughey, Comparison of turbulence modeling strategies for indoor flows, ASME J. of Fluids Engineering 131(2009), 051402-1-18. Y. Li, P.V. Nielsen, CFD and ventilation research, Indoor Air 21(2011), 442-453. G. Baker, F.A. McRobie, J.M.T. Thompson, Implications of chaos theory for engineering science, IMechE J. of Mechanical Engineering Sci Part C 211 (1997) 349-363. Z. Zheng, P. Xu, Y. Li, G. Chen, Nonlinear resonance and quasi-periodic solutions for ventilation flows in a single opening enclosure, Int. J. Bifurcation and Chaos 15 (2005) 18011808. L. Yang, P. Xu, Y. Li, Nonlinear Dynamic Analysis of natural ventilation in a two-zone building: Part A-Theoretical analysis, ASHRAE HVAC&R Research 12 (2006) 231-255. L. Yang, P. Xu, Y. Li, Nonlinear Dynamic Analysis of natural ventilation in a two-zone building: Part B-CFD simulations, ASHRAE HVAC&R Research 12 (2006) 257-278. A.W.M. van Schijndel, Chaotic behavior of the airflow in a ventilated room, in: Proceedings of the COMSOL Conference, Milan,2009, pp. -.(http://www.comsol.nl/paper/chaotic-behavior-ofthe-airflow-in-a-ventilated-room-6885 (Accessed in Nov. 18, 2013) Q. Chen, Int. Energy Agency, Energy conservation in building and community systems, Annex 20: Air flow patterns within buildings, Subtask 1: Room Air and Contaminant Flow, Research Item 1.46 Simulation of simple test cases, Report No AN20.1-CH-90-ETHZ15, 4 March 1991. A.D. Lemaire, Int. Energy Agency, Energy conservation in building and community systems, Annex 20: Air flow patterns within buildings, Subtask 1: Room Air and Contaminant Flow, Research Item No 1.46NL Simulation of simple test cases 2D1,2D2, Report No AN20.1-NL-90TNO-TPD10, May 1991. G. Xin, W.K. Chow, S. Liu, Windows and routes to chaos in mixed convection in confined spaces, Chaos Solitons Fractals 15 (2003) 543-558. H. Dillon, A. Emery, A. Mescher, Analysis of chaotic natural convection in a tall rectangular cavity with non-isothermal walls, Frontiers in Heat Mass Transfer 4 (2013) 1-9. H.C. Yee, J.R. Torczynski, S.A. Morton, M.R. Visbal, P.K. Sweby, On spurious behavior of CFD simulations, Int. J. For Num. Methods in Fluids 30 (1999) 675-711. Y. Li, A. Delsante, Natural ventilation induced by combined wind and thermal forces, Build. and Environ. 36 (2001) 59-71. Y. Li, A. Delsante, Z. Chen, M. Sandberg, A. Andersen, M. Bjerre, P. Heiselberg, Some examples of solution multiplicity in natural ventilation, Build. and Environ. 36 (2001) 851-858. J. Gong, Y. Li, Solution Multiplicity of smoke flows in a simple building, in: Fire Safety Science Proceedings of the 9th International IAFSS Symposium, Karlsruhe, 2008, pp. 895-906. R. Ezzouhri, P. Jaubert, F. Penot, S. Mergui, Large eddy simulation of turbulent mixed convection in a 3D ventilated cavity: Comparison with existing data, Int. J. of Therm. Sciences 48 (2009) 2017-2024.
Ac ce
[9]
12 Page 12 of 23
Ac ce
pt
ed
M
an
us
cr
ip t
[31] P. Heiselberg, Y. Li, A. Andersen, M. Bjerre, Z. Chen, Experimental and CFD evidence of multiple solutions in a naturally ventilated building, Indoor Air 14 (2004) 43-54. [32] J. van Schijndel, Simulation and experimental validation of chaotic behavior of airflow in a ventilated room, in: Proceedings of the 9th Nordic Symposium on Building Physics NSB 2011, Tampere,2011, vol.2, pp. 235-242.(http://alexandria.tue.nl/openaccess/Metis250706.pdf) (Accessed in Nov. 18, 2013) [33] I. Olmedo, P.V. Nielsen, Analysis of the IEA 2D test. 2D, 3D, steady or unsteady airflow? Aalborg University Department of Civil Engineering DCE Technical Report No.106, 2010. [34] B.E. Launder, D.B. Spalding, Lectures in Mathematical Models of Turbulence, Academic Press, London, 1972. [35] V. Yakhot, S.A. Orszag, Renormalization group analysis of turbulence I. Basic theory, J. Sci. Comput. 1 (1986) 1-51. [36] D.C. Wilcox, Reassessment of the scale determining equation for advance turbulence models, AIAA J. 26 (1988) 1299-1310. [37] D.C. Wilcox, Turbulence Modeling for CFD, DCW Industries Inc,La Canada, CA, 1994. [38] L. Davidson, P.V. Nielsen, Large Eddy Simulations of the flow in a three dimensional ventilated room, in: Proceedings of the 5th International Conference on Air Distribution in Rooms Roomvent’96, Yokohama,1996, pp. 161-2168. [39] W.P. Klinzing, E.M. Sparrow, Evaluation of turbulence models for external flows, Numerical Heat Transfer Part A Applications 55 (2009) 205-228. [40] Q. Chen, Comparison of different k-ε models for indoor airflow computation, Numerical Heat Transfer Part B Fundamentals 28 (1995) 353-369. [41] S-H. Peng, L. Davidson, S. Holmberg, The two equation turbulence k-ω model applied to recirculating ventilation flows, Thermo and Fluid Dynamics, Chalmers University of Technology, Göteborg: Report 96/13. [42] ANSYS CFX-Solver Theory Guide, 2009. [43] W.C. Lasher, D.B. Taulbee, The effect of the near-wall treatment on the behavior of ReynoldsStress turbulence models for reattaching flows, in: Seventh Symposium on Turbulent Shear Flows, Stanford University, U.S.A., August 21-23, 1989, pp. 10.5.1-10.5.6. [44] W.K. Chow, Y. Gao, On using the wall function and the staggered grid arrangement at upwind sharp convex corners in simulating approach airflow, Int. J. Comput. Fluid Dyn. 19 (2005) 381398. [45] B. Mohammadi, G. Medic, A critical evaluation of the classical k-ε model and wall-laws for unsteady flows over bluff bodies, Int. J. Comput. Fluid Dyn. 10 (1998) 1-11. [46] T.J. Craft, S.E. Gant, A.V. Gerasimov, H. Iacovides, B.E. Launder, Development and application of wall-function treatments for turbulent forced and mixed convection flows, Fluid Dyn. Res. 38 (2006) 127-144. [47] ANSYS CFX R13 Modelling Guide, 2010. [48] T. Esch, F.R. Menter, Heat transfer predictions based on two-equation turbulence models with advanced wall treatment, in: K. Hanjalic, Y. Nagano, M. Tummers (Eds.), Turbulence, Heat and Mass Transfer 4, K. Hanjalic, Y. Nagano and M. Tmmers (Editors), Begell House Inc.Antalya, Turkey, 2003, pp. 633-640. [49] B.A. Kader, Temperature and concentration profiles in fully turbulent boundary layers, Int. J. of Heat Mass Transfer 24 (1981) 1541-1544. [50] F.P. Incropera, D.P. DeWitt, Fundamentals of Heat and Mass Transfer, John Wiley and Sons Inc, New York, 1985. [51] H. Hanzawa, A.K. Melikow, P.O. Fanger, Airflow characteristics in the occupied zone of ventilated spaces, ASHRAE Transactions 93-1 (1987) 524-539. [52] H.A. Ersan, Numerical investigation of the effects of the outer turbulence in fluid flow and heat transfer characteristics, M.Sc. thesis from School of Natural and Applied Sciences, Mechanical Engineering Department, Uludag University, 2012. [53] Q. Chen, Ventilation performance prediction for buildings: A method overview and recent applications, Build. and Environ. 44 (2009) 848-858. [54] R. Yousaf, D. Wood, M. Cook, T. Yang, S. Hodder, D. Loveday, M. Passmore, CFD and PIV based investigation of indoor air flows dominated by buoyancy effects generated by human
13 Page 13 of 23
[59] [60] [61] [62]
ip t
[58]
cr
[57]
us
[56]
an
[55]
occupancy and equipment, in: Proceedings of BS2011 12th Conference of Int. Building Performance Simulation Association, Sydney, Australia, November 14-16, 2011, pp. 14651472. G.R. Iyer, S. Yavuzkurt, Comparison of low Reynolds number k-ε models in simulation of momentum and heat transport under high free stream turbulence, Int. J. of Heat Mass Transfer 42 (1999) 723-737. M. Cehlin, B.Moshfegh, Numerical modeling of a complex diffuser in a room with displacement ventilation, Build. and Environ. 45 (2010) 2240-2252. G. Cao, M. Ruponen, R. Paavilainen, J. Kurnitski, Modelling and simulation of the near-wall velocity of a turbulent ceiling attached plane jet after its impingement with the corner, Build. and Environ. 46 (2011) 489-500. M. Loomans, The measurement and simulation of indoor air flow, Ph.D. thesis, Eindhoven University of Technology, Eindhoven, The Netherlands, 1998. E. Lee, C.E. Feigley, J. Khan, An Investigation of inlet velocity in simulating the dispersion of indoor contaminants via Computational Fluid Dynamics, Ann. Occup. Hyg. 46(2002) 701-712. P.E. Hancock, P. Bradshaw, The effect of free-stream turbulence on turbulent boundary layers, ASME J. of Fluids Engineering 105 (1983), 284-289. P.V. Nielsen, A. Restivo, J.H. Whitelaw, Buoyancy-affected flows in ventilated rooms, Numerical Heat transfer 2 (1979) 115-127. F. Durrani, M.J. Cook, J.J. McGuirk, T. Chenvidyakarn, Simulating multiple steady states in naturally ventilated enclosures using Large Eddy Simulation, in: Proceedings of BS2013 13th Conference of Int. Building Performance Simulation Association, Chambery, France, August 26-28, 2013, pp. 275-281.
Ac ce
pt
ed
M
[63] S.P. O’Halloran, M.H. Hosni, B.T. Beck, T.P. Gielda, Scaled room three-dimensional computational fluid dynamics simulations, in: Proceedings of IMECE’02 ASME Int. Mechanical Engineering Congress and Exposition, New Orleans, Louisiana, November 17-22, 2002, pp. 1-5.
14 Page 14 of 23
Table(s) with Caption(s)
TABLES Table 1 Investigated computational cases. Tu
0.01
0.04
0.40
0.00156
0.0156
0.156
0.0156
0.00156
0.0156
0.156
Cases
1
2
3
4
5
6
7
ip t
LS
Ac
ce pt
ed
M
an
us
cr
Room Geometry Full Scale Model Annex20 Room
Page 15 of 23
Figure(s)
an
us
cr
ip t
FIGURES
Ac
ce pt
ed
M
Fig. 1.
Page 16 of 23
ip t cr
Ac
ce pt
ed
M
an
us
Fig. 2.
Page 17 of 23
ip t cr us an M ed Ac
ce pt
Fig. 3.
Page 18 of 23
ip t cr us
Ac
ce pt
ed
M
an
Fig. 4.
Page 19 of 23
ip t cr us an M ed ce pt Ac
Fig. 5.
Page 20 of 23
ip t cr us an M ed ce pt Ac
Fig. 6.
Page 21 of 23
List of Figure Captions
FIGURE CAPTIONS Fig. 1. Geometry and measurement line in the center plane (z = W/2) of IEA Annex20 room. Fig. 2. Mesh structure of IEA Annex20 room and dimensionless velocity profile comparison for different number of nodes with Nielsen’s measurements [8].
ip t
Fig. 3. Velocity vectors and isotherms for k-ε models with scalable wall function (SWF) and k-ω models with automatic near wall treatment (ANWT) a) IEA Annex20 Simulation by Lemaire (1991) [23] b) c) d) e) This study (Ar = 0.173).
cr
Fig. 4. Comparison of predicted streamlines a) Nielsen’s experimental study [8] b) RNG k-ε model (This study).
us
Fig. 5. Effects of inlet turbulence intensity (Tu) and length scale (LS) on velocity and temperature distribution for RNG k-ε model with scalable wall function (SWF) (I) Tu=0.01, (II) Tu=0.4 a) LS=0.00156, b ) LS=0.0156, c) LS=0.156.
Ac
ce pt
ed
M
an
Fig. 6. Comparison of double solution of mean velocity fields for RANS with RNG k-ε model with scalable wall function (SWF) (This study) c) d), and double solution for RANS with Std. k-ε model with wall function (WF) (Lemaire, 1991) [23] b), double solution for LES with dynamic approach and different geometry (Ezzouhri et. al 2009) [30] a).
Page 22 of 23
*Highlights (for review)
HIGHLIGHTS
Solution multiplicity encountered in ventilation is related to chaos theory.
-
Inlet turbulence parameters (Tu, LS) can cause solution multiplicity in ventilation.
-
An increase in Tu causes solution multiplicity at low and medium Length Scales (LSs).
-
LS variation has no solution multiplicity effect at a low turbulence intensity (Tu).
-
Possibility of solution multiplicity must be considered in any validation study.
Ac
ce pt
ed
M
an
us
cr
ip t
-
Page 23 of 23
S0378-7788(15)00094-8 http://dx.doi.org/doi:10.1016/j.enbuild.2015.01.067 ENB 5673
To appear in:
ENB
Received date: Revised date: Accepted date:
12-8-2014 12-1-2015 30-1-2015
Please cite this article as: E. Pulat, H.A. Ersan, Numerical Simulation of Turbulent Airflow in a Ventilated Room: Inlet Turbulence Parameters and Solution Multiplicity, Energy and Buildings (2015), http://dx.doi.org/10.1016/j.enbuild.2015.01.067 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Numerical Simulation of Turbulent Airflow in a Ventilated Room: Inlet Turbulence Parameters and Solution Multiplicity Uludag University, Faculty of Engineering, Mechanical Engineering Department, Gorukle Campus, TR-16059 Bursa, Turkey 2 Uludag University, School of Natural and Applied Sciences, Mechanical Engineering Department, Gorukle Campus, TR-16059 Bursa, Turkey * Corresponding Author: Tel.: +90 224 4428176, Fax: +90 224 4428021, E-Mail: [email protected]
ip t
1
Erhan PULAT1,* and Hıfzı Arda ERSAN2
cr
Abstract
ed
M
an
us
In this study, airflow and temperature distributions in the well-known International Energy Agency (IEA) Annex 20 room are predicted numerically to investigate the effects of the inlet turbulence intensity and the length scale on the flow characteristics, while considering the possibility of solution multiplicity related to chaos theory. The flow is considered to be turbulent, steady, incompressible, and two-dimensional. Computations are performed using the standard k-ε, RNG k-ε, standard k-ω, and SST k-ω turbulence models with scalable and automatic wall functions, and the results are compared with numerical and experimental results from the literature. The validated turbulence model is then used to investigate the effects of the turbulence intensity and the length scale. At a low inlet turbulence intensity value (Tu=0.01), the length scale variation has no influence on the flow pattern. However, the length scale affects the flow pattern at a high inlet turbulence intensity value (Tu=0.4). At constant low and medium length scale values, an increase in the inlet turbulence intensity from 0.01 to 0.4 affects the flow pattern. However, an increase in the turbulence intensity has no influence on the flow pattern at a constant high length scale value.
Ac ce
1. Introduction
pt
Keywords: Room Ventilation, Eddy Viscosity Turbulence Models, Turbulence Intensity and Length Scale, Solution Multiplicity, CFD
On average, people spend 90% of their life indoors [1]. Good airflow is important in ventilated spaces to provide good air quality and thermal comfort for the occupants, to maintain specified thermal and flow conditions for industrial processes and to ensure the efficient use of energy. The importance of ventilation related energy use is increasing and may represent up to 50% of the total energy use of a building, particularly for certain typologies such as office buildings [1]. Computational fluid dynamics (CFD) is one of the primary methods used to assess indoor airflow [2, 3, 4], among other methods such as the Building Energy Simulation (BES) approach and zonal models. In recent years, CFD has taken a prominent role in the simulation of indoor environment airflow problems [5, 6]. However, indoor airflow is a complex system, and new modeling approaches such as time series modeling can be found in the literature [7]. Pioneering work using the CFD technique to simulate room air motion was performed by Nielsen [8]. The flow inside a building environment is generally turbulent [2], and the IEA Annex 20 room, described in Nielsen [9], is one of the experimental setups for studying the suitability of room turbulent air distribution models [10, 11]. Aside from these studies, which were mainly focused on the comparison of turbulence models, and although the effects of turbulent boundary conditions in ventilated enclosures were identified in Reynolds-averaged 1 Page 1 of 23
Ac ce
pt
ed
M
an
us
cr
ip t
Navier-Stokes (RANS) simulations by Cao and Meyers [12], turbulent conditions have been largely neglected in simulations in the indoor ventilation literature [13, 14]. As observed by Abdilganie et al. [15] and Cao and Meyers [12], further investigation is required. Solution multiplicity refers to the existence of more than one solution for the overall flow patterns in a room under the same boundary conditions. Different initial conditions can lead to different solutions [16]. This phenomenon is in essence related to the chaotic behavior in nonlinear dynamic systems that can be encountered in all fields of science, including engineering sciences [17]. Zheng et al. [18] studied the nonlinear equations for ventilation flows in an enclosure with a single opening. Theoretical and computational nonlinear dynamic analysis of natural ventilation in a two-zone building was performed by Yang et al. [19, 20]. Van Schijndel [21] performed numerical research of the chaotic behavior of the airflow in an ordinary ventilated room. Chen [22] and Lemaire [23] reported two different airflow pattern results from simulations of the IEA Annex 20 room. The numerical study of Xin et al. [24] was intended to reveal some of the possible sequences of bifurcations leading to chaotic mixed convection in confined spaces. Recently, Dillon et al. [25] investigated chaotic natural convection in a tall rectangular cavity with non-isothermal walls. In addition to the above studies, Yee et al. [26] paid particular attention to isolating the different nonlinear behaviors and spurious dynamics caused by some of the numerical uncertainties observed in practical CFD computations. The multiplicity of airflow in naturally ventilated buildings [27, 28], smoke flows [29], and mixed convection in 3D ventilated cavities [30] has been confirmed by experiments (Müllejans, 1966 in [16], [30, 31, 32]). As seen in the above CFD and experimental studies, different initial conditions in CFD studies can lead to different solutions in both steady and transient calculations, and one stable solution can shift to another when there is sufficient perturbation. For example, in natural ventilation, the airflow pattern and the solution multiplicity are dependent on air change parameters [27, 28, 31]. In the case of a ventilated room, chaotic computational behavior is observed when very small changes are made in the supply air temperature [21, 32]. In mixed convection in confined spaces, multiplicity is dependent on two main system parameters, the Re and Ar numbers [22, 23, 24]. Although turbulent inlet parameters were studied from the point of view of flow and heat transfer in some of the studies mentioned above [12, 13 14, 15], solution multiplicity was not considered. In addition to providing a comparison of the popular RANS-based eddy viscosity turbulence models incorporating the scalable wall function and automatic near-wall treatment approaches, the main objective of this study is to investigate the effects of the inlet turbulence intensity and the length scale on the flow characteristics of the IEA Annex 20 room while considering the possibility of solution multiplicity. After a brief explanation of turbulence models and near-wall treatment methods is provided, a validated RNG k-ε turbulence model with a scalable wall function is used for this objective. 2. Materials and methods 2.1. Geometry The geometry of the IEA Annex 20 room described by Nielsen [9] is shown in Figure 1. The dimensions of the two-dimensional room are H=3 m, L=9 m, h=0.168 m and t=0.48 m, and it is equipped with an inlet and outlet to avoid numerical problems, as in Voigt’s study [10]. In a recent study, Olmedo and Nielsen [33] checked the two dimensionality of the test case, and another aim of their study is to analyze the results obtained considering the case as steady or unsteady. It has been concluded that the IEA Annex20 case solved with time
2 Page 2 of 23
dependent equations and k-epsilon model gives a solution similar to the solution obtained with steady state three dimensional and the two dimensional equations. 2.2. Turbulence models and near-wall treatment methods
Ac ce
pt
ed
M
an
us
cr
ip t
Computations were performed using the standard k-ε [34], RNG k-ε [35], standard k-ω [36] and SST k-ω [37] turbulence models. These two-equation models with eddy viscosity are widely used, as they offer a good compromise between numerical effort and computational accuracy. The advantage of the standard k-ε model from a numerical point of view is that it is robust and reliable. Although the physics is treated in a simplistic manner, the model works surprisingly well in many types of flows [38]. In a recent study, Klinzing and Sparrow [39] compared five turbulence models commonly offered in commercially available CFD software, including FLUENT, ANSYS CFX, and STAR-CD. They found that the most suitable turbulence model for external flows was the standard k-ε model with enhanced nearwall treatment. Chen [40] previously evaluated the performance of five k-ε models for the simulation of forced-convection, natural-convection, and mixed-convection room airflow. The results from the RNG k-ε model were slightly better than the standard k-ε model, and the other models, including the low Reynolds number k-ε, the two-layer k-ε, and the two-scale kε models, were not stable. Therefore, Chen recommended the RNG k-ε model for indoor airflow simulations. The k-ω model was compared with the standard and low Reynolds number k-ε models to model recirculating flows in a ventilated space by Peng et al. [41]. To improve the prediction accuracy of the standard k-ω model, they proposed some modifications. The SST k-ω model is a modification of the k-ω model, but it predicts a recirculation zone larger than that measured in PIV experiments in a water scale model of the Annex 20 room [10]. The k-ε and k-ω models seek to solve a modified set of transport equations by introducing mean and fluctuating components ( U i U i u i ). The continuity equation is not altered, and the momentum and scalar transport equations contain turbulent flux terms in addition to the molecular diffusion fluxes. These are the Reynolds stresses, u i u j . These terms arise from the nonlinear convective term in the unaveraged equations. The resulting steady-state RANS equations, assuming air is Newtonian and incompressible, and dropping the bar for averaged quantities, are as follows: Continuity: U j 0 x j Momentum: U iU x j
j
(1)
p ij u i u xi x j
j
S
(2)
M
Energy: T ( U j h tot ) ( u j h ) U i ( ij u i u j ) S E x j x j x j
(3)
where is the molecular stress tensor. S M is a source term added for buoyancy calculations:
3 Page 3 of 23
S M ,buoy ( ref ) g
(4)
The density difference is evaluated using the Boussinesq model because the density variation is driven only by small temperature variations in the room ventilation. The buoyancy source term is approximated as
ref ref (T Tref )
(5)
1 T
(6)
p
cr
ip t
where is the thermal expansivity
us
The mean total enthalpy is given by
an
1 htot h U iU i k 2
(7)
The total enthalpy contains a contribution from the turbulent kinetic energy, k, given by 1 2 ui 2
M
k
(8)
ed
Eddy viscosity models use the gradient diffusion hypothesis to relate the Reynolds stresses to the mean velocity gradients and the eddy (turbulent) viscosity (the Boussinesq approach) [41, 42]. 2 ij k t U k 3 x k
pt
U i U j ui u j t x xi j
(9)
Ac ce
where t is the eddy viscosity that is linked to the turbulent kinetic energy and its dissipation rate via the relation
t C
k2
(10)
The main advantage of this approach from the observation that the representation of Reynolds stress term is of exactly the same form as that of the diffusion terms in the original equations. They can be combined if an effective viscosity is defined as the sum of the laminar viscosity and the turbulent viscosity,
eff t
(11)
In two-equation models, the turbulence velocity scale is computed from the turbulent kinetic energy, which is obtained from the solution of its transport equation. The turbulence length scale is estimated from two properties of the turbulence field, usually the turbulent kinetic 4 Page 4 of 23
energy and its dissipation rate. The dissipation rate of the turbulent kinetic energy is obtained from its transport equation. Standard versions of four models used have been explained briefly as follows. Standard k-ε model
ip t
The k-ε model introduces two new variables into the system of equations. In this model, C is a constant, 0.09. The values of k and ε come directly from the differential transport equations for the turbulent kinetic energy and the turbulence dissipation rate [34].
cr
RNG k-ε model
us
The RNG k-ε model is based on a renormalization group analysis of the Navier-Stokes equations. The transport equations for turbulence generation and dissipation are the same as those for the standard k-ε model, but the model constants differ, with the constant C 1 replaced by the function C 1RNG [35].
an
Standard k-ω model
k
(12)
ed
t
M
One of the advantages of the k-ω formulation is the near-wall treatment for low Reynolds number computations. The k-ω models assume that the turbulence viscosity is linked to the turbulent kinetic energy and turbulent frequency via the relation
pt
The k-ω model solves two transport equations, one for the turbulent kinetic energy, k, and one for the turbulent frequency, ω [36]. The Shear Stress Transport (SST) k-ω model
Ac ce
The SST k-ω model accounts for the transport of the turbulent shear stress and gives highly accurate predictions of the onset and amount of flow separation under adverse pressure gradients [37]. The standard k-ω model does not account for the transport of the turbulent shear stress. This results in an overprediction of the eddy viscosity. The proper transport behavior can be obtained by incorporating a limiter in the formulation of the eddy viscosity as follows:
t a1 k max(a1 , SF2 )
(13)
F2 is a blending function that restricts the limiter to the wall boundary layer, because the underlying assumptions are not correct for free shear flows. S is an invariant measure of the strain rate. The blending function is critical to the success of the method. Its formulation is based on the distance to the nearest surface and on the flow variables:
5 Page 5 of 23
F2 tanh arg 22
with
2 k 500 arg 2 max ' , 2 y y
(14)
Near-wall treatment methods
Ac ce
pt
ed
M
an
us
cr
ip t
When computing wall-bounded turbulent flows it is necessary to alter the turbulence model in the vicinity of the wall. Two basic approaches are generally used: the so-called wall function approach and the low Reynolds number model [43]. Near-wall treatment is crucial in any turbulence model, and the most popular approach is the wall function approach, due to its simplicity. Basically, the method requires that the first grid point near the wall be located in the fully turbulent region. In the low Reynolds number approach, the basic idea is to integrate the equations up to the wall. This approach requires the clustering of grid points near the wall, which increases the total number of grid points required for an accurate solution. In flows with steep pressure gradients or swirl, the wall function representation of the near-wall region is not sufficiently universal to allow accurate predictions [43]. However, according to Chow and Gao [44], the standard k-ε model provides steady numerical results without fluctuations and might be a suitable method for representing the separated and reattached flows if the wall functions are applied properly. Additionally, Mohammadi and Medic [45] have shown that the careful application of numerical methods and adequate numerical dissipation lead to a substantial improvement in the results of the classical k-ε model. A useful discussion of wall function treatments, including the mixed convection effect, is provided by Craft et al. [46]. An extension of this method is used in ANSYS CFX as follows. The main problem in the wall function approach is the danger of excessive mesh refinement, so that the physical assumptions are especially problematic in flows at lower Reynolds numbers in the near-wall mesh. To avoid this problem, the scalable wall function (SWF) formulation in ANSYS CFX is used in the k-ε model. The basic idea behind this approach is to limit the y+ value used in the logarithmic formulation by a lower bound of y͂+ = max(y+, 11.06). The computed y͂+ is not allowed to fall below this limit. However, in flows at low Reynolds numbers, this can cause an error in the displacement thickness, and it is therefore desirable to offer a formulation that will automatically switch from a wall function to a low Reynolds number near-wall formulation as the mesh is refined. The main idea behind this formulation, called the automatic near-wall treatment (ANWT), is to blend the wall value for ω between the logarithmic and the near-wall formulation. This method has been used in k-ω models [47]. A blending function that depends on y+ can therefore be defined that can be used to derive a formulation describing the behavior of ω in linear and logarithmic near-wall regions as follows [42, 47, 48]: 2 1 y vis y log2 y
(15)
A similar formulation is used to describe the relation between the velocity at the first grid point and the wall shear stress:
u
u 4 uVis
4
log 4
(16)
In both the wall function and the automatic near-wall treatment approaches, in the treatment of the energy equation near the wall, an algebraic formulation is required to link the temperature and the heat flux. The formulation of Kader [49] is used in ANSYS-CFX: 6 Page 6 of 23
Pr y e 2.12 ln 1 y Pr e
1
(17)
2
with
(18) (19)
ip t
1 where (Pr) 3.85 Pr 3 1.3 2.12 ln(Pr) 4 0.01(Pr y ) , 1 5 Pr 3 y
and the non-dimensional temperature is defined as
2.3. Boundary conditions and computational procedures
cr
qw Tw T , T c p u T
(20)
us
ghT0 u0
pt
Ar
ed
M
an
The modeled inlet conditions were a uniform velocity of 0.455 m/s, an inlet turbulence intensity of 4 % and a temperature of 20 ºC, consistent with Nielsen’s test case [9]. The turbulence length scale used was LS=0.0156, which is very close to Nielsen’s test case value of 0.0168. The outlet condition was a constant pressure of zero. The velocity was set to zero on the room walls and the inlet and outlet walls to impose a no-slip condition. Based on the conditions in the supply opening, the inlet Reynolds number was 5000. Air was treated as an ideal gas, and its density and viscosity were assumed to be 1.116 kg/m3 and 1.725 Pa·s [50]. The left and right wall and the ceiling were assumed to be adiabatic. A constant heat flux of 63.08 W/m2 was added along the floor in the non-isothermal case, similar to Lemaire’s study [23]. The effect of buoyancy is commonly expressed in the ventilation literature by the Archimedes number, Ar: (21)
Ac ce
where ΔT0 is the temperature difference between the inlet and outlet, and h and u0 are the height and velocity of the inlet, respectively. Archimedes number was 0.173 in the present study. The characteristics of turbulent flow are also commonly expressed by the turbulence intensity, the turbulence length scales, and the turbulent kinetic energy and its dissipation rate. The turbulence intensity, Tu, is the standard deviation divided by the mean velocity ( u ): Tu
u' x100% u
(22)
where the velocity fluctuations, u ' , are in the main direction of the flow. The turbulence length scales are the integral scale and the microscale. It is assumed that the turbulent motion consists of the superposition of eddies of various sizes. The integral scale identifies the average size of the largest eddies, whereas the microscale is a measure of the smallest eddies, which are mainly responsible for dissipation [51]. The ANSYS CFX R13 code uses the finite volume method, and pressure and velocity are coupled using the SIMPLEC algorithm with Rie-Chow interpolation. Buoyancy is modeled 7 Page 7 of 23
using the Boussinesq approach. A mesh independency study was performed by comparing three different mapped mesh structures (coarse, medium, and fine), corresponding to 12168, 19506 and 29054 nodes, using the standard k-ε model. The mesh structure with 29054 nodes was chosen by comparing the u/uo velocity profiles for the different numbers of nodes, as seen in Figure 2. Further information about mesh structures can be found in [52]. 3. Results and discussion
Ac ce
pt
ed
M
an
us
cr
ip t
In Figure 3, the velocity vectors and isotherms obtained under the specified boundary conditions from the popular RANS-based eddy viscosity turbulence models using the relatively new near-wall modeling approaches, namely the scalable wall function (SWF) and the automatic near-wall treatment (ANWT), are compared with the results of Lemaire’s [23] numerical study, because that study includes the non-isothermal case for the IEA Annex 20 room. As shown in Figure 3, the flow patterns obtained from the standard and RNG k-ε models are consistent with Lemaire’s flow pattern, that is, the main vortex is in the clockwise direction. In the case of the standard k-ω and SST k-ω models, however, the flow pattern changes completely: the cold jet falls down, causing a reversal of the main vortex. Due to this different flow pattern, the isotherms in the k-ω models are also not consistent with Lemaire’s isotherms. Additionally, the isotherms are predicted more accurately by the RNG k-ε model than the standard k-ε model when compared with Lemaire’s results, as shown in Figure 3a, 3b, and 3c. The streamlines in the RNG k-ε model are also qualitatively consistent with Nielsen’s [8] experimental results, as seen in Figure 4. For this reason, the effects of the inlet turbulence intensity and the length scale on the flow pattern and temperature distribution have been investigated with the RNG k-ε model. Additionally, as mentioned previously, this model was recommended for indoor airflow simulations by Chen [40]. Although one model might perform well for one flow but poorly for another, the performance of the RNG k-ε model was relatively stable [53]. Yousaf et al. [54] used ANSYS CFX to apply the SST k-ω and RNG kε models, using the same near-wall modeling approach as the present study, to more complex indoor airflows dominated by buoyancy effects generated by human occupancy and equipment. They reported differences of 6% and 8% in the temperature values near the inlet in the SST and RNG models, respectively; at the other sensor points, the RNG model predictions differed by less than 2.6%. The SST model appears to overestimate the temperature field. To investigate the effects of the turbulence intensity (Tu) and length scale (LS) on the flow pattern and temperature distribution, two Tu values (Tu = 0.01 and 0.4) and three LS values (LS = 0.00156, 0.0156 and 0.156) are considered, as shown in Figure 5. The seven cases shown in Table 1 are investigated using the RNG k-ε model. Case 4 is the comparison condition, representing the test case boundary conditions described in section 2.3. If the turbulence intensity is greater than 0.05, it is commonly assumed that the flow is in the high free-stream turbulence condition [55]. Although Tu value of 0.4 can be treated as very high for ventilation flows, high Tu values are encountered in air supply devices such as diffusers. For example, Cehlin and Moshfegh [56] was reported that turbulence intensity in the openings was around 0.1-0.3% while the intensity outside the diffuser in the near zone was up to 0.6. According to Cao et al. [57], in a ventilated room, due to wall effects, the magnitude of turbulence intensity can grow to a value greater than 0.2 from 0.02. In addition, the levels of inlet turbulence intensities investigated by Cao and Meyers [12] were 0.05 and 0.3, by Loomans [58] were 0.1 and 0.35, and by Joubert et al. [13] were 0.04 and 0.374. So high Tu values can be treated as reasonable, and as mentioned by van Hooff et al. [3], the flow pattern encountered in room ventilation can be challenging due to the presence of transitional flow, turbulence anisotropy, and adverse pressure gradient. Really, as mentioned by Cehlin and 8 Page 8 of 23
Ac ce
pt
ed
M
an
us
cr
ip t
Moshfegh [56] prediction of room airflows is a challenging task because of the complex nature of turbulence. Increase of turbulent intensity at the inlet also leads to an increased turbulent intensity in the rest of domain as shown by Cao and Meyers [12]. As known turbulence increases heat transfer rate and high Tu values may cause high heat transfer coefficients in ventilation studies together with room heating. From built environment point of view, understanding airborne contaminant transport and distribution in buildings is essential for addressing issues such as home and office air quality, worker exposure in manufacturing workrooms, product protection in clean rooms, accident response planning in nuclear facilities, avoiding fatalities in confined spaces and infection control in medical settings [59], and as mentioned by Cao and Meyers [12], when pollutant dispersion is also considered, the inlet turbulent quantities both turbulent intensity and length scales have an important effect on pollutant dispersion. As shown in Figure 5 (I), there is no influence of the length scale on the flow pattern or the temperature distribution at the low inlet turbulence intensity value (Tu=0.01). Although the length scale varies by a factor of a hundred (from 0.00156 to 0.156), the flow pattern and temperature distribution are not affected by this increase. However, at the high turbulence intensity value (Tu=0.4), at the low and medium length scale values (LS=0.00156 and 0.0156), the flow pattern changes completely, producing flow patterns similar to those obtained by the standard and SST k-ω models for the test case, as shown in Figure 3d and e. The cold jet falls down, causing a reversal of the main vortex. At the high length scale value (LS=0.156), the clockwise flow pattern is obtained again, and due to this change in the flow pattern, the temperature distribution is also clockwise. A similar flow pattern change was reported by Chen [22] and Lemaire [23] as a result of a small change in the Archimedes number (Ar). Their transition Archimedes numbers are different due to their different inlet sizes and Reynolds numbers [22]. This phenomenon was also reported by Heiselberg et al. [31] and Ezzouhri et al. [30] for different room dimensions, and was ascribed to a hysteresis effect. Ezzouhri et al. [30] concluded that the numerical bifurcation between the two flow structures was highly dependent on the turbulence level at the inlet. In the present study, at fixed small and medium length scale values, as the inlet turbulence intensity is increased from 0.01 to 0.4, the clockwise flow pattern is converted into a counterclockwise flow. At the high fixed inlet turbulence intensity value, if the length scale is decreased from 0.156 to 0.0156 or 0.00156, the clockwise flow pattern is also converted into a counterclockwise flow, as shown in Figure 5. These results are consistent with the findings of Hancock and Bradshaw [60], who showed for the first time that the free stream turbulence length scale was also an important parameter, in addition to the free stream turbulence intensity [55], and showed that the free stream turbulence effect decreased as the length scale increased. Really, if the length scale is high (LS=0.156) as seen from Figure 5c, flow pattern is not affected from increase in the free stream turbulence intensity from 0.01 to 0.4. A comparison of the present results with Lemaire’s [23] results for the same geometry using the standard k-ε model and Ezzouhri et al.’s [30] results for a different geometry using the LES approach is given in Figure 6. Additionally, Ezzouhri et al. [30] observed that the dynamic LES model not only correctly predicts the mean flow characteristics of the flow but is also able to correctly reproduce the flow bifurcation and the hysteresis effect. However, in the present study, as well as in Chen’s study [22] using the Lam-Bremhorst low Reynolds number k-ε model, Lemaire’s study [23] using the standard k-ε model, and Nielsen et al.’s study [61] using the standard k-ε model, it has been shown that RANS simulations have capability comparable to LES simulations. Similarly, in the recent study by Durrani et al. [62], the RNG k-ε model was used in addition to the LES model to represent turbulence in an unsteady RANS (URANS) simulation. They concluded that the URANS method was unable to capture the details of the flow structures
9 Page 9 of 23
pt
4. Conclusions
ed
M
an
us
cr
ip t
predicted by LES, but by its nature, the LES approach required approximately 500% more computation time than URANS. A better understanding of solution multiplicity in ventilation flows is very important because researchers and practitioners may encounter this phenomenon. For example, in the study of O’Halloran et al. [63], two turbulence models, the RNG k-ε model and the Reynolds stress model (RSM), were evaluated in simulations and compared with experimental particle image velocimeter (PIV) data. In their study, although qualitatively the RNG k-ε model and RSM results compared well with each other and with the PIV results, one difference was apparent in the RNG k-ε model result. They expressed this difference as follows: “In the recirculation area near the inlet of the jet, from approximately x=0 to 150 mm, the direction of flow is opposite to that of the PIV results and the RSM model results. The flow reversal appears to result from a horseshoe vortex that wraps around the inlet flow. The numerical reason for this difference could not be determined, neither by the researchers nor by the software manufacturer, and further investigation is needed to understand the cause of the differences.” This difference might be related to solution multiplicity. As shown in Figure 3, different turbulence models produce different solutions under the same boundary conditions. Additionally, as shown in Figure 5, a single turbulence model can also produce different solutions under different boundary conditions. However, as reported by Nielsen et al. [61] and Li and Nielsen [16], the existence of two solutions under the same boundary conditions is not only found experimentally, as in the results of Müllejans for the same Ar and Re numbers [Müllejans, 1966 in [16]]; two solutions can be obtained numerically from a single turbulence model for the same Ar and Re numbers [54, 61]. Consequently, in any validation process in the evaluation of a suitable turbulence model, one must consider whether the flow is expected to have multiple solutions. For example, in this study, the flow pattern predictions from the standard and SST k-ω models producing counterclockwise main vortices represent the second solution in each case, and the predictions of these models must be compared using the measurements from this second solution, if available, before making a final decision about their suitability.
Ac ce
Airflow and temperature distributions in the well-known IEA Annex 20 room have been predicted numerically to investigate the effects of the inlet turbulence intensity and the length scale on the flow characteristics while considering solution multiplicity. For this purpose, the RNG k-ε model is preferred because of its slightly better prediction of the temperature distribution, although similar flow patterns are obtained from the standard k-ε model. The inlet turbulence parameters have a great influence on the flow pattern, from the viewpoint of solution multiplicity, at similar Reynolds and Archimedes numbers, and the free stream turbulence effect decreases as the length scale increases. But to determine a threshold value of Tu that the flow pattern changes from cw direction to ccw direction, it is required consideration of a few additional Tu values. Such computations is planned as a future study. In addition, it has been shown that present steady RANS simulations have capability comparable to LES simulations required more computation time than Unsteady RANS in literature. Consequently, in any indoor airflow process or validation study, one must consider whether the flow is expected to have multiple solutions.
Nomenclature β
thermal expansivity [1/K]
cp
constant pressure specific heat [J/kgK] 10 Page 10 of 23
dissipation rate of turbulent kinetic energy [m2/s3]
g
gravitational acceleration [m/s2]
h
enthalpy [J/kg]
k
turbulent kinetic energy [m2/s2]
λ
thermal conductivity [W/mK]
μ
viscosity [kg/ms]
ν
kinematic viscosity [m2/s]
p
pressure [Pa]
qw
wall heat flux [W/m2]
ρ
density [kg/m3]
S
source term
T
temperature [K]
ω
turbulent frequency (ε / k) [1/s]
uτ
dimensionless friction velocity
U
velocity [m/s]
y+
dimensionless distance from the wall
cr us an M
ed
pt
References
ip t
ε
Working Group 16, Ventilation, good indoor air quality and rational use of energy, European Commission Joint Research Centre, Institute For Health and Consumer Protection, Physical and Chemical Exposure Unit, Report No 23, EUR 20741 EN, 2003.
[2]
Q. Chen, Computational Fluid Dynamics for HVAC: Successes and failures, ASHRAE Transactions 103-1 (1997) 178-187. T. van Hoof, B. Blocken, G.J.F. van Heijst, On the suitability of steady RANS CFD for forced mixing ventilation at transitional slot Reynolds numbers, Indoor Air 23 (2013) 236-249. E. Djunaedy, J.L.M. Hensen, M.G.L.C. Loomans, Development of a guideline for selecting a simulation tool for airflow prediction, in: Proceedings of the 8th IBPSA conference BuildingSimulation’2003, Eindhoven, 2003, pp. 267-274. R.A. Pitarma, J.E. Ramos, M.E. Ferreira, M.G. Carvalho, Computational fluid dynamics an advanced active tool in environmental management and education, Manag. of Environ. Qual.: An Int. J. 15 (2004) 102-110. P.V. Nielsen, Computational fluid dynamics and room air movement, Indoor Air 14-Suppl 7 (2004) 134-143. M. Maciejewska, A. Szczurek, G. Sikora, A.Wylomanska, Diffusive and subdiffusive dynamics of indoor microclimate: A time series modeling, Phys. Rev. E 86 (2012) 031128-1-11. P.V. Nielsen, Flow in air conditioned rooms, PhD thesis from the Technical University of Denmark, 1974.
[3] [4] [5] [6] [7] [8]
Ac ce
[1]
11 Page 11 of 23
[15] [16] [17] [18] [19] [20] [21] [22] [23]
[24] [25] [26] [27] [28] [29] [30]
ip t
cr
us
[14]
an
[13]
M
[12]
ed
[11]
pt
[10]
P.V. Nielsen, Specification of a two-dimensional test case – International Energy Agency, energy conservation in buildings and community systems, Annex20: Air flow pattern within buildings, Aalborg University Technical Report, 1990. L.K. Voigt, J.N. Sorensen, J.M. Pedersen, M. Brons, Review of four turbulence models using topology, in: Proceedings of the 8th IBPSA conference BuildingSimulation’2003, Eindhoven, 2003, pp. 1325-1331. M. Omri, N. Galanis, Evaluation of confined natural and forced convection predictions by different turbulence models, Int. J. of Numer. Methods for Heat and Fluid Flow 19 (2009) 5-24. S-J. Cao, J. Meyers, Influence of turbulent boundary conditions on RANS simulations of pollutant dispersion in mechanically ventilated enclosures with transitional slot Reynolds number, Build. and Environ. 59 (2013) 397-407. P. Joubert, A. Sandu, C. Beghein, F. Allard, Numerical study of the influence of inlet boundary conditions on the air movement in a ventilated enclosure, in: Proceedings of the 5th International Conference on Air Distribution in Rooms Roomvent’96, Yokohama,1996, pp. 235-242. S.A. Al-Sanea, M.F. Zedan, M.B. Al-Harbi, Effect of supply Reynolds number and room aspect ratio on flow and ceiling heat-transfer coefficient for mixing ventilation, Int. J. of Therm. Sciences 54 (2012) 176-187. A.M. Abdilghanie, L.R. Collins, D.A. Caughey, Comparison of turbulence modeling strategies for indoor flows, ASME J. of Fluids Engineering 131(2009), 051402-1-18. Y. Li, P.V. Nielsen, CFD and ventilation research, Indoor Air 21(2011), 442-453. G. Baker, F.A. McRobie, J.M.T. Thompson, Implications of chaos theory for engineering science, IMechE J. of Mechanical Engineering Sci Part C 211 (1997) 349-363. Z. Zheng, P. Xu, Y. Li, G. Chen, Nonlinear resonance and quasi-periodic solutions for ventilation flows in a single opening enclosure, Int. J. Bifurcation and Chaos 15 (2005) 18011808. L. Yang, P. Xu, Y. Li, Nonlinear Dynamic Analysis of natural ventilation in a two-zone building: Part A-Theoretical analysis, ASHRAE HVAC&R Research 12 (2006) 231-255. L. Yang, P. Xu, Y. Li, Nonlinear Dynamic Analysis of natural ventilation in a two-zone building: Part B-CFD simulations, ASHRAE HVAC&R Research 12 (2006) 257-278. A.W.M. van Schijndel, Chaotic behavior of the airflow in a ventilated room, in: Proceedings of the COMSOL Conference, Milan,2009, pp. -.(http://www.comsol.nl/paper/chaotic-behavior-ofthe-airflow-in-a-ventilated-room-6885 (Accessed in Nov. 18, 2013) Q. Chen, Int. Energy Agency, Energy conservation in building and community systems, Annex 20: Air flow patterns within buildings, Subtask 1: Room Air and Contaminant Flow, Research Item 1.46 Simulation of simple test cases, Report No AN20.1-CH-90-ETHZ15, 4 March 1991. A.D. Lemaire, Int. Energy Agency, Energy conservation in building and community systems, Annex 20: Air flow patterns within buildings, Subtask 1: Room Air and Contaminant Flow, Research Item No 1.46NL Simulation of simple test cases 2D1,2D2, Report No AN20.1-NL-90TNO-TPD10, May 1991. G. Xin, W.K. Chow, S. Liu, Windows and routes to chaos in mixed convection in confined spaces, Chaos Solitons Fractals 15 (2003) 543-558. H. Dillon, A. Emery, A. Mescher, Analysis of chaotic natural convection in a tall rectangular cavity with non-isothermal walls, Frontiers in Heat Mass Transfer 4 (2013) 1-9. H.C. Yee, J.R. Torczynski, S.A. Morton, M.R. Visbal, P.K. Sweby, On spurious behavior of CFD simulations, Int. J. For Num. Methods in Fluids 30 (1999) 675-711. Y. Li, A. Delsante, Natural ventilation induced by combined wind and thermal forces, Build. and Environ. 36 (2001) 59-71. Y. Li, A. Delsante, Z. Chen, M. Sandberg, A. Andersen, M. Bjerre, P. Heiselberg, Some examples of solution multiplicity in natural ventilation, Build. and Environ. 36 (2001) 851-858. J. Gong, Y. Li, Solution Multiplicity of smoke flows in a simple building, in: Fire Safety Science Proceedings of the 9th International IAFSS Symposium, Karlsruhe, 2008, pp. 895-906. R. Ezzouhri, P. Jaubert, F. Penot, S. Mergui, Large eddy simulation of turbulent mixed convection in a 3D ventilated cavity: Comparison with existing data, Int. J. of Therm. Sciences 48 (2009) 2017-2024.
Ac ce
[9]
12 Page 12 of 23
Ac ce
pt
ed
M
an
us
cr
ip t
[31] P. Heiselberg, Y. Li, A. Andersen, M. Bjerre, Z. Chen, Experimental and CFD evidence of multiple solutions in a naturally ventilated building, Indoor Air 14 (2004) 43-54. [32] J. van Schijndel, Simulation and experimental validation of chaotic behavior of airflow in a ventilated room, in: Proceedings of the 9th Nordic Symposium on Building Physics NSB 2011, Tampere,2011, vol.2, pp. 235-242.(http://alexandria.tue.nl/openaccess/Metis250706.pdf) (Accessed in Nov. 18, 2013) [33] I. Olmedo, P.V. Nielsen, Analysis of the IEA 2D test. 2D, 3D, steady or unsteady airflow? Aalborg University Department of Civil Engineering DCE Technical Report No.106, 2010. [34] B.E. Launder, D.B. Spalding, Lectures in Mathematical Models of Turbulence, Academic Press, London, 1972. [35] V. Yakhot, S.A. Orszag, Renormalization group analysis of turbulence I. Basic theory, J. Sci. Comput. 1 (1986) 1-51. [36] D.C. Wilcox, Reassessment of the scale determining equation for advance turbulence models, AIAA J. 26 (1988) 1299-1310. [37] D.C. Wilcox, Turbulence Modeling for CFD, DCW Industries Inc,La Canada, CA, 1994. [38] L. Davidson, P.V. Nielsen, Large Eddy Simulations of the flow in a three dimensional ventilated room, in: Proceedings of the 5th International Conference on Air Distribution in Rooms Roomvent’96, Yokohama,1996, pp. 161-2168. [39] W.P. Klinzing, E.M. Sparrow, Evaluation of turbulence models for external flows, Numerical Heat Transfer Part A Applications 55 (2009) 205-228. [40] Q. Chen, Comparison of different k-ε models for indoor airflow computation, Numerical Heat Transfer Part B Fundamentals 28 (1995) 353-369. [41] S-H. Peng, L. Davidson, S. Holmberg, The two equation turbulence k-ω model applied to recirculating ventilation flows, Thermo and Fluid Dynamics, Chalmers University of Technology, Göteborg: Report 96/13. [42] ANSYS CFX-Solver Theory Guide, 2009. [43] W.C. Lasher, D.B. Taulbee, The effect of the near-wall treatment on the behavior of ReynoldsStress turbulence models for reattaching flows, in: Seventh Symposium on Turbulent Shear Flows, Stanford University, U.S.A., August 21-23, 1989, pp. 10.5.1-10.5.6. [44] W.K. Chow, Y. Gao, On using the wall function and the staggered grid arrangement at upwind sharp convex corners in simulating approach airflow, Int. J. Comput. Fluid Dyn. 19 (2005) 381398. [45] B. Mohammadi, G. Medic, A critical evaluation of the classical k-ε model and wall-laws for unsteady flows over bluff bodies, Int. J. Comput. Fluid Dyn. 10 (1998) 1-11. [46] T.J. Craft, S.E. Gant, A.V. Gerasimov, H. Iacovides, B.E. Launder, Development and application of wall-function treatments for turbulent forced and mixed convection flows, Fluid Dyn. Res. 38 (2006) 127-144. [47] ANSYS CFX R13 Modelling Guide, 2010. [48] T. Esch, F.R. Menter, Heat transfer predictions based on two-equation turbulence models with advanced wall treatment, in: K. Hanjalic, Y. Nagano, M. Tummers (Eds.), Turbulence, Heat and Mass Transfer 4, K. Hanjalic, Y. Nagano and M. Tmmers (Editors), Begell House Inc.Antalya, Turkey, 2003, pp. 633-640. [49] B.A. Kader, Temperature and concentration profiles in fully turbulent boundary layers, Int. J. of Heat Mass Transfer 24 (1981) 1541-1544. [50] F.P. Incropera, D.P. DeWitt, Fundamentals of Heat and Mass Transfer, John Wiley and Sons Inc, New York, 1985. [51] H. Hanzawa, A.K. Melikow, P.O. Fanger, Airflow characteristics in the occupied zone of ventilated spaces, ASHRAE Transactions 93-1 (1987) 524-539. [52] H.A. Ersan, Numerical investigation of the effects of the outer turbulence in fluid flow and heat transfer characteristics, M.Sc. thesis from School of Natural and Applied Sciences, Mechanical Engineering Department, Uludag University, 2012. [53] Q. Chen, Ventilation performance prediction for buildings: A method overview and recent applications, Build. and Environ. 44 (2009) 848-858. [54] R. Yousaf, D. Wood, M. Cook, T. Yang, S. Hodder, D. Loveday, M. Passmore, CFD and PIV based investigation of indoor air flows dominated by buoyancy effects generated by human
13 Page 13 of 23
[59] [60] [61] [62]
ip t
[58]
cr
[57]
us
[56]
an
[55]
occupancy and equipment, in: Proceedings of BS2011 12th Conference of Int. Building Performance Simulation Association, Sydney, Australia, November 14-16, 2011, pp. 14651472. G.R. Iyer, S. Yavuzkurt, Comparison of low Reynolds number k-ε models in simulation of momentum and heat transport under high free stream turbulence, Int. J. of Heat Mass Transfer 42 (1999) 723-737. M. Cehlin, B.Moshfegh, Numerical modeling of a complex diffuser in a room with displacement ventilation, Build. and Environ. 45 (2010) 2240-2252. G. Cao, M. Ruponen, R. Paavilainen, J. Kurnitski, Modelling and simulation of the near-wall velocity of a turbulent ceiling attached plane jet after its impingement with the corner, Build. and Environ. 46 (2011) 489-500. M. Loomans, The measurement and simulation of indoor air flow, Ph.D. thesis, Eindhoven University of Technology, Eindhoven, The Netherlands, 1998. E. Lee, C.E. Feigley, J. Khan, An Investigation of inlet velocity in simulating the dispersion of indoor contaminants via Computational Fluid Dynamics, Ann. Occup. Hyg. 46(2002) 701-712. P.E. Hancock, P. Bradshaw, The effect of free-stream turbulence on turbulent boundary layers, ASME J. of Fluids Engineering 105 (1983), 284-289. P.V. Nielsen, A. Restivo, J.H. Whitelaw, Buoyancy-affected flows in ventilated rooms, Numerical Heat transfer 2 (1979) 115-127. F. Durrani, M.J. Cook, J.J. McGuirk, T. Chenvidyakarn, Simulating multiple steady states in naturally ventilated enclosures using Large Eddy Simulation, in: Proceedings of BS2013 13th Conference of Int. Building Performance Simulation Association, Chambery, France, August 26-28, 2013, pp. 275-281.
Ac ce
pt
ed
M
[63] S.P. O’Halloran, M.H. Hosni, B.T. Beck, T.P. Gielda, Scaled room three-dimensional computational fluid dynamics simulations, in: Proceedings of IMECE’02 ASME Int. Mechanical Engineering Congress and Exposition, New Orleans, Louisiana, November 17-22, 2002, pp. 1-5.
14 Page 14 of 23
Table(s) with Caption(s)
TABLES Table 1 Investigated computational cases. Tu
0.01
0.04
0.40
0.00156
0.0156
0.156
0.0156
0.00156
0.0156
0.156
Cases
1
2
3
4
5
6
7
ip t
LS
Ac
ce pt
ed
M
an
us
cr
Room Geometry Full Scale Model Annex20 Room
Page 15 of 23
Figure(s)
an
us
cr
ip t
FIGURES
Ac
ce pt
ed
M
Fig. 1.
Page 16 of 23
ip t cr
Ac
ce pt
ed
M
an
us
Fig. 2.
Page 17 of 23
ip t cr us an M ed Ac
ce pt
Fig. 3.
Page 18 of 23
ip t cr us
Ac
ce pt
ed
M
an
Fig. 4.
Page 19 of 23
ip t cr us an M ed ce pt Ac
Fig. 5.
Page 20 of 23
ip t cr us an M ed ce pt Ac
Fig. 6.
Page 21 of 23
List of Figure Captions
FIGURE CAPTIONS Fig. 1. Geometry and measurement line in the center plane (z = W/2) of IEA Annex20 room. Fig. 2. Mesh structure of IEA Annex20 room and dimensionless velocity profile comparison for different number of nodes with Nielsen’s measurements [8].
ip t
Fig. 3. Velocity vectors and isotherms for k-ε models with scalable wall function (SWF) and k-ω models with automatic near wall treatment (ANWT) a) IEA Annex20 Simulation by Lemaire (1991) [23] b) c) d) e) This study (Ar = 0.173).
cr
Fig. 4. Comparison of predicted streamlines a) Nielsen’s experimental study [8] b) RNG k-ε model (This study).
us
Fig. 5. Effects of inlet turbulence intensity (Tu) and length scale (LS) on velocity and temperature distribution for RNG k-ε model with scalable wall function (SWF) (I) Tu=0.01, (II) Tu=0.4 a) LS=0.00156, b ) LS=0.0156, c) LS=0.156.
Ac
ce pt
ed
M
an
Fig. 6. Comparison of double solution of mean velocity fields for RANS with RNG k-ε model with scalable wall function (SWF) (This study) c) d), and double solution for RANS with Std. k-ε model with wall function (WF) (Lemaire, 1991) [23] b), double solution for LES with dynamic approach and different geometry (Ezzouhri et. al 2009) [30] a).
Page 22 of 23
*Highlights (for review)
HIGHLIGHTS
Solution multiplicity encountered in ventilation is related to chaos theory.
-
Inlet turbulence parameters (Tu, LS) can cause solution multiplicity in ventilation.
-
An increase in Tu causes solution multiplicity at low and medium Length Scales (LSs).
-
LS variation has no solution multiplicity effect at a low turbulence intensity (Tu).
-
Possibility of solution multiplicity must be considered in any validation study.
Ac
ce pt
ed
M
an
us
cr
ip t
-
Page 23 of 23