Numerical Experiments in

Turbulent Natural Convection

Using Two-Equation
X. Albets-Chico
A. Oliva1
Eddy-Viscosity Models
This work is focused on the simulation and prediction of turbulent natural convection
C. D. Pérez-Segarra flows by means of two-equation eddy-viscosity models. In order to show the generality,
precision, and numerical issues related to these models under natural convection, three
Centre Tecnològic de Transferència de Calor,
different buoyancy-driven cavities have been simulated: a tall cavity with a 30:1 aspect
Universitat Politècnica de Catalunya (UPC),
ratio, a cavity with a 5:1 aspect ratio, and, finally, a 4:1 aspect ratio cavity. All cases are
ETSEIAT, C/ Colom 11,
solved under moderate and/or transitional Rayleigh numbers (2.43 ⫻ 1010, 5 ⫻ 1010, and
08222 Terrassa, Catalunya, Spain
1 ⫻ 1010, respectively) and all simulations are compared to experimental and/or direct
numerical simulation data available in literature. These different situations allow to
check the applicability of two-equation eddy-viscosity models in buoyancy-driven flows,
giving criteria on computational effort/precision and their physical behavior.
关DOI: 10.1115/1.2907432兴

Keywords: turbulence, buoyancy, natural convection, two equation, k-⑀, k-␻

1 Introduction bulent buoyancy-driven differentially heated cavities with differ-

ent aspect ratios 共30:1, 5:1, and 4:1兲 at different Rayleigh num-
Turbulent flows, driven or affected by buoyancy, are an impor-
bers, as shown henceforth.
tant type of flows that arise in many engineering applications,
The first cavity 共also called tall cavity because of its high
such as heating, ventilating, and air conditioning 共HVAC兲, venti-
height/width relation兲 presents a characteristical structure. This
lated facades, electronic equipment, solar energy collectors, cool-
flow is completely developed, i.e., hot and cold thermal boundary
ing of nuclear reactors, material processing, etc. Heat transfer and layers interact and show a turbulent core and laminar sublayers
fluid flow related to these cases are strongly affected by the tran- near the walls. The influence of the corners is negligible, and the
sition point from laminar to turbulent flow and by the Rayleigh temperature field is strongly affected by convection. The second
number. case exhibits a different behavior. Because of the lower aspect
Direct CFD solution of the Navier–Stokes equations demands ratio, two different thermal boundary layers are developed, and
high CPU-time requirements due to the unsteadiness and three- the lower part of the hot wall and the upper of the cold, e.g., the
dimensional nature of turbulent situations. Consequently, direct core of the cavity, remain laminar. In this situation, the ability of
numerical simulations 共DNS兲 关1兴 are only available, at the mo- different models to describe the transition point will be an impor-
ment, in low-Rayleigh number and simple geometry cases. For tant aspect to achieve credible predictions. This aspect is also
more general situations, it is necessary to use other kinds of tools related to different problems, which affect the convergence and
to obtain numerical solutions in turbulent flows with a much lower asymptotic behavior of the results, as has been indicated in litera-
CPU-time effort. ture 关2,8兴. Finally, the third cavity, because of the even smaller
Among other types of models 关1兴, Reynolds averaged Navier– aspect ratio and the lower Rayleigh number, amplifies previous
Stokes 共RANS兲 techniques have been chosen in this work. RANS perspective and brings clarity to the limiting skills and drawbacks
techniques are based on a time-averaged filtering of the governing of RANS models to predict turbulent natural convection.
equations. Historically, these models have been designed and fit-
ted under turbulent forced convection situations. Because of this,
different researchers 共Davidson 关2,3兴, Hanjalić 关4–6兴, Ince and 2 Mathematical Models
Launder 关7兴, etc.兲 have carried out an additional effort to enable Different two-equation eddy-viscosity models have been tested
these models to solve turbulent natural convection cases with a to predict fluid flow and heat transfer in turbulent natural convec-
relative precision and low computational cost. tion cavities. As it is well known, such kind of models are de-
Since there are a lot of industrial situations where information signed over a time-averaged Navier–Stokes set of equations
on turbulent natural convection is required, the main objective of 共RANS兲. The models use the “eddy” viscosity as an intermediate
this work is to show the possibilities of RANS two-equation eddy- variable to model turbulent flow stresses and heat transfer. This
viscosity models as a tool to describe them. To achieve this task in “intermediate” variable is obtained using two different turbulent
a moderately high-Rayleigh number situations, different k-⑀ and quantities: turbulent kinetic energy and a variable related to its
k-␻ models have been used. dissipation. Both variables are obtained from the resolution of
These models allow to obtain numerical solutions in a relative their own transport equation.
short CPU time, with reasonably accurate results on different tur-
2.1 Governing Equations. The time-averaged governing
equations 共continuity, momentum and energy兲 of the fluid flow
1
Corresponding author. assuming fluid Newtonian behavior, constant thermophysical
Contributed by the Heat Transfer Division of ASME for publication in the JOUR-
NAL OF HEAT TRANSFER. Manuscript received March 5, 2007; final manuscript re-
properties, nonparticipant radiant medium and negligible heat fric-
ceived October 17, 2007; published online May 19, 2008. Review conducted by tion, and influence of pressure on temperature are written in tensor
Gautam Biswas. notation as

Journal of Heat Transfer Copyright © 2008 by ASME JULY 2008, Vol. 130 / 072501-1

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/28/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 Table 1 Terms in turbulence models used in this work

Term IL model GPC model WX88 model WX93 model PDH+ D model

␲k ␳共˜⑀ + D兲 ␳⑀ C␮␳␻k C␮ f k␳␻k C␮ f k␳␻k

Sk Gk Gk Gk Gk Gk
Pz c1 Pk共˜⑀ / k兲 c1 f 1 Pk共⑀ / k兲 c1 Pk共␻ / k兲 c1 f 1 Pk共␻ / k兲 c1 f 1 Pk共␻ / k兲
␲z c2 f 2␳共˜⑀2 / k兲 c2 f 1␳共⑀2 / k兲 c2␳␻2 c2␳␻2 c2␳␻2
Sz c1c3Gk共˜⑀ / k兲 + E + YAP c1c3 f 1Gk共⑀ / k兲 + f 1YAP — — C␻

⳵ūi tional to both the square root of the turbulent kinetic energy 共k1/2兲
=0 共1兲 and the turbulent length scale 共l兲, which is given by l = k3/2 / ⑀ for
⳵xi
k-⑀ models and by l = k1/2 / ␻ for k-␻ models.
⳵ūi ⳵ūi ⳵ p̄ ⳵␶ij The high- and low-Reynolds number k-␻ models by Wilcox
␳ + ␳ū j =− + − ␳␤共T̄ − T0兲gi 共2兲 共WX88 关9兴, WX93 关10兴兲 and by Peng–Davidson–Holmberg
⳵t ⳵x j ⳵xi ⳵x j 共PDH+ D兲 with Gk corrections 关2兴 have been used in this work.
Among k-⑀ models, the Ince–Launder 共IL兲 model 关7兴 and the
⳵T̄ ⳵T̄ 1 ⳵q̇i Goldberg–Peroomian–Chakravarthy 共GPC兲 model 关11兴 have been
␳ + ␳ūi =− 共3兲
⳵t ⳵xi c p ⳵xi employed. A summary of the terms used by these models is pre-
sented in Table 1.
where
The generic dissipative variable 共z兲 adopts three different forms

␶ij = ␮ 冉 ⳵ūi ⳵ū j

