Evaluation of the eddy viscosity turbulence models for the simulation of convection–radiation coupled heat transfer in indoor environment

Energy & Buildings 184 (2019) 8–18

Contents lists available at ScienceDirect

Energy & Buildings journal homepage: www.elsevier.com/locate/enbuild

Evaluation of the eddy viscosity turbulence models for the simulation of convection–radiation coupled heat transfer in indoor environment Xiangdong Li a, Jiyuan Tu a,b,∗ a

School of Engineering, RMIT University, PO Box 71, Bundoora, VIC 3083, Australia Key Laboratory of Ministry of Education for Advanced Reactor Engineering and Safety, Institute of Nuclear and New Energy Technology, Tsinghua University, PO Box 1021, Beijing 100086, China

b

a r t i c l e

i n f o

Article history: Received 27 August 2018 Revised 21 November 2018 Accepted 22 November 2018 Available online 7 December 2018 Keywords: Eddy viscosity turbulence models Thermal radiation Thermal buoyancy flows Contaminant dispersion Indoor environment

a b s t r a c t This paper presents a numerical evaluation of the eddy viscosity turbulence models (the zero-equation, one-equation, k-ε group and k-ω group models) in terms of CFD modelling of convection-radiation coupled heat transfer in indoor environment. The results show that all the models except the zero-equation model have acceptable accuracies to predict the bulk air temperature. However, the zero-equation and one-equation models seriously overpredict the convective heat transfer coefficient at the heater surfaces while the k-ε group models severely underpredict it, ultimately resulting in the incorrect predictions of the thermal buoyancy flows and contaminant dispersion. On the contrary, the k-ω group models perform best to calculate the convection-radiation heat split, which is critical to the successful predictions of the near-wall flow field and distribution patterns of the contaminants. This study also demonstrates the importance of simultaneously considering the convective and radiative heat transfer mechanisms when evaluating the turbulence models for indoor environment simulations. © 2018 Elsevier B.V. All rights reserved.

1. Introduction Turbulence modelling is an indispensable and critical part of the computational fluid dynamics (CFD) simulation of indoor airflows. Turbulence is featured with the energy-carrying eddies with a wide range of length and time scales. In indoor environment, the eddies play a vital role in transporting heat and mass. The accuracy of turbulence modelling can largely determine the effectiveness of CFD-based evaluations of air distribution, thermal comfort and indoor air quality. The approaches to model turbulence can be classified into two big categories, as shown in Fig. 1. The first category relies on the direct solution of the Navier–Stokes equations. Depending on the spectrum of the turbulence eddy scales in which the Navier–Stokes equations are directly solved, there are direct numerical simulation (DNS), large eddy simulation (LES) and detached eddy simulation (DES). However, the direct solution of the Navier–Stokes equations requires very high mesh resolutions, making it unpractical to simulate realistic indoor environments due to the extremely high computational cost. Alternatively, the Reynolds averaged Navier–Stokes (RANS) approach only gives time-averaged solutions to the Navier–Stokes equations and does

∗ Corresponding author at: School of Engineering, RMIT University, PO Box 71, Bundoora, VIC 3083, Australia. E-mail address: [email protected] (J. Tu).

https://doi.org/10.1016/j.enbuild.2018.11.043 0378-7788/© 2018 Elsevier B.V. All rights reserved.

not request ultra-fine mesh to capture the features of turbulence eddies. According to a study by Hooff et al. [1] who simulated the cross-ventilation flows of an isolated building using the RANS and LES models, respectively, the computational demands of RANS is 80–100 times smaller than that of LES although the latter performs better in reproducing the transient features. The solution of the RANS equations requires the supplementation of additional terms to describe the nonlinear Reynolds stress term, which has led to the development of the Reynolds stress models and eddy viscosity models. The Reynolds stress model [2] relies on the solution of the Reynolds stress transport equations while the eddy viscosity models involve the modelling of the eddy viscosity or turbulence viscosity. The eddy viscosity models are the most prevailing turbulence models for indoor airflow simulations. Depending on the approach to model the turbulence viscosity, there are the zero-equation model [3], one-equation model [4], two-equation models (e.g., the k-ε and k-ω group models [5,6]) and hybrid models (e.g., the “two-layer” model [7] that combines a one-equation and k-ε models together), etc. The standard k-ε model was widely employed in the early CFD studies of indoor airflows [8], though it exhibited problems with treating natural convective heat transfer [9]. To improve its capability to capture the low-Reynolds number effects, the ReNormalization Group (RNG) k-ε model was developed [5]. The model uses an analytically derived differential equation to calculate

X. Li and J. Tu / Energy & Buildings 184 (2019) 8–18

9

Fig. 1. The approaches to model turbulence.