+
⳵x j ⳵xi
冊
− ␳ui⬘u⬘j and qi = − ␭
⳵T̄
⳵xi
+ c p␳ui⬘T⬘
depending on the model. GPC model 关11兴 uses the standard dis-
sipation rate 共⑀兲, while the IL model 关7兴 uses the “corrected” dis-
sipation rate 共˜⑀兲. On the other hand, WX88 关9兴, WX93 关10兴, and
As previously mentioned, as a consequence of the time-averaging PDH+ D 关2兴 models employ the specific dissipation rate 共␻兲. The
process of the Navier–Stokes equations, new terms have appeared corrected dissipation rate 共˜⑀兲 is a numerical correction needed to
in momentum and energy equations. They are the Reynolds stress achieve a zero dissipation value at the wall 共˜⑀w = 0兲. Consequently,
tensor 共␳ui⬘u⬘j 兲 and turbulent heat flux 共␳ui⬘T⬘兲. Using eddy- a D term in k-equation has to be added, which is the asymptotic
viscosity models, these terms are modeled by analogy with the value of ⑀ at wall: D = 2␯⳵k1/2 / ⳵x j. Hence, for the IL model ⑀ at
Stokes viscosity law and the Fourier heat conduction through the wall is set to 0 while for the GPC model is set to ⑀w = 2␯k p / y 2p
eddy viscosity 共␮t兲: 共being k p the turbulent kinetic energy of the first inner node and y p

␳ui⬘u⬘j = − ␮t 冉 ⳵ui ⳵u j
+
⳵x j ⳵xi
冊
2
+ ␳k␦ij
3
共4兲
the distance to the wall兲. Finally, for the k-␻ models, ␻ = 6␯ / ␤y 2p
has been prescribed at the first inner node, as the value at the wall
is theoretically infinity 共being y p the distance to the wall兲.
Models constants and damping functions have not been here
␮t ⳵T̄ included for brevity. They are clearly stated in original references:
␳ui⬘T⬘ = − 共5兲
␴T ⳵xi Ref. 关11兴 for the GPC model, Ref. 关7兴 for the IL model, Ref. 关9兴
for the WX88 model, Ref. 关10兴 for the WX93 model, and, finally,
2.2 Turbulence Models. Two-equation eddy-viscosity mod- Ref. 关2兴 for the PDH+ D model.
els are based on the resolution of the turbulent viscosity 共␮t兲 by 2.3 Buoyancy Effects. The treatment of the buoyancy term in
means of two different turbulent quantities: the turbulent kinetic the turbulent transport equations needs to be especially observed,
energy 共k兲 and another quantity related to its dissipation 共e.g., the since the behavior of turbulence models in natural convection is
dissipation rate of turbulent kinetic energy 共⑀兲, the specific turbu- strongly affected by this term. Table 2 indicates the different Gk
lent kinetic energy dissipation 共␻兲, the characteristic time of tur- terms adopted for the models. Simple gradient diffusion hypoth-
bulent kinetic energy dissipation 共␶兲, etc.兲. In this paper, k-⑀ and esis 共SGDH兲 has been adopted for GPC, WX88, and WX93 mod-
k-␻ models are employed. Both turbulent quantities are obtained els in this work. Since these models were designed under forced
from the resolution of their own transport equations: Eq. 共6兲 for convection cases, present adaptation to natural convection flows
the turbulent kinetic energy 共k兲 and Eq. 共7兲 for its dissipation rate has been carried out using the simplest method. On the contrary,
共z = ⑀ or z = ␻兲. Hence, the additional set of equations that closes IL and PDH+ D models show specific treatments for buoyancy
the problem is term, while IL model uses the generalized gradient diffusion hy-

␳
⳵k
⳵t
+ ␳ūi
⳵k
=
⳵
⳵xi ⳵xi
冋冉 冊 册␮+
␮␶ ⳵k
␴k ⳵x j
+ Pk − ␲k + Sk 共6兲
Table 2 Gk term for all models

␳
⳵z
⳵t
+ ␳ūi
⳵z
=
⳵
⳵xi ⳵xi
冋冉 冊 册␮+
␮t ⳵z
␴z ⳵x j
+ Pz − ␲z + Sz 共7兲
Gk

␮ t = C ␮ f ␮␳ k l
1/2
共8兲
IL model
共
−g j␤␳ −
3 C␮ f ␮
2 ␴T
k / ⑀ ui⬘uk⬘
⳵T
⳵xk 兲
共GGDH兲

where Pk = −␳ui⬘u⬘j 共⳵ui / ⳵x j兲, being the generation of turbulent ki-

netic energy due to the shear stress, which plays an important role
GPC/ WX88/ WX93 model
共
−g j␤␳ −
␮t ⳵T
␴T ⳵x j 兲
共SGDH兲

共 兲关共
␮t ⳵T
关共 兲 兴兲共 兲兴
3
in combination with the generation/destruction of turbulent kinetic PDH+ D model Rt 10
−g j␤␳ 1 − exp 1+
energy due to the buoyancy effects 共Gk兲, as explained in this ␴T ⳵x j 12 R3.25
t
paper. 共SGDH+ Damp. f.兲
As exposed, the turbulent or eddy viscosity is directly propor-

072501-2 / Vol. 130, JULY 2008 Transactions of the ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/28/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 Table 3 Cavity 1 „A = 30, RaH = 2.43Ã 1010…. WX93 convergence-estimators using five meshes: m1, m2, m3, m4, m5 „10Ã 10 CV,
20Ã 20 CV, 40Ã 40 CV, 80Ã 80 CV, 160Ã 160 CV…. Acceptable converged solution „p between 1 and 2, high percentage of Richard-
son nodes….

v* = v̄ / uref T* = T̄ / ⌬Tref k* = k / uref

2

Grid Rn GCI* Rn GCI* Rn GCI*

n3 / n2 / n1 共%兲 p 共%兲 共%兲 p 共%兲 共%兲 p 共%兲

m1 m2 m3 83 1.6 2.4E − 01 88 0.5 3.2E + 00 70 1.1 7.2E − 03

m2 m3 m4 93 1.7 7.2E − 02 90 1.0 5.7E − 01 87 1.0 3.6E − 03
m3 m4 m5 91 1.4 3.6E − 02 95 1.0 3.0E − 01 94 1.2 1.4E − 03

pothesis 共GGDH兲 to account for temperature gradients, which are the grid convergence index 共GCI兲 proposed by Roache 关20,21兴.
orthogonal to the gravity vector. PDH+ D applies the usual SGDH This technique is useful to establish criteria on the sensitivity of
adding a damping function, which improves numerical behavior the simulation to the computational model parameters, which ac-
and model accuracy in transitional natural induced flows. See de- count for the discretization: the mesh spacing and the order of
tails in Refs. 关2,7兴. accuracy of the numerical solutions 共observed order of accuracy
As it can be inferred from general formulation, current k-␻ p兲. The error band, where the independent grid solution is ex-
models do not present any specific buoyancy term for ␻-equation. pected to be contained 共uncertainty due to discretization, GCI兲,
In fact, this is in agreement with different authors, e.g., Markatos and criteria on the credibility of these estimations are also
et al. 关12,13兴 共who stated there was evidence that buoyancy effect achieved by this process. Local estimators of the GCI and p are
should be reflected only in the k-equation兲, Heindel 关14兴 共who calculated at the grid nodes where their evaluation is possible.
detected a negligible effect of constant c3兲, or Peng and Davidson These grid nodes are called Richardson nodes. Global values of
关2兴 共who directly fixed c3 = 0 for PDH+ D model兲. Additionally, GCI and p are obtained by means of volumetric averaging. An
different computations have been carried out by present work and estimation is considered credible when the global observed order
a negligible effect of the c3 constant has been obtained for k-␻ of accuracy p approaches its theoretical value, and when the per-
models. On the contrary, for IL and GPC models, the use of c3 centage of Richardson nodes is high enough.
= tanh兩 u / v兩 has been maintained, following suggestions made by More concretely, present work has generally employed a series
Henkes 关17兴 in order to achieve Rodi’s reasonings 关18兴 共c3 should of five meshes 共m1 , . . . , m5兲. The order of accuracy p, the percent-
be close to unity in the vertical boundary layer and close to 0 in age of Richardson nodes, and the normalized error estimation
the horizontal boundary layer兲 关2兴. have been checked for each group of three consecutive meshes.
This is to say, three consecutive meshes define a monotonic nu-
merical convergence behavior. Consequently, an estimation for the
3 Numerical Aspects final solution can be obtained if the percentage of Richardson
nodes is high enough and the reached order of accuracy is close to
3.1 Numerical Procedure. The set of governing partial dif-
its theoretical value 共which depends, essentially, on the discretiza-
ferential equations 共Eqs. 共1兲–共3兲 and Eqs. 共6兲 and 共7兲兲 have been
tion for the convective and the diffusive terms兲. The finest mesh
solved using finite volume techniques. Two different structured
solution 共for each group of meshes兲 is compared to the respective
grids have been employed: staggered meshes 共where velocities are
estimated independent grid solution. This difference is normalized
calculated on the control volume faces兲 and collocated meshes
by a reference value and presented in Tables 3 and 5–7 as a
共where velocities are calculated on the control volume nodal
percentage.
points using the Rhie–Chow interpolation method兲. Fully implicit
As seen, this technique is able to detect nonmonotonic numeri-
time integration has been applied and a pressure based method of
cal convergence 共when the percentage of Richardson nodes is low,
the Semi-implicit method for pressure-linked equations 关19兴
or when the order of accuracy 共p兲 is far from the expected value兲.
共SIMPLE兲 family is used to solve the velocity-pressure field cou-
pling. Central differences are employed for the evaluation of dif- On the other hand, when both parameters are placed between
fusion terms and convective terms are discretized by means of acceptable values, a credible error band where the independent
UPWIND and/or SMART schemes for all test cases. A multigrid grid solution is expected to be contained 共uncertainty due to dis-
iterative solver is used to solve the algebraic linear system of cretization, GCI兲 is reliable. If the error band is relatively small,
the final solution is acceptable. Otherwise, one can continue the
equations. For k and ⑀ /˜⑀ or ␻ equations, source terms have been
h-refinement until the error band is between acceptable limits. As
linearized in the usual form 关19兴 to adding inconditional positive
shown afterwards, nonmonotonic numerical convergence could
values avoiding numerical instabilities. All simulations have been
reveal numerical problems, which do not allow the model to con-
carried out using a pseudotransient iterative algorithm, applying
verge toward an independent grid solution.
the highest time step, which has ensured convergence, which, gen-
Present work has encountered all these situations. This numeri-
erally, depends on the case and on the model. The coarsest grid of
cal tool has been very effective to corroborate numerical prob-
10⫻ 10 CV volumes has been used, doubling number of CV as
lems, which were already pointed in diverse precedent publica-
grid refinement is applied, and the finest mesh of 160⫻ 160 CV
tions 共i.e., Refs. 关2,8兴兲 and very helpful to ensure the numerical
has been chosen.
accuracy of present final solutions. See Refs. 关20–22兴 for details.
In order to correctly reproduce the rapid gradients, which take
place near walls, the adopted meshes have been designed using a
hyperbolic tangent stretching 共detailed in Ref. 关8兴兲 close to the 4 Test Cases and Results
isothermal and adiabatic boundaries. For all presented cavity
flows, a factor of 3.0 has been applied for the vertical isothermal Results computed with previously exposed turbulence models
walls, and a factor of 2.5 for the horizontal adiabatic ones. are presented in this section. They are discussed through compari-
son with experimental and/or DNS data. Various cavity configu-
3.2 Numerical Solution Verification. Numerical solutions rations have been used in several experimental and numerical
are verified using a postprocessing procedure based on the gener- works available in literature. Three different test cases have been
alized Richardson extrapolation for h-refinement studies and on chosen in order to investigate the behavior of different models

Journal of Heat Transfer JULY 2008, Vol. 130 / 072501-3

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/28/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 when solving the turbulent natural convection. The differentially Table 4 Cavity 1 „A = 30, RaH = 2.43Ã 1010…. The difference „in
heated cavities of aspect ratios 30:1, 5:1, and 4:1, in combination percentage… with respect to the experimental data †15‡ „Nu
with their based-on-height Rayleigh numbers 共RaH = 2.43⫻ 1010, = 149, ␮t* at center= 30.4….
5 ⫻ 1010, and 1 ⫻ 1010, respectively兲 have been selected as they
present different turbulent natural convection phenomena, which Turbulence model Nu ␮t* at center
is interesting in order to test turbulence models’ abilities and limi-
IL model 150.65共1%兲 27.11共−12% 兲
tations. For all cases, the vertical walls are isotherms 共at different
GPC model 131.31共−12% 兲 22.56共−25% 兲
temperatures兲 while the top and bottom walls are adiabatic. As WX88 model 135.87共−9 % 兲 24.33共−20% 兲
shown, major problems take place when sudden relaminarization WX93 model 133.67共−10% 兲 21.83共−28% 兲
occurs in the core of the 4:1 and 5:1 cavities, and when turbulence PDH+ D model 152.75共2.5%兲 24.69共−19% 兲
regime is moderate or low. Hence, cavities of aspect ratios 5:1 and
4:1 present the most difficult aspects for turbulence models as
detailed in this study.