the effective viscosity in the near-wall region and gives better predictions of indoor airflows than the standard k-ε model. For a long period of time, the RNG k-ε model was regarded as the best turbulence model to simulate indoor airflows and a number of successful CFD simulations using the model have been reported [10–12]. Based on the benchmark experimental data of the International Energy Agency (IEA) Annex20 room, Yuce and Pulat [13] studied natural, forced and mixed convection, respectively. They stated that the RNG k-ε model achieved satisfactory predictions of the air velocity and temperature fields in all the cases. It is noticed that in order to achieve controllable experimental conditions, most of the validation experiments were conducted under isolated conditions using small-scale cavities and chambers [14,15] in which the velocity and temperature fields were strictly controlled and quite simple. However, the realistic indoor space is a complex environment involving various heat and mass transfer processes. The comprehensive CFD simulation involving multiple transfer processes in realistic indoor spaces has been found to be very sensitive to the selected turbulence model. For example, Liu et al. [16] simulated the air distribution in a MD-82 airliner cabin using the RNG k-ε model. The results revealed that the accuracy of the simulation was acceptable only when the cabin was unoccupied and the flow was quasi-isothermal. When the cabin was packed with heat-generating thermal manikins, the RNG k-ε model yielded significant predictive errors. However, a recent study by their co-workers [17] showed that the RNG k-ε model with a modified turbulent Schmidt number performed very well in predicting the contaminant dispersion in an airliner cabin mockup. It was proposed that the multi-transfer processes can significantly change the features of air turbulence and thus the turbulence model needs to be carefully formulated in order to capture their interactions. Unfortunately, the turbulence models have been rarely evaluated in terms of modelling multiple transfer processes and the selection of turbulence models is still on a case-by-case base. Convection and thermal radiation are two dominating mechanisms of heat transfer in common indoor spaces where the radiative mechanism generally accounts for 30–70% [18] of the total heat exchange. The study of Kogawa et al. [19] demonstrated that the surface radiation can enhance the velocity boundary layer in both vertical and horizontal directions and significantly alter the near-wall turbulence features. At the meantime, the gas radiation by the participating components such as H2 O and CO2 in the air can further enhance the air turbulence and thermal diffusion [20]. Apparently, the evaluation of turbulence models should be based on the comprehensive consideration that accounts for the effects of both convective and radiative heat transfer. However, the turbulence models have been rarely evaluated under the convectionradiation coupled indoor conditions as most of the previous studies simply neglected the radiative effects. It was until recent years that the importance of turbulence modelling in simulating convection-radiation coupled heat trans-

fer was realised. However, there are only a few studies [21–23] on this topic available in the literature. The common conclusion rising from these studies is that the SST k-ω model performs relatively better than the RNG k-ε model when both convective and radiative heat transfer exist. However, as mentioned before, the realistic indoor environment involves complex multi-transfer processes, which necessitates model validation in a more comprehensive circumstance. Among them, the dispersion of contaminants is an important consideration when assessing the indoor air quality. In particular, the emission of formaldehyde and volatile organic compounds (VOCs) by building materials [24] is a primary contaminant source in common built environment. As thermal radiation can considerably change the temperature of solid walls and turbulence features in the air [19,20], it will inevitably alter the emission and dispersion characteristics of indoor contaminants. It is proposed that an effective evaluation of the turbulence models should take into account the multi-transfer processes and their interactions. Thus, this study aims to evaluate turbulence models in terms of comprehensive modelling of convectionradiation coupled heat transfer in conjunction with contaminant emission and dispersion. A number of eddy viscosity turbulence models, including the zero-equation model, one-equation model, standard k-ε and k-ω models, SST model, RNG k-ε model and SST k-ω model, are evaluated using the experimental data reported in the literature. The models are also compared against each other using a realistic office room.

2. The theoretical models 2.1. The radiation model The indoor air flow field is modelled using the incompressible RANS equations coupled with a Boussinesq approximation for the thermal buoyance flows. The model equations have been fully validated in our previous studies [25–27] and will not be repeated here. Thermal radiation is assumed to happen only between solid surfaces as the concentration of participating components in common built environments is very low, so that the air is optically transparent and the absorption and emission by the air can be safely ignored [28]. In this case, thermal radiation only affects the airflow field through heating or cooling the solid surfaces. Surface radiation is modelled as a boundary condition of the energy equation. The Discrete Transfer model [29] is used owing to its proven accuracy and high cost-efficiency [30], in which the thermal radiation leaving a surface in a certain range of solid angle is approximated using a representative ray. The ray is then traced through the domain until it reaches another surface. To calculate the radiative heat exchange, all the surfaces are assumed to be grey and diffuse. For the kth (k = 1 - N) surface mesh element in the computational domain, the following energy balance correla-

10

X. Li and J. Tu / Energy & Buildings 184 (2019) 8–18

tion exists.





qext + qin,k = 1 − εr,k qin,k + εr,k σ

Tk4

+ hk (Tk − Ta )

(1)

where, qext is the externally applied heat flux, qin,k is the incident heat flux from the surroundings, ε r,k is the local emissivity, σ is the Stefan–Boltzmann constant, hk , Tk and Ta is the local convective heat transfer coefficient, surface temperature and air temperature, respectively. Eq. (1) means that the total energy exerted on the mesh element by external heating and radiation is dissipated through reflection, emission and convection. The local incident heat flux is calculated by integrating the radiation intensity approaching the element.



qin,k =

 



 

Iin,k  s · nd

(2)



s ·n

where, Ω is the hemispherical solid angle, Iin,k is the intensity of 



the coming ray, s and n are the ray direction vector and normal direction pointing out of the domain, respectively. For a diffuse surface, the local radiation intensity emitted from a surface mesh element is calculated by

Iout,k =

  εr,k σ Tk4 + 1 − εr,k qin,k π

(3)

2.2.3. The two-equation models In the two-equation models, the turbulence viscosity is modelled as the production of a turbulent velocity and length scale. Both the velocity and length scales are solved using separate transport equations. The turbulent velocity scale is computed from the turbulent kinetic energy while the turbulence length scale is estimated from two property parameters of the turbulence fields. Depending on the parameters selected to calculate the length scale, there are the k-ε and k-ω group models, respectively. The k-ε group models assume that the turbulence viscosity is correlated to the turbulent kinetic energy k and dissipation rate ε by

μt = Cμ ρ

k2

(9)

ε

where, Cμ is a constant. k and ε is modelled separately using their own transport equations, respectively. The k-ε group models have many variants depending on the approach of formulation. When the Navier–Stokes equations are analysed through the renormalization group [5], the transport equations for k and ε take the following forms (namely, the RNG k-ε model)

   ∂ ∂  ∂ μ ∂k μ+ t + Pk − ρε + PkB ρU j k = (ρ k ) + ∂t ∂xj ∂xj σk ∂ x j

2.2. The turbulence models

(10)

When turbulence modelling is based on the eddy viscosity hypothesis, the critical job is to formulate the eddy viscosity or turbulence viscosity μt , which is essential to the closure of the RANS equations through the effective viscosity