4.1 Natural Turbulent Convection in a 30:1 Heated agreement with the experimental data. This was expected as they
Cavity. This case corresponds to a differentially heated tall cavity are models designed for forced convection cases. The PHD+ D
with an aspect ratio of 30 共H / W兲 and a Rayleigh number based- model provides these improvements by redesigning model con-
stants, adding a cross-diffusion term for ␻-equation and a damp-
on-height cavity H 共RaH = g␤⌬TrefH3Pr/ ␯2兲 of 2.43⫻ 1010. This is
ing function in Gk calculus.
a natural convection cavity flow generated with two isothermal
However, all models have reasonably given accurate results.
walls 共east and west兲, and two adiabatic walls 共north and south兲.
The reason could be probably found in the geometry of the case.
The work of Daffa’alla and Betts 关15兴 presents experimental re-
More concretely, the tending-to-infinity cavity aspect ratio en-
sults, which are used to check accuracy of the models.
sures, as has been previously introduced, the collapse of both
This case has been solved for all models using a collocated
turbulent boundary layers 共from hot and cold walls兲, showing a
discretization 共able to solve complex geometries兲 and UPWIND
complete turbulent core of the cavity. As two-equation turbulent
as a numerical scheme for the convective terms. The influence of
models have been deeply suited to forced convection cases, typi-
the spatial discretization and the numerical scheme for convective
cal transition from laminar to turbulent flow due to the wall-
terms have been tested through comparison with two additional
damping effects 共i.e., the resolution of viscous sublayer and the
simulations: First, IL model over a staggered grid and, second,
normal-to-wall increasing turbulence intensity兲 has been qualita-
SMART for convective terms. Based on a mesh of 160⫻ 160 tively well captured for all low-Reynolds number models, as well
control volumes, the first comparison 共collocated versus stag-
as for the k-␻ WX88 model, which in spite of its high-Reynolds
gered兲 shows that the maximum relative difference 共using
number model condition is able to provide both viscous and iner-
nondimensional velocities, temperatures, and lengths兲 关8兴 are of
tial layer solutions, thanks to the use of ␻ as dissipation quantity.
0.1% for ū, 0.7% for v̄, and 0.4% for T̄, while the relative differ- Next cases, however, will show different geometries and, con-
ences in turbulent quantities are even smaller. The second com- sequently, different flow structures. Two-equation turbulence
parison 共UPWIND versus SMART兲 has shown maximum differ- models will be checked under relaminarization produced by the
ences of 0.4% for ū, 4.0% for v̄, and 4.5% for T̄ and once again geometry, as now shown.
smaller differences for the turbulent quantities, which was ex-
pected. Since the differences are relatively small and the flow is 4.2 Natural Turbulent Convection in a 5:1 Heated Cavity.
buoyancy affected, UPWIND scheme seems to be good enough to This case corresponds to a differentially heated rectangular cavity
study the behavior of turbulence models to predict turbulent natu- with an aspect ratio of 5 共H / W兲 and a Rayleigh number based-on-
ral convection in the present case. height cavity H 共being RaH = g␤⌬TrefH3Pr/ ␯2兲 of 5 ⫻ 1010. The
The influence of grid refinement is tested in the evolution of boundary conditions of the case follows the same pattern than
global parameters as Nu. Moreover, the results of the local esti- presented for 30:1 cavity: Two vertical isothermal walls and adia-
mators p and GCI on the convergence show that the numerical batic top and bottom walls. With reference to accuracy of the
results are credible using 160⫻ 160 CV for all models, spatial models, the experimental data from Cheesewright et al. 关16兴 have
discretizations, and numerical schemes. As an example, the results been used to analyze turbulence models’ predictions.
of the verification process for WX93 model using UPWIND In order to test the influence of the numerical scheme, another
scheme over a collocated grid are shown in Table 3 for v*, T*, and comparison 共UPWIND versus SMART兲 has been done, this time
k *. using WX93 model. Higher differences 共2% for ū, 6.75% for v̄,
As shown in Table 3, the order of accuracy 共p兲 has the expected and 4.6% for T̄兲 than the first cavity problem 共30:1 aspect ratio兲
value when using UPWIND scheme 共approximately 1兲, the GCI have been obtained, probably due to the recirculating effects. In
values are low enough in last mesh, and the percentage of Rich- spite of this, UPWIND has been the numerical scheme used in all
ardson nodes is high enough. All these estimators conclude that simulations. Their results are acceptable to predict asymptotic
160⫻ 160 CV are sufficient to obtain a credible numerical solu- transition point while this scheme allows faster and more stable
tion. The convergence behavior of the other simulations have been numerical convergence because of its dissipative essence 关23兴.
similar. PDH+ D model has passed the verification process only Historically, this case has reported different numerical problems
with four meshes, since it has not been possible to reach conver- when using low-Reynolds number turbulence models. Some re-
gence with the coarsest mesh 共10⫻ 10 CV兲, as this model uses a searchers like Peng and Davidson 关2兴, Pérez-Segarra et al. 关8兴,
cross-diffusion term in ␻-equation, which implies the need of a etc., found that two-equation eddy-viscosity models were inca-
finer initial grid. A comparison of different turbulence models and pable to reproduce this case with realistic behavior, due to the
the experimental results reported by Daffa’Alla and Betts 关15兴 is encounter of numerical problems.
presented in Table 4. More concretely, Peng and Davidson 关2兴 detected and corrected
The heat transfer rate 共mean Nusselt number兲 is underpredicted a grid degradation of buoyancy effects 共Gk兲. In fact, as grid is
for all models, except IL and PDH+ D. IL is the model that pre- refined, SGDH approach is related to an increasing overestimation
sents better behavior, as its results are relatively close to the ex- of Gk, which 共since this case presents a stable stratification in the
perimental ones 共see Fig. 1共b兲兲. Its special treatment of the Gk core of the cavity兲 becomes a sink term 共gravity vector and tem-
term 共GGDH兲 improve predictions. It is interesting to say that the perature gradient vector have the same direction and sign兲. There-
PDH+ D model notoriously enhances accuracy showed by the fore, the increasing Gk sink term provokes the creeping evolution
other k-␻ models 共WX88 and WX93兲, which present poorer to a laminar solution as the grid is refined, which obviously dis-

072501-4 / Vol. 130, JULY 2008 Transactions of the ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/28/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 IL
500 0.1 GPC
WX88
WX93
PDH+D
400 IL 0.05 Exp.[15]
GPC
WX88
WX93
Nu

300

v*
PDH+D 0

200
-0.05

100
-0.1

0.25 0.5 0.75 0 0.25 0.5 0.75 1

(a) Y/H (b) X/W

Fig. 1 Cavity 1 „A = 30, RaH = 2.43Ã 1010…. Nu prediction „a… and v̄ at y = H / 2 prediction „b… for all the models.

torts accuracy and disables computation credibility. This effect predict such natural convection cavity.
can be checked in Fig. 2, where the WX93 model with SGDH and Unfortunately, there are some other problems related to this
the PDH+ D model 共which uses SGDH with a damping function case. Another “numerical” laminarization 共turbulence intensity
f g兲 are compared. As it is visible, normalized buoyancy effect suppression due to numerical effects兲 takes place when using dif-
共i.e., Nkb = Gk / C␮k␻兲 obtained using SGDH becomes degraded ferent low-Reynolds number k-⑀ models. Despite the fact that
when refining. However, when using f g damping function other researchers have reported interesting results for such models
asymptotic buoyancy effects are obtained, and a final turbulent when solving 5:1 cavity 共see Refs. 关7,14,25,26兴兲, this present work
prediction is achieved. It is interesting to note that the WX88 has encountered critical numerical problems trying to solve the
model does not present this grid degradation, and its verification mentioned cavity with Ince–Launder model 关7兴 共besides other low
estimators are placed in the desired range. Hence, SGDH grid
Reynolds k-⑀ models not reported here, as the Launder–Sharma
sensitivity seems to be closely related to damping functions
model 关27兴兲. There are some differences between this numerical
共which are applied on WX93 and GPC models兲, and do not appear
in the WX88 model. This is in agreement with the work of Peng problem and the previous SGDH grid-sensitivity drawback. First,
and Davidson 关2兴, who encountered the same problem using an- this laminarization process presents a complete turbulence reduc-
other low-Reynolds number model 共PDH兲. However, and despite tion 共␮t = 0兲. Second, it affects all the cavity 共␮t = 0 ∀ x , y兲. Third it
the relative advantage of the WX88 model, its design for high- abruptly appears while convergence process is running. After dif-
Reynolds number phenomena disables it to correctly predict the ferent attempts, computations have shown that by removing
transition, as has been said, and it is not an acceptable model to damping part of dissipation from ˜⑀-equation 共f 2 function兲, conver-

4 4

2
2

0
0

-2
-2
Nks, Nkb
Nks, Nkb

-4 Nks m5
-4
Nks m4
-6 Nks m5 Nks m3
Nks m4 -6 Nks m2
Nks m3 Nkb m5
-8 Nkb m5 Nkb m4
Nkb m4 -8 Nkb m3
-10 Nkb m3 Nkb m2
-10
-12
-12
0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2
(a) x/W (b) x/W

Fig. 2 Cavity 2 „A = 5, RaH = 5 Ã 1010…. Comparison of the normalized buoyancy and shear production effects for WX93 and PDH
+ D model. While for WX93 model, Nkb increases as grid is refined „giving an asymptotic laminar cavity… „a…: PHD+ D model is able
to give an asymptotic prediction, thanks to the stable Nkb prediction for m5 mesh „b….

Journal of Heat Transfer JULY 2008, Vol. 130 / 072501-5

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/28/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 4 4

2 2

0 0

-2 -2

Nks, Nkb
Nks, Nkb

Nks m5
-4 Nks m3 -4 Nks m4
Nks m2 Nks m3
Nks m1 Nkb m5
-6 -6 Nkb m4
Nkb m3
Nkb m2 Nkb m3
-8 Nkb m1 -8