μe f f = μ + μt

(4)

Depending on the approach to formulate μt , different turbulence models have been developed, as briefed below. It should be noted that only the rationale of the turbulence models is briefly introduced. The details of the constants and coefficients can be found elsewhere [29] but are omitted here for concision. 2.2.1. The zero-equation model In the zero-equation model, no transport equation is solved to formulate μt and this is where the model has got its name. Instead, a constant turbulence viscosity is calculated as a production of a turbulent velocity scale Ut and a turbulence length scale lt .

μt = ρ fμUt lt

(5)

where, fμ is a proportionality constant. Ut is the maximum velocity in the fluid domain and lt is simply correlated to the domain volume VD by

lt =

 3

VD /7

(6)

Apparently, the zero-equation model has little physical foundation. 2.2.2. The one-equation model In the one-equation model [31], the turbulence viscosity is correlated to a dynamic kinematic eddy viscosity v˜ t by

μt = ρ v˜ t

   ∂ ∂  ∂ μt ∂ε ρU j ε = μ+ (ρε ) + ∂t ∂xj ∂xj σεRNG ∂ x j ε + (Cε1RNG Pk − Cε2RNG ρε + Cε1RNG PεB ) k

where, PkB is the buoyancy production term accounting for the buoyancy-induced turbulence.

PkB =

μt ∂T ρβ gi ρσρ ∂ xi









∂ρ v˜ t ∂ρU j v˜ t v˜ t ρ v˜ t ∂ v˜ + = c1 D1 ρ v˜ t S − c2 ρ S1e + μ + ∂t ∂xj LvK σt ∂ x j (8) where, σ t is a constant and LvK is the von Karman length scale [29].

(12)

Cε 1 RNG is the RNG coefficient defined by

Cε1RNG = 1.42 −

η (1 − η/4.38) 1 + η3 βRNG

(13)

The k-ω group models correlate the turbulence viscosity to the turbulence kinetic energy k and turbulence frequency ω via

μt = ρ

k

ω

(14)

Similar to the k-ε group models, the k-ω group models have a number of variants. To enhance the model performance in the near-wall regions, Menter [6] introduced a blending function F1 and developed the SST k-ω model based on the standard k-ω model of Wilcox [32]. The model consists of a k-ω model near the surface and a transformation of the k-ε model to a k-ω formulation in the outer region, where the transport equations of k and ω are, respectively

   ∂ ∂  ∂ μ ∂k μ+ t + Pk − β  ρ kω + PkB ρU j k = (ρ k ) + ∂t ∂xj ∂xj σk ∂ x j

(7)

The turbulent kinematic eddy viscosity is modelled using a transport equation.

(11)

(15)

   ∂ ∂  ∂ μ ∂ω μ+ t ρU j ω = (ρω ) + ∂t ∂xj ∂xj σω ∂ x j 1 ∂ k ∂ω ω + (1 − F1 )2ρ + α3 Pk − βρω2 + PωB σω 2 ω ∂ x j ∂ x j k

(16)

where, the blending factor F1 is equal to one near the surface and decreases to zero outside the boundary layer. In addition, the

X. Li and J. Tu / Energy & Buildings 184 (2019) 8–18

11

Fig. 2. The computational cases for model validation, (a) Case 1 (The experiments by Li et al. [33]), (b) Case 2 (The experiment by Xu et al. [34]). Table 1 The experimental conditions of Li et al. [33]. E (W)

ACR (h−1 )

Tin (°C)

300 1.0 16.0 Apparent HTC (W m−2 K−1 )

Tos (°C) Low-X

High-X

Low-Y

High-Y

Low-Z

High-Z

23.2 0.15

22.7 0.36

20.2 0.36

23.0 0.36

21.9 0.36

20.3 0.36

Note: E = applied heat power; ACR = air change rate; Tin = air temperature Tos = temperature of external wall surface, HTC = heat transfer coefficient.

transport of turbulent shear stress is restricted by a limiter to the formulation of the eddy viscosity.

μt a1 k = ρ max (a1 ω, SF2 )

(17)

Where F2 is a blending function restricting the limiter to the wall boundary layer.

at

inlet;

using the multi-component Eulerian equations [29]. Two separate continuity equations are solved for the air and CO2 , respectively.



 ∂ ( ρ m f a ) + ∇ · ρm f a U m = 0 ∂t

 ∂ (ρm fCO2 ) + ∇ · ρm fCO2U m − ρm Dk (∇ fCO2 ) = 0 ∂t

(18) (19)

3. Model validation and discussion

where, ρ m is the mixture density. fa and fCO 2 and the mass frac-

3.1. The validation cases

tion of air and CO2 , respectively (fa + fCO 2 = 1.0). U m is the massaveraged mixture velocity and Dk is the kinematic diffusivity of CO2 in the air.

The experimental data of Li et al. [33] and Xu et al. [34] are employed in this study for model validation. Li et al.’s experiments [33] (Case 1) were performed using a test chamber (4.2 m × 3.6 m × 2.75 m) containing a rectangular heat source (0.4 m × 0.3 m × 0.3 m), as illustrated in Fig. 2(a). Conditioned air (16 °C) was supplied through an inlet close to the floor and exhausted via a vent close to the ceiling, forming a displacement ventilation scheme. The airflow rate was 41.6 m3 /h, which was equivalent to an air change rate of 1.0 h−1 . The power of the heater was 300 W which was assumed to be equally distributed on the heater surface. The heater was made of aluminium (ε r = 0.1) and the interior walls of the chamber were painted black (ε r = 0.95). Details of the experimental conditions are shown in Table 1. Xu et al.’s [34] experiments (Case 2) were conducted in a chamber which had a similar ventilation scheme and dimensions (4.0 m × 3.88 m × 2.17 m) to those of Li et al. [33], as shown in Fig. 2(b). The chamber contained two cylindrical heaters (0.4 mdiameter × 1.06 m-height) each with a power of 100 W. On the top of each heater was a perforated diffuser, through which pure carbon dioxide (CO2 ) gas was released at a flow rate of 6.6 L/min. The interior chamber walls were coated with polyurethane foam in order reduce the heat leak from the ambience. According to the measured air temperatures at the inlet and outlet, respectively, the total heat leak was calculated to be 26 W, which was assumed to be equally distributed on the chamber floor, ceiling and walls. Details of the experimental conditions are shown in Table 2. In Xu et al.’s [34] experiments pure CO2 was injected into the chamber. Given the high CO2 concentration near the diffusers (up to 100%), the effects of CO2 on the air flow field are accounted for