-10 -10

-12 -12
0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2
(a) x/W (b) x/W

Fig. 3 „a… Nks, Nkb prediction for m1, m2, and m3 meshes for IL model under Cavity 2 „A = 5, RaH = 5 Ã 1010…. „b… Nks, Nkb prediction
for m5, m4, and m3 meshes for GPC model under Cavity 3 „A = 4, RaH = 1 Ã 1010….

gence has been possible and a turbulent solution has been lence destruction due to the buoyancy and suppresses the turbu-
achieved, obviously presenting a distorted prediction. Therefore, it lence energy diffusion from the boundary layer to the neighboring
seems plausible to think that 共because of the relaminarization, outer region, avoiding delayed prediction of transition onset, as
which takes place in the core of the cavity兲 dissipation of ˜⑀ is aforementioned 关2兴. Additionally, this damping function corrects
affected by low Ret numbers and c˜⑀2 is changed from 1.92 to the asymptotic behavior for Gk near the vertical wall 关2兴, since
1.344 共70% of high-Reynolds number version兲 in this region. Gk ⬀ x2 as x → 0 关28兴.
Since production of ˜⑀ is not dampened 共P˜⑀ = c⑀1˜⑀ / kPk兲, an increas- In Fig. 4, differences between WX93 and PDH+ D model are
ing sink effect on k-equation is produced through an increasing ˜⑀ shown. As to this point, obtained convergence estimators have to
prediction. As the grid is refined, this effect is intensified and, be commented: As to the results presented for 30:1 cavity 共where
finally, it returns to a sudden and unstable convergence point, all the models returned expected and verified convergence estima-
which “switches off” all turbulent variables. tors兲, for 5:1 cavity, only the PDH+ D model shows acceptable
With reference to the present work, three initial meshes have results. This is shown in Tables 5 and 6. Thus, this fact shows that
been converged for the IL model 共i.e., 10⫻ 10, 20⫻ 20, and 40 the nonasymptotic turbulent prediction is ratified by the conver-
⫻ 40 CVs兲, while for the fourth mesh 共80⫻ 80 CV兲, it has not gence estimators for WX93 and GPC models. The estimators for
been possible to obtain convergence with nontrivial turbulent val- PDH+ D model have the expected values, which is in agreement
ues, which was also pointed out in Ref. 关8兴. As to the grid- with its stable transition point prediction. The PDH+ D model
sensitivity problem, Fig. 3共a兲 共made with m1, m2, and m3 shows a good experimental agreement in this case. It is shown in
meshes兲 shows how the grid-sensitivity problem is not based on Fig. 4共b兲, where the experimental Nu from Cheesewright et al.
Gk prediction 共there is no degradation of Nkb with grid refine- 关16兴 is plotted.
ment兲, as the IL model is using the GGDH, which is not related to
the SGDH grid-sensitivity problem, as has been explained in Refs. 4.3 Natural Turbulent Convection in a 4:1 Heated Cavity.
关2,7兴. Present IL Gk evolution ratifies this aspect. Additionally, Finally, an analysis of a 4:1 cavity at moderate Rayleigh number
detailed effect of f 2 damping function is in agreement with sug- based-on-height cavity H 共RaH = g␤⌬TrefH3Pr/ ␯2 = 1 · 1010兲 has
gestions made by Davidson in Ref. 关3兴, where a correction on f 2 been carried out. This case could be so-called “moderate” since
damping function is proposed in order to modify the decay of first observed chaotic behavior is visible at RaH = 2.3⫻ 108 关29兴.
dissipation in transitional regions. The DNS results from Trias et al. detailed in Ref. 关24兴 have been
Hence, taking into account the grid degradation of the SGDH used to check abilities and accuracy of present models, as well as
for buoyancy effects and the sudden laminarization of some k-⑀ to understand limitations of turbulent natural convection modeling
turbulence models due to grid refinement, this work has not been by means of RANS eddy-viscosity techniques.
able to present converged solutions for the IL model. For other On the one hand, with reference to relaminarization and geom-
models 共GPC, WX88, WX93, PDH+ D兲, converged solutions etry aspect, this case presents, a priori similar phenomena to the
have been obtained. As exposed, they are in agreement with nu- 5:1 cavity, which could be interesting in order to ratify and/or
merical experiments carried out by other authors 关2兴. In fact, GPC refine comments on previously observed grid-sensitivity behavior.
and WX93 models return turbulent solutions but present nonas- On the other hand, the moderate Rayleigh number applied to this
ymptotic behavior due to the SGDH grid-sensitivity problem. The geometry will bring clarity to the abilities of the models and their
WX88 model provides an asymptotic and completely turbulent relationship with turbulence intensity 共and, consequently, will
cavity and, therefore, is not capable to predict transition. On the give information on the range of RaH where these kinds of models
contrary, the PDH+ D model is able to solve this problem due to can become accurate tools兲.
the damping function incorporated in Gk source term. This is the In fact, with reference to numerical problems 共grid-sensitivity
only model that gives required asymptotic numerical solutions. and numerical laminarization兲, this case has presented similar dif-
The model offers a reasonably accurate solution and acceptable ficulties than previously treated 5:1 cavity. The IL model has re-
numerical behavior. ported the same laminarization problem trying to solve even using
In fact, mentioned damping function 共f g兲 constrains the turbu- a coarser mesh 共40⫻ 40 CV兲, presumably due to the lower turbu-

072501-6 / Vol. 130, JULY 2008 Transactions of the ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/28/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 800 800

700 700

600 600 m5
m4
m3
500 500 m2
m1
m5 Exp. hot wall [16]
Nu

Nu
400 m4 400 Exp. cold wall [16]
m3
m2
300 m1 300

200 200

100 100

0.25 0.5 0.75 0.25 0.5 0.75

(a) Y/H (b) Y/H

Fig. 4 Cavity 2 „A = 5, RaH = 5 Ã 1010…. Comparison of the Nu evolution for WX93 and PDH+ D model. On the one hand, it is shown
that the WX93 model is not capable of giving asymptotic behavior in the prediction of transition point when refining the mesh „a….
On the other, the PDH+ D model is able to give an asymptotic transition point prediction „b….

lence intensity. Again, a turbulent solution has been achieved by stable at lower turbulence intensity. It is interesting to note that
blocking damping effects on dissipation of ˜⑀. This fact reinforces this effect has not been observed using the WX93 model, which is
suggestions made for 5:1 cavity with reference to the mathemati- paradigmatic of the very sensitive relationship between SGDH
cal problem of the damping on ˜⑀-equation in IL model. Further- and damping functions of the model.
more, SGDH grid-sensitivity drawback has been encountered The PDH+ D model has presented the most interesting numeri-
when using WX93 model, and a progressive degradation of the cal behavior, returning a turbulent solution and a laminar-to-
turbulent prediction has been observed. Regarding WX88 model, turbulent transition, which is in agreement with related verifica-
the grid-sensitivity problems have been avoided, but accuracy and tion estimators. These estimators have revealed acceptable orders
prediction of transition point from laminar to turbulent flow has of numerical accuracy 共p兲 and high percentages of Richardson
not been achieved due to the model limitations. nodes for all variables, as shown in Table 7. In fact, this is in
Unlike 5:1 cavity, the GPC model has returned a turbulent pre- agreement with previous conclusions observed in 5:1 cavity.
diction showing acceptable verification estimators. In this case, Hence, only the GPC and PDH+ D models have revealed ac-
the GPC model has avoided a degradation of the Gk term when ceptable numerical solutions under used grids and applied
grid refinement has been applied, as shown in Fig. 3共b兲. It seems pseudotransient numerical algorithm. The other models 共IL,
that the coupling between damping functions and SGDH depends WX88, and WX93兲 are related to different numerical and/or de-
on the turbulence intensity and on RaH number, and it becomes sign problems, which do not allow them to solve this case.