3.2. Comparison and discussions The computational domains shown in Fig. 2 are discretised using unstructured tetrahedral mesh with refined inflation layers applied close to all the solid surfaces for better depiction of nearwall airflows. Mesh independence is achieved at 1.6 and 2.8 million mesh elements for the domain shown in Fig. 1(a) and (b), respectively. The transport equations are discretised using a finitevolume method and solved using the commercial CFD code AnsysCFX 18.2. Convergence is achieved within 80 0 0 iterations when the residual of the continuity equation drops down to 1 × 10−6 . For Case 1, the predicted air velocity vectors and eddy viscosity in the plane X = 0 are shows in Fig. 3. The figure shows that all the turbulence models predict a rising thermal plume above the heater. However, these models yield different eddy viscosity fields. The zero-equation model, which estimates the eddy viscosity based on the maximum velocity and domain volume, yields a constant eddy viscosity in the entire domain. Other models, including the one-equation model, predict a lower eddy viscosity in the near-wall regions and a larger eddy viscosity in the thermal plume where the air velocity is higher. Comparatively, the k-ω group models predict the lowest eddy viscosity. Particularly, the kω group predict that the eddy viscosity in the bulk regions is two to three magnitude orders lower than that predicted with the kε group models. In addition, the results also show that the eddy viscosity fields yielded from the standard k-ω and SST k-ω models are very close to each other (Fig. 3(e) and (f)). The magnitude of the eddy viscosity predicted by all the models is between 10−8

12

X. Li and J. Tu / Energy & Buildings 184 (2019) 8–18 Table 2 The experimental conditions of Xu et al. [34]. Air supply rate (m3 /h) Inlet temperature (°C) CO2 injection rate (L/min) CO2 injection temperature (°C)

162 20 6.6 20

Heat input to each heater (W) Heat leak (W/m2 ) Emissivity of heaters Emissivity of walls

100 −0.389 0.91 0.95

Fig. 3. The eddy viscosity contours predicted with different turbulence models (Case 1, Plane X = 0).

Fig. 5. The predicted thermal conditions on the heater surface (Li et al.’ case [33]). Fig. 4. The comparison of the predicted air temperature against the experimental data (Case 1).

and 10−2 Pa.s, which is a wide range compared to the dynamic viscosity of the air (1.81 × 10−5 Pa s). Fig. 3 shows that the eddy viscosity is significantly larger than or at least comparable to the dynamic viscosity in the vast majority of the domain. This will have a remarkable effect on the prediction of other parameters. The predicted air temperature profiles along a selected vertical line are shown in Fig. 4, where the line (Line 1) is located at the low-X and high-Y corner. Fig. 4 shows that only the zero-equation model significantly fail to predict the air temperature, while all the other models achieve acceptable agreement with the experimental data. In particular, the ceiling temperature is predicted to be lower than the air immediately below it and the floor is hotter than the air above it although minor errors are still observed in some local

regions. Among the models, it seems that the standard k-ε model has the lowest accuracy, the RNG k-ε model performs slightly better than the standard k-ε model and the k-ω group models perform best to predict the air temperature. It seems according to Fig. 4 that the accuracy of the oneequation model is comparable to that of the k-ω group models in terms of predicting the air temperature. However, the numerical results of the surface temperature and radiative heat flux tell another story. Fig. 5 shows the area-averaged convective heat transfer coefficient (“Convective HTC” in the figure) and fraction of the radiative heat flux (“Fraction of RHF” in the figure) at the heater surface. The results reveal that the zero-equation and one-equation model generate every large convective heat transfer coefficients (113.5 and 113.4 W/(m2 K), respectively) while the standard k-ε and RNG k-ε models predict very low convective heat transfer coefficients (3.6 and 3.3 W/(m2 K), respectively). The standard k-ω

X. Li and J. Tu / Energy & Buildings 184 (2019) 8–18

13

Fig. 6. The comparison of the numerical results against the experimental data, (a) temperature increase (T - Tin ), (b) CO2 concentration C∗ = ((C - Cin )/(Cout - Cin )).

and SST k-ω models produce medium heat transfer coefficients of 29.1 - and 29.2 W/(m2 K), respectively. The lower convective heat transfer coefficient means less convective heat flux at the heater surface, which result in an elevated fraction of the radiative heat flux. As shown in the Fig. 5, the zero-equation model predicts that the radiative mechanism only accounts for 2.2% of the total heat

dissipation while the RNG k-ε predicts that the fraction is 29.4%. Strangely, the one-equation model generates a proportion of 14.6% which is comparable to the predictions of the k-ω group model when it predicts a very high convective heat transfer coefficient that is close to that of the zero-equation model. This prediction is strange and cannot be properly explained.

14

X. Li and J. Tu / Energy & Buildings 184 (2019) 8–18

Fig. 7. The predicted thermal flow field and CO2 concentration field in the X = 0 Plane.

Fig. 8. The thermal flow field predicted with the pure-convection assumption, (a) RNG k-ε , (b) SST k-ω.

However, the actual convective heat transfer coefficient was not measured in Li et al.’s experiments [33], which is estimated here using the empirical correlation [35].

Nu = hl /λ =c (Gr · P r )n

(20)

where, Nu, Gr and Pr are dimensionless Nusselt number, Grashof number and Prandtl number, respectively. h, l and λ are the convective heat transfer coefficient, characteristic length and thermal conductivity, respectively. c and n are empirical constants recommended by Lienhard (c = 0.59 and n = 1/3) [35]. Based on the CFD results, it is calculated that the area-averaged convective heat transfer coefficient of the rectangular heater is between 13.5 and 27.8 W/(m2 K), which is close to the values predicted with the k-ω group models (29.1 and 29.2 W/(m2 K)). Therefore, it is proposed that the k-ω group models perform best when describing the convective heat transfer on the heater surfaces. The distinctly different predictive capabilities of the k-ε group and k-ω group models in terms of predicting the convectionradiation heat flux split on the heater surface are also found in Case 2. Fig. 6 shows the comparison between the predicted air temperature and CO2 concentration profiles along 6 selected vertical lines (see Fig. 5 for the exact locations). The zero-equation and one-equation models are omitted here due to their severe over-predictions of the convective heat transfer. Fig. 6(a) shows the both groups of two-equation turbulence models have acceptable accuracy to predict the air temperature, however, Fig. 6(b) indicates that the k-ω group models perform clearly better than the k-ε group models to predict the CO2 concentration. The k-ω group models predict that the region with higher CO2 concentration is

located above the height of the CO2 diffusers (Z = 1.06 m), which is in a good agreement with the experimental measurements. However, the k-ε group models fail to capture the CO2 concentration profiles and incorrectly predict the high-concentration region is located lower than the diffuser height. In addition, Fig. 6 once again suggests that the standard k-ω and SST k-ω models have very close accuracies in predicting the temperature and concentration profiles. The predicted thermal flow and CO2 concentration fields in the plane X = 0 are illustrated in Fig. 7, which can help understand the deviations shown in Fig. 6(b). Here, only the CFD results of the RNG k-ε and SST k-ω model are shown as an illustration because those of the standard k-ε and RNG k-ε models are similar and those of the standard k-ω and SST k-ω models are very close. The results show that when the RNG k-ε model is used, no visible thermal plume is predicted above the heaters. Instead, remarkable descending flows of CO2 following the heater surface is observed (Fig. 7(a)). Consequently, the high-CO2 -concentration regions are predicted to be in the medium-to-lower levels around the heaters (Fig. 7(c)). On the contrary, the SST k-ω model predicts significant rising thermal plumes (Fig. 7(b)), which form “air cages” that trap the CO2 gas in a region above the heaters. As a result, the highconcentration regions are predicted to be located above the heaters (Fig. 7(d)). Similar to Case 1, the computations of Case 2 reveal that the k-ε group models predict much lower convective heat transfer coefficient and larger fraction of radiative heat flux at the heater surfaces than the k-ω group models. The area-averaged convective heat transfer coefficient predicted by the standard k-ε and RNG k-ε

X. Li and J. Tu / Energy & Buildings 184 (2019) 8–18

15

Fig. 9. The office room model.

Fig. 10. The thermal flow fields and contaminant dispersion predicted with different turbulence models.

models are 2.3 and 2.0 W/(m2 K), respectively, while the standard k-ω and SST k-ω model predict it to be one magnitude-order larger (26.3 and 26.4 W/(m2 K), respectively). Correspondingly, the fraction of the radiative heat flux yielded from the turbulence models is 51.1% (standard k-ε ), 53.% (RNG k-ε ), 41.9% (standard k-ω) and 42.2% (SST k-ω), respectively. It is suggested that the descending gas flows predicted by the k-ε group models (Fig. 6(a) and (c)) are due to the insufficient accounting for the convective heat transfer as buoyancy flow is dominantly driven by the convective heat transfer between the heater surface and the adjacent air. On the other hand, many studies of the turbulence models reported in the literature were conducted under the “pure convec-

tion” assumption [13,16]. In this case, thermal radiation is ignored and all the heat released by the heater is assumed to be solely transferred by the convective mechanism. The computations of this study reveal that when thermal radiation is ignored, the both models generate similar air temperature fields. As shown in Fig. 8, when the pure-convection assumption was used, both the RNG k-ε and SST k-ω models predict clear thermal plumes above the heaters, which are distinctly different from the descending airflows shown in Fig. 7. As the thermal plumes are solely driven by convective heat transfer, it confirms that the RNG k-ε model seriously underestimates the convective heat transfer when thermal radiation co-exists.

16

X. Li and J. Tu / Energy & Buildings 184 (2019) 8–18

Fig. 11. The surface temperature predicted with different turbulence models, (a) RNG k-ε , (b) SST k-ω.

4. Case study

Fig. 12. The area-averaged surface temperature of selected surfaces.

Although the experimentally visualised airflow field in the chamber is unavailable to fully validate the turbulence models, numerous experimental observations have proven the existence of thermal plumes above heating objects in low-momentum indoor spaces. In particular, Lincina et al. [36] experimentally studied the micro-environment of a number of thermal manikins with the heat input ranging from 49 W to 89 W and observed clear thermal plumes above all the manikins. The experimental observations by Zukowska et al. [37] further proved that the exhalation jets of 6.0 L/min from the mouth and nose does not significantly affect the thermal plume. In the experiments of Xu et al. [34], the heat input of each heater was 100 W and the CO2 flow rate was 6.6 L/min, which were comparable to the experimental conditions of Lincina et al. [36] and Zukowska et al. [37]. Therefore, it is reasonable to expect that obvious thermal plumes should exist above the heaters. In addition, considering the better agreement between the measured CO2 concentration and the predictive results (Fig. 6(b)), it can be safely concluded that the SST k-ω model performs better. Some existing studies suggested that the SST k-ω model performs slightly better than the RNG k-ε models [21,38] to simulate the convection-radiation coupled heat transfer. However, it should be noted that the comparisons were limited to the air temperature, just analogous to that shown in Figs. 4 and 6(a). In this case, the convection-radiation split and near-wall airflow patterns were not examined. On the contrary, this study suggests that the advantage of the k-ω group models could be significant if one looks at the heater surfaces.