Table 5 Cavity 2 „A = 5, RaH = 5 Ã 1010…. PDH+ D convergence estimators using five meshes: m1, m2, m3, m4, m5 „20Ã 20 CV, 40
à 40 CV, 80à 80 CV, 160à 160 CV, 320à 320 CV…. Acceptable converged solution „p between 1 and 2, high percentage of Richard-
son nodes….

v* = v̄ / uref T* = T̄ / ⌬Tref k* = k / uref

2

Grid Rn GCI* Rn GCI* Rn GCI*

n3 / n2 / n1 共%兲 p 共%兲 共%兲 p 共%兲 共%兲 p 共%兲

m1 m2 m3 81 0.3 2.2E + 00 94 −0.1 3.3E + 01 88 1.4 5.7E − 03

m2 m3 m4 82 0.8 9.4E − 01 95 0.2 1.5E + 01 95 1.2 7.8E − 03
m3 m4 m5 86 1.4 1.8E − 01 98 1.3 8.1E − 01 92 2.1 1.3E − 03

Table 6 Cavity 2 „A = 5, RaH = 5 Ã 1010…. WX93 convergence estimators using five meshes: m1, m2, m3, m4, m5 „10Ã 10 CV, 20
à 20 CV, 40à 40 CV, 80à 80 CV, 160à 160 CV…. The converged solution has not passed the verification process „p is not between
1 and 2, small percentage of the Richardson nodes….

v* = v̄ / uref T* = T̄ / ⌬Tref k* = k / uref

2

Grid Rn GCI* Rn GCI* Rn GCI*

n3 / n2 / n1 共%兲 p 共%兲 共%兲 p 共%兲 共%兲 p 共%兲

m1 m2 m3 66 0.3 7.8E + 00 86 0.0 1.2E + 03 83 1.2 9.5E − 03

m2 m3 m4 72 1.6 2.3E − 01 95 0.6 4.6E + 00 90 2.4 1.4E − 03
m3 m4 m5 65 0.1 2.0E + 01 89 −1.4 1.6E + 01 64 1.4 8.8E − 03

Journal of Heat Transfer JULY 2008, Vol. 130 / 072501-7

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/28/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 Table 7 Cavity 3 „A = 4, RaH = 1 Ã 1010…. PDH-D convergence estimators using five meshes: m1, m2, m3, m4, m5 „10Ã 10 CV, 20
à 20 CV, 40à 40 CV, 80à 80 CV, 160à 160 CV…. Acceptable converged solution „p between 1 and 2, high percentage of Richardson
nodes….

v* = v̄ / uref T* = T̄ / ⌬Tref k* = k / uref

2

Grid Rn GCI* Rn GCI* Rn GCI*

n3 / n2 / n1 共%兲 p 共%兲 共%兲 p 共%兲 共%兲 p 共%兲

m1 m2 m3 79 0.4 1.1E + 00 96 0.3 1.0E + 01 89 1.2 5.1E − 03

m2 m3 m4 82 0.0 2.0E + 02 94 −0.3 1.6E + 01 94 0.5 2.2E − 02
m3 m4 m5 86 1.7 7.2E − 02 98 1.6 5.1E − 01 93 2.4 6.8E − 04

With reference to accuracy, several grades can be commented intensity. With reference to y = 15H / 16 plane, turbulence intensity
thanks to the detailed information provided by the DNS data: first, is larger and differences are higher. Performance is admissibly
averaged values 共Nu兲; second, velocity and temperature fields; accurate. On the other hand, PDH+ D model offers worse predic-
third, turbulent statistics as k and ⑀ or C␮␻k. As was expected, the tions due to an overprediction of turbulence intensity. The results
deeper the grade, the higher the differences. As far as the Nusselt for y = H / 2 were especially distorted.
number is concerned, Fig. 5共b兲 shows the results for GPC and Finally, and as to turbulence statistics, two aspects have to be
PDH+ D models in comparison to DNS data. While the PDH
addressed. On the one hand, turbulent kinetic energy 共k兲 predic-
+ D model was well suited for previous cavities 共5:1 cavity at
RaH = 5 ⫻ 1010 and 30:1 cavity at RaH = 2.43⫻ 1010兲, this 4:1 cav- tions can be observed in Fig. 8. As shown, both GPC and PDH
ity at RaH = 1 ⫻ 1010 seems to be too “transitional” 共or too low + D models quantitatively present results comparable to the DNS
turbulent兲 to be well predicted by this model. In fact, it is easily data. As to their quality, GPC model offers a transition point quite
observed how transition from laminar to turbulent flow is moved close to DNS predictions, while PDH+ D predicts a very early
to an extremely lower height. On the other hand, GPC model transition point and presents a strongly larger turbulent boundary
共which reported acceptable accurate predictions for 30:1 cavity, layer. As can be observed in Fig. 8共a兲, it implies a strong overpre-
and which presented numerical problems due to SGDH grid- diction in y = H / 2 plane. This explains the distorted predictions for
sensitivity problems when solving 5:1 cavity兲 has shown a reason- Nu, and specifically, for u*, v*, and T* in y = H / 2 plane. However,
ably accurate behavior solving this 4:1 cavity, obtaining a Nu both models relatively present acceptable turbulent boundary-
prediction relatively close to the DNS data. Transition point is layer normal-to-wall height predictions. On the other hand, with
slightly underpredicted. It is interesting to note that the GPC pre- reference to ⑀ predictions, both models provide inaccurate predic-
diction is very close to the Nu obtained from a previous DNS 2D
tions 共as can be seen in Fig. 9兲, which was expected due to the
data carried out by Trias et al. 关24兴, which could reveal three-
dimensional effects, which are impossible to be observed by limitations of ⑀ / ␻-equation. GPC model offers a qualitatively
RANS-EVM models. good behavior for both y = H / 2 and y = 15H / 16 planes, while
Related to the second grade of comparison, Figs. 5–7 present PDH+ D model suffers a strong overprediction for both planes
the GPC and PDH+ D performances in comparison to DNS data once again. See Fig. 9.
关24兴 for temperature and velocities at y = H / 2 and y = 15H / 16 Hence, an accurate resolution of k has an immediate effect on
planes. In agreement with which was observed for Nu prediction, both temperature and velocity fields. This allows to obtain accept-
GPC model offers interesting predictions 共very close to DNS able transition point prediction, which is of great interest for en-
data兲, especially for y = H / 2 where there is a very small turbulence gineering purposes.

1 450

0.9
400
0.8
350
0.7
300
0.6 DNS [24]
250 PDH+D
Nu

0.5
T*

GPC
200
0.4

T* DNS [24] 150

0.3
T* PDH+D
T* GPC 100
0.2

0.1 50

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x/W y/H

(a) (b)

Fig. 5 Cavity 3 „A = 4, RaH = 1 Ã 1010…. „a… Dimensionless temperature profiles: y = H / 2. „b…Nu number prediction for GPC and
PDH+ D models. Sudden rising of Nu number indicates transition from laminar to turbulent flow.

072501-8 / Vol. 130, JULY 2008 Transactions of the ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/28/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 0.2
0.06

0.04

0.1
0.02

*
0
v*

v
-0.02

-0.04 v* DNS [24]

v* DNS [24] v* PDH+D
-0.1 v* PDH+D -0.06 v* GPC
v* GPC
-0.08

-0.2
-0.1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1
x/W x/W

(a) (b)

Fig. 6 Cavity 3 „A = 4, RaH = 1 Ã 1010…. Dimensionless v profiles: „a… y = H / 2; „b… y = 15H / 16.

5 Concluding Remarks shown an acceptable behavior when solving 30:1 cavity. All buoy-
Abilities and drawbacks related to two-equation eddy-viscosity ancy models 共SGDH, GGDH, and GGDH+ f g兲 have presented
models have been analyzed under turbulent natural convection correct numerical behaviors. This fact seems to be in agreement
cases. A detailed formulation based on k-⑀ and k-␻ models has with the original design of two-equation turbulence models, which
been presented and different rectangular cavities have been ad- are mainly devised to predict laminarization produced due to the
dressed with exposed models, offering several numerical and presence of walls under forced convection. On the contrary, solv-
physical aspects. Treated cavities expose different phenomena due ing 5:1 and 4:1 cavities, diverse numerical problems have ap-
to their diverse aspect ratio and turbulence intensity. In order to peared. Currently, these cases present an additional laminarization
quantify model accuracy, experimental and DNS data extracted phenomenon, which is generated by cavity aspect ratio 共in com-
from literature have been used. bination with moderate turbulent intensity兲, instead of wall-
As a result from computations, this work has basically revealed damping effect 共which is the only laminarization process in 30:1
two different subjects. First, numerical problems related to predic- cavity兲. This phenomena lead to two different numerical prob-
tion of turbulent natural convection. Second, skills and abilities of lems: On the one hand, the IL model suffers a sudden laminariza-
such models to accurately predict heat transfer and fluid flow for tion presumably due to a damping effect on dissipation of
the studied cases. ˜⑀-equation, which collapses all turbulent variables while refining
Related to the first item, low-Reynolds number two-equation the mesh. On the other hand, SGDH presents a numerical degra-
models seem to work reasonably well when solving turbulent dation with grid refinement in combination with damping func-
boundary layers under increasing turbulence intensity from wall to tions, which was also highlighted by Peng et al. 关2兴. This fact does
cavity core. In fact, even high-Reynolds number WX88 model has not allow to obtain an asymptotic prediction. Both effects are

2.0E-03 1.0E-02
9.0E-03
1.5E-03 8.0E-03
7.0E-03
1.0E-03 6.0E-03
5.0E-03
5.0E-04
*

4.0E-03
u

3.0E-03
0.0E+00 2.0E-03
1.0E-03
-5.0E-04 u* DNS [24] u* DNS [24]
0.0E+00
u* PDH+D u* PDH+D
u* GPC -1.0E-03 u* GPC
-1.0E-03 -2.0E-03
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x/W x/W

(a) (b)

Fig. 7 Cavity 3 „A = 4, RaH = 1 Ã 1010…. Dimensionless u profiles: „a… y = H / 2; „b… y = 15H / 16.

Journal of Heat Transfer JULY 2008, Vol. 130 / 072501-9

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/28/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 10-3

6.0E-04

-4
10
4.0E-04
k*

k*
-5
10
k* DNS [24] 2.0E-04 k* DNS [24]
k* GPC
k* GPC
k* PDH+D
k* PDH+D

10-6 0.0E+00
0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1
x/W x/W

(a) (b)

Fig. 8 Cavity 3 „A = 4, RaH = 1 Ã 1010…. Dimensionless k profiles: „a… y = H / 2; „b… y = 15H / 16.

related to the models’ inabilities to capture the laminarization pro- numerical limitations of the models. However, they are generally
duced in the core of the cavity. Hence, only the PDH+ D and GPC acceptable tools for different Rayleigh numbers and aspect ratio
models have returned verified numerical computations for 4:1 cavities, taking into account that, depending on the case, specific
cavity while solely PDH+ D model has achieved an acceptable models should be used.
numerical behavior for 5:1 cavity.
Related to the second subject, accuracy is clearly related to the Acknowledgment
flow structure of the case. For 30:1 cavity, all models have shown This work has been financially supported by the Ministerio de
an interesting behavior to predict Nu, and also vertical velocity Eduación y Ciencia, Secretaría de Estado de Universidades e In-
v*, despite an important underprediction of turbulent viscosity in vestigación, Spain 共ref. ENE2006-14247/ALT and ENE2007-
the core of the cavity. On the other hand, regarding the 5:1 cavity 67185兲.
Nu number and transition point, they are well captured by the
PDH+ D model, as well as k distribution along the cavity. Finally,
with reference to the 4:1 cavity, the GPC model has reported very Nomenclature
interesting results, relatively close to the DNS data regarding di- A ⫽ aspect ratio 共A = H / W兲
verse fields. While on the contrary, the PDH+ D model has not C ␮ , c 1 , c 2 , c 3 , c r␻ ⫽ turbulent model constants
been able to reproduce this lower turbulent case and has strongly cp ⫽ specific heat at constant pressure
over predicted turbulent variables. CV ⫽ control volume
Hence, turbulent natural convection prediction by means of f1 , f2 , fk , f␮ ⫽ damping functions of low-Reynolds number
two-equation turbulence models is associated with physical and models

180

102 160

140
1
10
120
0
10 100
ε∗

ε∗

80
10-1
ε∗ DNS [24] ε * DNS [24]
ε∗ GPC 60 ε * GPC
10
-2 ε∗ PDH+D ε * PDH+D
40

10-3 20

0
0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1
x/W x/W

(a) (b)

Fig. 9 Cavity 3 „A = 4, RaH = 1 Ã 1010…. Dimensionless ⑀ profiles: „a… y = H / 2; „b… y = 15H / 16.