It is well known that the room geometry significantly affects the distribution of radiative energy and airflow pattern [39]. It is hence necessary to evaluate the models under realistic conditions. Therefore, a typical office room is designed, as shown in Fig. 9. The room is equipped with a mixing ventilation system where both the supply and exhaust vents are located close to or in the ceiling. Two manikins representing the office occupants (downloaded at http://www.ie.dtu.dk/manikin [40]) are sitting in the room. The heat load of each manikin is 80 W, which represents the average human metabolic heat exclusive of the loss through breathing and sweating [41]. Other heat sources in the room include two desktop personal computers (PCs, 100 W each) and two lights (80 W each). A summer season scenario is assumed with an external atmospheric temperature of 37 °C. Heat invasion from the ambience is assumed to occur across the low-Y and high-Y walls which have an equivalent heat transfer coefficient of 1.5 W/m2 [42]. The other surfaces, including the floor, ceiling, desks, shelf and low-X and high-X walls, are all assumed to be adiabatic. The ventilation parameters are carefully calculated according to the ventilation standard ASHRAE 62.1-2007 [42] with a target indoor temperature of 24 °C, which results in an air supply rate of 492 m3 /h (equivalent to an air exchange rate of 12.8 h−1 ) and inlet temperature of 18 °C. In order to evaluate the effects of turbulence modelling on the prediction of contaminant release and dispersion, it is assumed that Manikin_2 is releasing some kind of gaseous contaminant (e.g., odour) and the shelf is releasing formaldehyde. Among the two selected contaminant sources, one is a heat source while the other is an adiabatic surface. For common indoor gaseous contaminants, their concentrations are very low so that the existence of the contaminants has not detectable influence on the air flow field. Therefore, the dispersion of the contaminants is modelled using a transportable scalar [29].



∂ (ρC ) + ∇ · ρCU −  ∇C = 0 ∂t

(21)

where, C is a scalar representing the contaminant concentration,  is the turbulent diffusivity of C. The same numerical procedures as aforementioned (e.g, the unstructured tetrahedral mesh with inflation layers and finite volume method, etc) are used for the discretisation and model solution. Emphasis here is put on the comparison between the RNG k-ε and SST k-ω models according to the outcomes of the above model validation. The predicted thermal flow fields in the plane cutting through Manikin_2 (Y = 2.3 m) are shown in Fig. 10. The results show that the air temperature fields predicted by the both models are very

X. Li and J. Tu / Energy & Buildings 184 (2019) 8–18

17

Fig. 13. The volumetric distribution of formaldehyde concentration (C∗ = C / Cmax ), (a) RNG k-ε , (b) SST k-ω.

close to each other. However, due to the underestimation of convective heat transfer by the RNG k-ε model, it fails to predict the thermal plumes above the manikin head (Fig. 10(a)) while the SST k-ω model succeeds to do so (Fig. 10(b)). Apparently, the different airflow fields in the vicinity of the manikins will inevitably lead to different predictions of the dispersion and distribution of contaminants if the contaminants are released by the manikin. As shown in Fig. 10(c) and (d), the RNG k-ε and SST k-ω models predict different odour distribution pattern in the air, despite the odour concentration at the surface of Manikin_2 is a constant. The SST k-ω model predicted a faster contaminant dispersion in the air, which is due to the thermal plume resulted from the stronger convective heat transfer predicted at the manikin surface. To evaluate the effects of turbulence modelling on the prediction of the convection-radiation interactions, the temperature of the solid surfaces and their area-averaged values are shown in Figs. 11 and 12, respectively. Similar to the validation cases, the RNG k-ε model predicts a higher surface temperature than the SST k-ω model on all the heat-releasing and heat-gaining surfaces. In particular, the temperature of the high-Y wall, which is close to the heat-releasing PCs, is predicted to be significantly higher by the RNG k-ε model than by the SST k-ω model. In addition, the RNG k-ε model generates a higher temperature at the manikin surfaces, which will certainly result in different evaluation of thermal comfort. The higher surface temperature of the heat sources predicted by the RNG k-ε model is attributed to the underestimation of convective heat transfer, which ultimately results in the larger fraction of radiative heat flux at the heat-releasing surfaces and higher temperature at the radiation-gaining surfaces, as shown in Fig. 12. However, it seems according to Fig. 12 that the relative error between the RNG k-ε and SST k-ω models is sensitive to the airflow conditions near the surface which are jointly shaped by the room layout and ventilation setup. For example, the temperature of Mankin_1 who sits in a low-velocity region is not that sensitive to the selected turbulence model while the thermal conditions of Mankin_2 are strongly affected by the turbulence model. As the turbulence model has a strong effect on the prediction of radiative heat transfer and surface temperature, it is expected it can affect the evaluation of contaminant omission by building materials as it is strongly correlated to the material temperature. In this study, the shelf is assumed to be newly painted and omitting formaldehyde. The formaldehyde concentration at the shelf surface is linearly correlated to the surface temperature in the range of 20– 27 °C according to the experimental measurements of Xiong et al [24]. The predicted 3-dimensional distributions of the formaldehyde concentration in the room are shown in Fig. 13, where the concentration is normalised in terms of the maximum value at the shelf surface. As shown in the figure, the RNG k-ε model predicts higher formaldehyde concentrations than the SST k-ω model, both at the

shelf surface and in the bulk air. However, the distribution pattern of formaldehyde predicted by the models are significantly different. The RNG k-ε model yields an evenly distributed formaldehyde concentration in the entire room while the SST k-ω model predicts a remarkable “formaldehyde plume” above the shelf in which the concentration is 3–4 times higher than that predicted with the RNG k-ε model. The figures once again confirm the inefficiency of the RNG k-ε model to model the near-wall and thermal buoyancy flows. Although the RNG k-ε model predicts a higher temperature and formaldehyde concentration at the shelf surface, it fails to depict the thermal flow and convective heat transfer in the near-wall regions. Consequently, it cannot capture the rising airflow near the vertical walls, which ultimately results in an ineffective prediction of contaminant dispersion. 5. Conclusions The viscosity eddy turbulence models including the zeroequation, one-equation, standard k-ε , RNG k-ε , standard k-ω and SST k-ω models are evaluated in terms of modelling the convection-radiation coupled heat transfer in conjunction with contaminant dispersion in indoor environment. The CFD results yielded from the models are compared against the experimental measurements of air temperature and contaminant concentration reported in the literature. If one only looks at the air temperature, all the models except the zero-equation model have acceptable accuracy and the one-equation model even looks performing better than the k-ε group models. However, a quantitative analysis of the convective and radiative heat fluxes at the heater surface suggests that the models performs distinctly different. A magnitudeorder analysis using the empirical correlation demonstrates that the zero-equation and one-equation models severely overpredict while the k-ε group models seriously underestimate the convective heat transfer coefficient. On the other hand, the k-ω group models achieve satisfactory agreement with both the experimental data and empirical calculations. The underestimation of the convective heat transfer by the RNG k-ε model leads to the inadequate calculation of the near-wall thermal buoyancy flows, which ultimately results in the incorrect predictions of the dispersion path and distribution pattern of contaminants in the air. Among the models studied, the standard k-ω and SST k-ω models performs best in modelling the convection-radiation coupled heat transfer. In should be noted that thermal radiation accounts for 30–70% [18] of the total heat exchange in common indoor spaces. The incident radiative energy onto solid surfaces is transferred to the air through convection, which can significantly change the airflow patterns in the near-wall regions although its effect on the temperature of the bulk air is found to be insignificant. Accurate description of the near-wall airflows is vital to the evaluations of air distribution, thermal comfort and contaminant dispersion. However,

18

X. Li and J. Tu / Energy & Buildings 184 (2019) 8–18

ignoring the radiative mechanism and numerically enforcing all the heat was solely transferred by the convective mechanism, just as many existing studies have assumed [13,16], could result in significant predictive errors of the near-wall airflows. In fact, many contaminants in indoor environment are released from solid materials and the omission is strongly related to the material temperature, therefore, a proper turbulence model that is able to effectively calculate the radiation-convection heat flux split is vital to the effective evaluation of the overall thermal comfort and indoor air quality. This study demonstrates that for the comprehensive CFD simulations that involve both convective and radiative heat transfer, the SST k-ω model performs significantly better than the widely accepted RNG k-ε model. Conflict of interest statement Xiangdong Li and Jiyuan Tu declare that: we have no proprietary, financial, professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript titled: Evaluation of the eddy viscosity turbulence models for the simulation of convection-radiation coupled heat transfer in indoor environment. Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.enbuild.2018.11.043. References [1] T. van Hooff, B. Blocken, Y. Tominaga, On the accuracy of CFD simulations of cross-ventilation flows for a generic isolated building: comparison of RANS, LES and experiments, Build. Environ. 114 (2017) 148–165. [2] M.F. King, C.J. Noakes, P.A. Sleigh, M.A. Camargo-Valero, Bioaerosol deposition in single and two-bed hospital rooms: a numerical and experimental study, Build. Environ. 59 (2013) 436–447. [3] Q. Chen, W. Xu, A zero-equation turbulence model for indoor airflow simulation, Energy Build. 28 (2) (1998) 137–144. [4] J. Taghinia, M.M. Rahman, T. Siikonen, Numerical simulation of airflow and temperature fields around an occupant in indoor environment, Energy Build. 104 (2015) 199–207. [5] V. Yakhot, S.A. Orszag, S. Thangam, T.B. Gatski, C.G. Speziale, Development of turbulence models for shear flows by a double expansion technique, Phys. Fluids A 4 (7) (1992) 1510–1520. [6] F.R. Menter, Two-equation eddy-viscosity turbulence models for engineering applications, AIAA J. 32 (8) (1994) 1598–1605. [7] W. Xu, Q. Chen, A two-layer turbulence model for simulating indoor airflow: part I. Model development, Energy Build. 33 (6) (2001) 613–625. [8] Q. Chen, J. Van Der Kooi, A methodology for indoor airflow computations and energy analysis for a displacement ventilation system, Energy Build. 14 (4) (1990) 259–271. [9] D. Risberg, M. Risberg, L. Westerlund, CFD modelling of radiators in buildings with user-defined wall functions, Appl. Therm. Eng. 94 (2016) 266–273. [10] P. Rohdin, B. Moshfegh, Numerical modelling of industrial indoor environments: a comparison between different turbulence models and supply systems supported by field measurements, Build. Environ. 46 (11) (2011) 2365–2374. [11] G.M. Stavrakakis, M.K. Koukou, M.G. Vrachopoulos, N.C. Markatos, Natural cross-ventilation in buildings: building-scale experiments, numerical simulation and thermal comfort evaluation, Energy Build. 40 (9) (2008) 1666–1681. [12] J.D. Posner, C.R. Buchanan, D. Dunn-Rankin, Measurement and prediction of indoor air flow in a model room, Energy Build. 35 (5) (2003) 515–526. [13] B.E. Yuce, E. Pulat, Forced, natural and mixed convection benchmark studies for indoor thermal environments, Int. Commun. Heat Mass Transf. 92 (2018) 1–14. [14] P.F. Linden, G.F. Lane-Serff, D.A. Smeed, Emptying filling boxes: the fluid mechanics of natural ventilation, J. Fluid Mech. 212 (2006) 309–335. [15] M.A. Menchaca-Brandan, F.A.D. Espinosa, L.R. Glicksman, The influence of radiation heat transfer on the prediction of air flows in rooms under natural ventilation, Energy Build. 138 (2017) 530–538.

[16] W. Liu, J. Wen, C.-H. Lin, J. Liu, Z. Long, Q. Chen, Evaluation of various categories of turbulence models for predicting air distribution in an airliner cabin, Build. Environ. 65 (2013) 118–131. [17] F. Li, J. Liu, J. Ren, X. Cao, Predicting contaminant dispersion using modified turbulent Schmidt numbers from different vortex structures, Build. Environ. 130 (2018) 120–127. [18] M.S. Owen, ASHRAE Handbook of Fundamentals, 30329, American Society of Heating, Refrigeration and Air-Conditioning Engineers, Atlanta, GA, USA, 2017. [19] T. Kogawa, J. Okajima, A. Sakurai, A. Komiya, S. Maruyama, Influence of radiation effect on turbulent natural convection in cubic cavity at normal temperature atmospheric gas, Int. J. Heat Mass Transf. 104 (2017) 456–466. [20] T. Kogawa, L. Chen, J. Okajima, A. Sakurai, A. Komiya, S. Maruyama, Effects of concentration of participating media on turbulent natural convection in cubic cavity, Appl. Therm. Eng. 131 (2018) 141–149. [21] S. Hussain, P.H. Oosthuizen, A. Kalendar, Evaluation of various turbulence models for the prediction of the airflow and temperature distributions in atria, Energy Build. 48 (2012) 18–28. [22] S. Gilani, H. Montazeri, B. Blocken, CFD simulation of stratified indoor environment in displacement ventilation: validation and sensitivity analysis, Build. Environ. 95 (2016) 299–313. [23] T. Kobayashi, K. Sugita, N. Umemiya, T. Kishimoto, M. Sandberg, Numerical investigation and accuracy verification of indoor environment for an impinging jet ventilated room using computational fluid dynamics, Build. Environ. 115 (2017) 251–268. [24] J. Xiong, P. Zhang, S. Huang, Y. Zhang, Comprehensive influence of environmental factors on the emission rate of formaldehyde and VOCs in building materials: correlation development and exposure assessment, Environ. Res. 151 (2016) 734–741. [25] X. Li, Y. Yan, Y. Shang, J. Tu, An Eulerian—Eulerian model for particulate matter transport in indoor spaces, Build. Environ. 86 (2015) 191–202. [26] X.D. Li, K. Inthavong, J.Y. Tu, Particle inhalation and deposition in a human nasal cavity from the external surrounding environment, Build. Environ. 47 (2012) 32–39. [27] X. Li, Y. Yan, J. Tu, The simplification of computer simulated persons (CSPs) in CFD models of occupied indoor spaces, Build. Environ. 93 (Part 2) (2015) 155–164. [28] G. El Hitti, M. Nemer, K. El Khoury, Reducing CPU time for radiative exchange and transient heat transfer analysis using zone method, Numer. Heat Transf. Part B 56 (1) (2009) 23–37. [29] Ansys Inc., Ansys CFX-solver theory guide (Release 15.0), Southpointe, Canonsburg, PA, 2013 15317. [30] X. Li, Y. Yan, J. Tu, Evaluation of models and methods to simulate thermal radiation in indoor spaces, Build. Environ. 144 (2018) 259–267. [31] F.R. Menter, Eddy viscosity transport equations and their relation to the k-epsilon model, J. Fluid Eng. – Trans. ASME 119 (4) (1997) 876–884. [32] D.C. Wilcox, Multiscale model for turbulent flows, AIAA J. 26 (11) (1988) 1311–1320. [33] Y. Li, M. Sandberg, L. Fuchs, Vertical temperature profiles in rooms ventilated by displacement: full-scale measurement and Nodal modelling, Indoor Air 2 (1992) 225–243. [34] M. Xu, T. Yamanaka, H. Kotani, Vertical profiles of temperature and contaminant concentration in rooms ventilated by displacement with heat loss through room envelopes, Indoor Air 11 (2) (2001) 111–119. [35] J.H.I. Lienhard, J.H.V. Lienhard, A Heat Transfer Textbook, fourth ed., Phlogiston Press, Cambridge, MA, USA, 2018. [36] D. Licina, J. Pantelic, A. Melikov, C. Sekhar, K.W. Tham, Experimental investigation of the human convective boundary layer in a quiescent indoor environment, Build. Environ. 75 (2014) 79–91. [37] D. Zukowska, A. Melikov, Z. Popiolek, Impact of personal factors and furniture arrangement on the thermal plume above a sitting occupant, Build. Environ. 49 (2012) 104–116. [38] T. Wu, C.W. Lei, On numerical modelling of conjugate turbulent natural convection and radiation in a differentially heated cavity, Int. J. Heat Mass Transf. 91 (2015) 454–466. [39] R. Zhuang, X. Li, J. Tu, CFD study of the effects of furniture layout on indoor air quality under typical office ventilation schemes, Build. Simul. 7 (3) (2014) 263–275. [40] D.N. Sørensen, L.K. Voigt, Modelling flow and heat transfer around a seated human body by computational fluid dynamics, Build. Environ. 38 (6) (2003) 753–762. [41] M. Luo, Z. Wang, K. Ke, B. Cao, Y. Zhai, X. Zhou, Human metabolic rate and thermal comfort in buildings: the problem and challenge, Build. Environ. 131 (2018) 44–52. [42] A.S. Committee, Ventilation for acceptable indoor air quality, American Society of Heating, Refrigerating and Air-Conditioning Engineers, Inc, 1791 Tullie Circle NE, Atlanta, GA 30329, 2007.