072501-10 / Vol. 130, JULY 2008 Transactions of the ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/28/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 H ⫽ cavity height ⫽ dimensionless specific dissipation rate 共␻*
k ⫽ turbulent kinetic energy = ␻H2 / ␯兲
k* ⫽ dimensionless turbulent kinetic energy
共k / uref
2
兲
lref ⫽ reference length 共lref = H兲 References
关1兴 Pope, S., 2000, Turbulent Flows, Cambridge University Press, Cambridge.
l ⫽ turbulent length scale 共l = k3/2 / ⑀ for k-⑀ 关2兴 Peng, S., and Davidson, L., 1999, “Computation of Turbulent Buoyant Flows
models兲 and 共l = k1/2 / ␻ for k-␻ models兲 in Enclosures With Low-Reynolds Number k-␻ Models,” Int. J. Heat Fluid
Nks ⫽ normalized shear stress production of k Flow, 20, pp. 172–184.
共Nks = Pk / ␳C␮k␻ or Pk / ␳⑀兲 关3兴 Davidson, L., 1990, “Calculation of the Turbulent Buoyancy-Driven Flow in a
Rectangular Cavity Using an Efficient Solver and Two Different Low Rey-
Nkb ⫽ normalized buoyancy production of k 共Nkb nolds Number k-⑀ Models,” Numer. Heat Transfer, Part A, 18, pp. 129–147.
= Gk / ␳C␮k␻ or Gk / ␳⑀兲 关4兴 Hanjalić, K., and Vasić, S., 1993, “Computation of Turbulent Natural-
Nu ⫽ Nusselt number Convection in Rectangular Enclosures With an Algebraic Flux Model,” Int. J.
Pk ⫽ production of k due to shear stresses Heat Mass Transfer, 36, pp. 3603–3624.
关5兴 Hanjalić, K., Kenjereš, S., and Durst, F., 1995, “Natural Convection in Parti-
Gk ⫽ production of k due to buoyancy effects tioned Two-Dimensional Enclosures at Higher Rayleigh Numbers,” Int. J. Heat
Pr ⫽ Prandtl number 共Pr= ␮c p / ␭兲 Mass Transfer, 39共7兲, pp. 1407–1427.
关6兴 Kenjereš, S., Gunarjo, S. B., and Hanjalić, K., 2005, “Contribution to Elliptic
p̄ ⫽ the mean pressure Relaxation Modelling of Turbulent Natural and Mixed Convection,” Int. J.
RaH ⫽ Rayleigh number based-on-height length Heat Fluid Flow, 26共7兲, pp. 569–586.
共RaH = g␤⌬TrefH3Pr/ ␯2兲 关7兴 Ince, N., and Launder, B., 1989, “Computation of Buoyancy-Driven Turbulent
Ret ⫽ turbulent Reynolds number 共Ret = k / ␯␻ or Flows in Rectangular Enclosures,” Int. J. Heat Fluid Flow, 10共1兲, pp. 110–
117.
k2 / ␯⑀兲 关8兴 Pérez-Segarra, C. D., Oliva, A., Costa, M., and Escanes, F., 1995, “Numerical
Sij ⫽ mean rate of strain tensor Experiments in Turbulent Natural and Mixed Convection in Internal Flows,”
Int. J. Numer. Methods Heat Fluid Flow, 5共1兲, pp. 13–33.
T̄ ⫽ mean temperature 关9兴 Wilcox, D. C., 1988, “Reassessment of the Scale-Determining Equation for
T⬘ ⫽ turbulent fluctuating temperature Advanced Turbulence Models,” AIAA J., 26, pp. 1299–1310.
Tref ⫽ reference temperature 共Tref = Tcold兲 关10兴 Wilcox, D. C., 1994, “Simulation of Transition With a Two-Equation Turbu-
lence Model,” AIAA J., 32, pp. 247–255.
T* ⫽ dimensionless temperature 关11兴 Goldberg, U., Peroomian, O., and Chakravarthy, S., 1998, “A Wall-Distance-
共T* = T̄ − Tref / ⌬Tref兲 Free k-⑀ Model With Enhanced Near-Wall Treatment,” ASME J. Fluids Eng.,
120, pp. 457–462.
⌬Tref ⫽ reference temperature 共⌬Tref = Thot − Tcold兲 关12兴 Markatos, N., Malin, M., and Cox, G., 1982, “Mathematic Modelling of
t ⫽ time Buoyancy-Induced Smoke Flow in Enclosures,” Int. J. Heat Mass Transfer,
tref ⫽ reference time 共tref = lref / uref兲 25, pp. 63–75.
关13兴 Markatos, N., and Pericleous, K., 1984, “Laminar and Turbulent Natural Con-
ui ⫽ mean velocity in the i-direction vection in an Enclosed Cavity,” Int. J. Heat Mass Transfer, 27, pp. 755–772.
ui⬘ ⫽ turbulent fluctuating velocity in the 关14兴 Heindel, T., Ramadhyani, S., and Incropera, F., 1994, “Assessment of Turbu-
lence Models for Natural Convection in an Enclosure,” Numer. Heat Transfer,
xi-direction Part B, 26, pp. 147–172.
ui⬘T⬘ ⫽ turbulent heat flux 关15兴 Dafa’alla, A., and Betts, P., 1996, “Experimental Study of Turbulent Natural
uref ⫽ reference velocity for 30:1 and 5:1 cavities Convection in a Tall Air Cavity,” Exp. Heat Transfer, 9, pp. 165–194.
共uref = 冑g␤⌬TrefLref兲
关16兴 Cheeswright, R., King, K., and Ziai, S., 1986, “Experimental Data for the
Validation of Computer Codes for the Prediction of Two-Dimensional Buoyant
uref ⫽ reference velocity for 4:1 cavity Cavity Flows,” Proceedings of Significant Questions in Buoyancy Affected
共uref = ␣冑RaH / H兲 Enclosure or Cavity Flows, pp. 75–81.
关17兴 Henkes, R., van Der Vlugt, F., and Hoogendoorn, C., 1991, “Natural Convec-
u ⫽ dimensionless velocity 共ū / uref兲
* tion Flow in a Square Cavity Calculated With Low Reynolds-Number Turbu-
v* ⫽ dimensionless velocity 共v̄ / uref兲 lence Models,” Int. J. Heat Mass Transfer, 34, pp. 377–388.
关18兴 Rodi, W., 1984, “Turbulence Models and Their Application in Hydraulics: A
u␶ ⫽ friction velocity 共u␶ = 冑␶w / ␳兲 State-of-the-Art Review,” University of Karlsrhue, Germany.
W ⫽ cavity width 关19兴 Patankar, S. V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere,
Washington, DC.
xi ⫽ Cartesian coordinate in the i-direction 关20兴 Cadafalch, J., Pérez-Segarra, C. D., Cónsul, R., and Oliva, A., 2002, “Verifi-
y + ⫽ dimensionless distance to the nearest wall cation of Finite Volume Computations on Steady State Fluid Flow and Heat
共y + = yu␶ / ␯兲 Transfer,” ASME J. Fluids Eng., 124, pp. 11–21.
关21兴 Roache, P., 1997, “Quantification of Uncertainty in Computational Fluid Dy-
Greek Symbols namics,” Annu. Rev. Fluid Mech., 29, pp. 123–160.
关22兴 Roache, P., 1998, Verification and Validation in Computational Science and
␣ ⫽ thermal diffusivity 共␣ = ␯ / Pr兲 Engineering, Hermosa, Albuquerque, NM.
␦ij ⫽ Kronecker delta 关23兴 Roache, P., 1972, “On Artificial Viscosity,” J. Comput. Phys., 10, pp. 169–
⑀ ⫽ dissipation rate of k 184.
关24兴 Trias, F. X., Soria, M., Oliva, A., and Pérez-Segarra, C. D., 2007, “Direct
⑀* ⫽ dimensionless dissipation rate Numerical Simulation of Two and Three-Dimensional Turbulent Natural Con-
共⑀* = ⑀ / 共␯ / tref
2
兲兲 vection Flows in a Differentially Heated Cavity of Aspect Ratio 4,” J. Fluid
␭ ⫽ thermal conductivity Mech., 586, pp. 259–293.
关25兴 Murakami, S., Kato, S., Chiakamoto, T., and Laurence, D., 1996, “New Low-
␮ ⫽ dynamic viscosity Reynolds-Number k-⑀ Model Including Damping Effect due to Buoyancy in a
␮t ⫽ eddy or turbulent viscosity Stratified Flow Field,” Int. J. Heat Mass Transfer, 39共16兲, pp. 3483–3496.
␮t* ⫽ dimensionless eddy or turbulent viscosity 关26兴 Liu, F., and Wen, J., 1999, “Development and Validation of an Advanced
Turbulence Model for Buoyancy Driven Flows in Enclosures,” Int. J. Heat
共␮t* = ␮t / ␮兲 Mass Transfer, 42, pp. 3967–3981.
␯ ⫽ kinematic viscosity 共␯ = ␮ / ␳兲 关27兴 Launder, B., and Sharma, B., 1974, “Application of the Energy Dissipation
Model of Turbulence to the Calculation of Flow Near a Spinning Disc,” Lett.
␳ ⫽ density Heat Mass Transfer, 1共2兲, pp. 131–138.
␴T ⫽ turbulent Prandtl number 关28兴 So, R., and Sommer, T., 1994, “A Near-wall Eddy Conductivity Model for
␴k , ␴z ⫽ model constants Fluids With Different Prandtl Numbers,” ASME J. Fluids Eng., 116, pp. 844–
␶ ⫽ turbulent time scale 854.
关29兴 Xin, S., and Quere, P., 1995, “Direct Numerical Simulations of Two-
␻ ⫽ specific dissipation rate of k Dimensional Chaotic Natural Convection in a Differentially Heated Cavity of
␻* Aspect Ratio 4,” J. Fluid Mech., 304, pp. 87–118.

Journal of Heat Transfer JULY 2008, Vol. 130 / 072501-11

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/28/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use