Excerpt from the Proceedings of the COMSOL Conference 2009 Boston
Benchmark Comparison of Natural Convection in a Tall Cavity
Heather E Dillon∗,1 , Ashley Emery1 , and Ann Mescher1
University of Washington
Corresponding author: Seattle, WA 98105, email@example.com
Abstract: A comparison of the commercial computational software .
code COMSOL is performed with the bench- For natural convection studies in a cavity
mark solutions provided by the literature for several parameters are used to classify the sys-
a tall, diﬀerentially heated rectangular cavity tems. The geometry of the cavity is repre-
for aspect ratios of 8, 15, 20, and 33. At small sented by the height H, the width L and the
Rayleigh numbers the ﬂow is dominated by aspect ratio (A).
conduction. As the Rayleigh number is in-
creased the ﬂow becomes unstable, ﬁrst result- H
ing in multicellular secondary ﬂow patterns, L
and then as the Rayleigh number is further For the computational work described in
increased becoming turbulent. For this work this paper the range of aspect ratios consid-
a range of Rayleigh numbers from 1e3 to 3e5 is ered was A = 8 − 33.
considered for a Prandtl number of 0.71. The Several dimensionless parameters are used
predicted wavenumber, critical Rayleigh num- to characterize the ﬂow. The Rayleigh num-
ber, and characteristic streamfunctions are ber (Ra) is deﬁned below and is often used
compared with experimental work reported in buoyancy ﬂows to characterize the transi-
in the literature. The eﬀect of COMSOL’s tion between conduction dominated ﬂow and
anisotropic and crosswind numerical diﬀusion convection dominated ﬂow. The Prandlt num-
values are also evaluated. ber (P r) is the ratio of the viscous diﬀusion
and thermal diﬀusion. For this work, only the
Keywords: benchmark, comsol, natural con- value of P r = 0.71 (representative of air) is
vection, tall cavity considered, and a range of Ra = 1e3 − 3e5 is
1 Introduction ρ2 gcp β(∆T )L3
Ra = (2)
To benchmark the natural convection mod- kµ
eling capabilities of COMSOL, comparisons cp is the heat capacity of the gas, g is the
were made for free convection in tall rectan- acceleration of gravity, β is the isobaric coef-
gular two dimensional cavities. The geometry ﬁcient of thermal expansion, µ is the dynamic
is a cavity of air with isothermal heated and viscosity, k is the thermal conductivity, and ρ
cold walls. The diﬀerentially heated walls cre- is the density.
ate a buoyancy driven convection inside the
cavity. Understanding the performance of the Pr = (3)
COMSOL Multiphysics numerical diﬀusion for k
this case will allow more complex geometries, The critical Rayleigh number (Rac ) is an
including high aspect ratio annular systems, to important value for comparisons to experi-
be modeled with greater conﬁdence. The focus mental work. It is deﬁned as the Rayleigh
is on correctly simulating the ﬂuid behavior for number at which the ﬂow transitions from one
cavities with high aspect ratios (greater than stability regime to another, which corresponds
15). A summary of the computational domain to the transitions from the conduction domi-
for work of this type is shown in Figure 1. nated regime to convection dominated regime.
The study of ﬂow in a thermally driven In this work the critical Rayleigh number char-
cavity is a classical problem that has been ex- acterizes the formation of cells with oscillatory
plored by many researchers since the 1950s. temperatures. Experimental studies often re-
The speciﬁc case of a square cavity has been port the wavenumber (n) at Rac based on the
used extensively as a numerical benchmark for wavelength (λ) of the cells.
Fig. 1: a) Computational domain for the tall cavity problem. b) Literature summary for the tall
cavity problem. Solid and dashed lines represent stability at the Boussinesq limit (-) and higher
temperature diﬀerences (- -). Dot-dash and dotted lines represent boundary layer regimes for the
Boussinesq limit (-.-) and the higher temperature diﬀerences (...). The scope of the current work is
shown for reference. Adapted from .
ties (A = 1 − 5) including de Vahl Davis .
2π For small aspect ratios the ﬂow characteristics
λ are not multi-cellular and are time invariant.
The heat transfer across the cavity is usu-
ally reported in terms of the Nusselt num-
ber (Nu), which is the ratio of the convective 2.1 Experimental Studies
heat transfer coeﬃcient to the conduction heat
transfer coeﬃcient. Several experimental studies have been per-
formed on natural convection in a tall cavity.
Table 1 summarizes the natural convection ex-
2 Literature periments performed. The experimental liter-
ature is more extensive for small aspect ratios
The literature has shown that the stability of
(A < 4).
the ﬂow in a cavity is governed by the Prandtl
number, the Rayleigh number and the geom- Most of the experimental work consisted
etry of the cavity. At small Rayleigh num- of ﬂow visualization. Most of the authors esti-
bers the ﬂow is dominated by conduction. As mated the observed critical Rayleigh number
the Rayleigh number is increased the ﬂow be- and the wavenumber.
comes unstable, ﬁrst resulting in multicellu-
lar secondary ﬂow patterns, and then as the 2.2 Computational Models
Rayleigh number is further increased the ﬂow
becomes chaotic. A high level summary of the The computational models for this problem
expected ﬂow behavior for Cartesian geometry are summarized in Table 1. All the compu-
is shown in Figure 1. tational studies unless otherwise noted used
The focus in most experimental and com- the Boussinesq approximation and none of the
putational work is on small aspect ratio cavi- models considered radiative heat transfer.
Author Year Ra Pr A Description
Elder  1965 1e3 − 1e5 Silicon 1-60 Streak photos showed
Vest and Arpaci  1969 3.7e5 0.71, 1000 20-33 Streak photos of the ﬂow.
Reported Rac .
Vest and Arpaci  1969 7e3 − 3e5 0.71, 1000 20-33 Galerkin method.
Korpela et al  1973 1e2 − 1e4 0-50 ∞ Reported Rac .
Bergholz  1978 1.6e3 − 3e4 0.73-1000 ∞ Galerkin method.
Lee and Korpela  1983 3e4 − 2e5 0-1000 15-40 Reported Nu and
Chenoweth and Paolucci  1986 1e5 − 1e6 0.71 1-10 Compare ideal gas
Chait and Korpela  1989 5e2 − 1.5e4 0.71, 1000 ∞ Pseudo spectral method.
Le Quere  1990 7e3 − 4e4 0.71 16 Explored return to
Liakopoulos et al.  1990 0.71 10-25 Included ﬂux wall conditions.
Facas  1992 9e3 − 6e4 0.71 15 Considered inclined
cavities (75o ).
Suslov and Paolucci  1995 6e3 − 1e4 0.71 ∞ Non-Boussinesq impact
on stability and considered
Rac with ∆T .
Xin and Le Quere  2002 3e5 − 5e5 0.71 8 Benchmark study reported Rac .
Christon et al.  2002 1e5 0.71 8 Comparison study of
methods, grids, etc.
Hayden et al.  2003 7e3 − 3e5 0.71 33 Commercial code FIDAP.
Present work 2009 1e4 − 6e5 0.71 8-33 Commercial code COMSOL.
Table 1: Summary of experimental and computational natural convection studies in a tall cavity.
The computational literature also demon- depending on the speciﬁcs of the spatial dis-
strates that for low Prandtl number ﬂows the cretization, grid resolution, governing equa-
primary instability is a thermal-shear insta- tions, and solver parameters. Other param-
bility which obtains kinetic energy primarily eters such as the time averaged velocity were
from shear production in the ﬂow. Larger calculated more precisely by the participants.
Prandtl number ﬂows experience a thermal- Most importantly, the critical Rayleigh num-
buoyant instability which obtains energy from ber and the predicted Hopf bifurcation were
the buoyant force. computed accurately by all methods.
One of the most interesting papers for the Some authors ( and ) predict an
rectangular geometry is the work of Christon oscillatory time-dependent instability for tall
et al.  which is a compilation of solutions cavities. Lee and Korpela and Liakopoulos
submitted by many participants with a vari-  predict that a stationary instability pre-
ety of solution methods. The authors found cedes the onset of oscillatory convection for
that no speciﬁc method (ﬁnite element, ﬁnite high aspect ratios. Suslov and Paolucci in-
volume, and ﬁnite diﬀerence) provided supe- dicate this observation is simply a product
rior results. The authors also note that the of temperature invariant air properties .
amplitude of the periodic temperature oscil- They also observe that two modes of instabil-
lations can vary by an order of magnitude ity are possible depending on the magnitude
of the temperature diﬀerence between the hot 3.2 Air Properties
and cold wall.
In a similar work Chenoweth and Paolucci The density of the ﬂuid for natural convection
 observed that as the temperature diﬀer- systems is often represented with the Boussi-
ence in the cavity increases, a lower critical nesq approximation. The Boussinesq approxi-
Rayleigh number is found. They conﬁrmed mation is used for computational problems of
that the nature of the instability changes with this type to simplify the formation of the cou-
increasing temperature diﬀerence. pled Navier-stokes equations. The approxima-
tion assumes that density variations are small
in the ﬂuid except in evaluating the buoy-
3 Solution Method ancy force (gravity multiplied by density). In
general, the Boussinesq approximation is only
3.1 Governing Equations valid when the temperature diﬀerences in the
The buoyancy driven ﬂow is modeled as a cou- system are small (less than 28o C).
pled system with ﬂuid motion (Navier Stokes), For the system considered in this work
convection and conduction heat transfer. The the Boussinesq approximation is valid, how-
governing equations are given in Equations 5- ever some authors did observe diﬀerent results
8. when the air was treated as an ideal gas (
The boundary conditions for the system and ). The experimental results compared
are speciﬁed as adiabatic on the upper and in this work involved larger temperature diﬀer-
lower cavity boundaries. The right (hot) and ences. For this reason the Boussinesq approx-
left (cold) walls are constrained as constant imation was not used and the air was treated
temperatures as shown in Figure 1. as an ideal gas.
The weakly compressible Navier Stokes Other ﬂuid properties (k, cp , µ) were all
equation for this system is given in simpliﬁed temperature dependent and were deﬁned by
form. using the polynomial ﬁt of Reeve .
∂u 3.3 Spatial Discretization
ρ + ρ(u · )u = (5)
∂t For this analysis a grid resolution study was
· (−pI + η( u + uT )) conducted for the A = 15 geometry. In the re-
−(2η/3 ·u)I + F ﬁnement study the density of the mesh was
increased and the magnitude and formation
of the temperature oscillations was recorded.
·u = 0 (6) The mesh was optimized when the formation
of the oscillations occurred with the same mag-
·(c · ψ) = 0 (7) nitude and roughly the same iteration time as
the most reﬁned mesh. The ﬁnal mesh was
In this form u represents the velocity vec- 576 triangular quadratic Lagrangian elements
tor of the ﬂuid, ρ is the ﬂuid density, η is the with 5513 degrees of freedom.
dynamic viscosity of the ﬂuid, and F is the
force applied to the ﬂuid. Where c is the dif-
fusion coeﬃcient and ψ is the streamfunction. 4 Numerical Diﬀusion in COMSOL
For the models discussed in this section the Multiphysics
ﬂuid force is speciﬁed in the y-direction for
buoyancy driven ﬂow as Fy = −ρg. Initially a small study was conducted to ex-
The heat transfer convection is governed plore each of the numerical diﬀusion param-
by Equation 8. eters discussed in this work. In the study
a baseline value for each numerical diﬀusion
parameter was taken to be the COMSOL de-
ρcp + · (−k T ) = −ρcp u · T (8) faults for weakly compressible Navier-Stokes.
Each diﬀusion parameter was then adjusted
For this equation k is the thermal conduc- and a qualitative observation about the eﬀect
tivity, cp is the heat capacity at constant pres- of altering the parameters on the temperature
sure, and T is the temperature. oscillations in the system was made. It was
found that to most achieve results most closely sive. SUPG has the largest eﬀect near
matching the experimental data required re- the boundary layers . Because SUGP
ducing slightly the anisotropic and crosswind and GLS are closely related the SUGP
diﬀusion. option was not examined in detail in the
The numerical diﬀusion parameters in baseline study.
COMSOL Multiphysics are outlined in the
COMSOL documentation , however a short • Galerkin least-squares (GLS).
overview of each type is included here along Galerkin least-squares (GLS) is a more
with the results of modifying the values. advanced version of SUPG and is of the
same order of accuracy. In COMSOL,
• Isotropic Diﬀusion. In COMSOL, GLS is the default streamline diﬀusion
isotropic diﬀusion is equivalent to adding mechanism and is part of the baseline
a term to the physical diﬀusion coeﬃ- solution. In the baseline study GLS
cient. The new term is a tuning param- was found to be a more robust solution
eter (δid ) with recommend values less method. In some cases the SUPG failed
than 0.5. The eﬀect of isotropic diﬀusion to converge at all without crosswind dif-
is to dampen oscillations and impede fusion enabled. For this reason GLS
propagation of oscillations . When was the streamline diﬀusion method of
isotropic diﬀusion was introduced in this choice for this work.
study it dampened the oscillations or
• Crosswind Diﬀusion. Crosswind dif-
halted them altogether as expected.
fusion introduces artiﬁcial diﬀusion in
• Anisotropic Diﬀusion. Anisotropic the orthogonal direction to the stream-
diﬀusion is similar to isotropic diﬀusion lines. Crosswind diﬀusion is most use-
method but occurs only in the direc- ful when undershoots and overshoots can
tion of the streamline. Like isotropic occur in the numerical solutions. Cross-
diﬀusion, anisotropic diﬀusion modiﬁes wind diﬀusion methods do not alter the
the equation with a new term. The equation but are non-linear . The ef-
COMSOL documentation suggests val- fect of adjusting the crosswind diﬀusion
ues for anisotropic diﬀusion (δsd ) should coeﬃcient (ck ) in COMSOL for the base-
be less than 0.5 . In the paramet- line study was signiﬁcant. Increasing
ric study conducted for this work the ck to 0.2 reduces the oscillations signiﬁ-
default values of δsd = 0.25 resulted cantly. Decreasing ck from the baseline
in values close to the baseline solution. smooths the oscillations but does not
Increasing the anisotropic diﬀusion re- change the amplitude. A small increase
duced the amplitude of the oscillations in the period of the oscillations occurs.
and had a slight increase of period. Removing crosswind diﬀusion caused the
Decreasing the anisotropic diﬀusion in- solution to require a much longer time to
creased the amplitude of the oscillations stabilize and also reduced the period of
and slightly reduced the period. oscillations.
• Streamline upwind Petrov-
Galerkin (SUGP). Streamline diﬀu-
sion introduces artiﬁcial diﬀusion in the
5.1 Numerical Diﬀusion
streamline direction. The Streamline
upwind Petrov-Galerkin (SUGP) does Based on the results of the study of the ef-
not perturb the original transport equa- fects of modifying the diﬀusion coeﬃcients,
tion. It is closely related to upwinding the subsequent work focused on adjustments
schemes in ﬁnite diﬀerence and ﬁnite to the crosswind diﬀusion parameter (ck ) and
volume methods. SUPG can be shown the anisotropic diﬀusion parameter (δsd ). This
to add a smaller amount of stability analysis was performed for a cavity with A =
than anisotropic diﬀusion but yields a 15. The results are shown in Figure 2-3.
higher accuracy than anisotropic diﬀu- For the anisotropic diﬀusion Figure 2 (a-
sion. SUPG is less stabilizing than GLS c) shows the changes in the streamfunction as
but is also computationally less expen- the tuning parameter δsd was changed. In this
(a) δsd = 0.05 (b) δsd = 0.15 (c) δsd = 0.25 (d) ck = 0.01 (e) ck = 0.05 (f) ck = 0.1
Fig. 2: (a-c) Eﬀect of variation of anisotropic diﬀusion for Ra = 1.4e5, A = 15 and ck = 0.01 on the
streamfunction at t = 100. (d-f) Eﬀect of variation of crosswind diﬀusion for Ra = 1.4e5, A = 15
and δsd = 0.05 on the streamfunction at t = 100.
case the structure of the main cell is best de- properties found for A = 15 did not give good
ﬁned for the lower values of δsd . Figure 3 docu- results for this aspect ratio, and no critical
ments the temperature oscillations that occur Rayleigh number could be determined. The
near the intersection of the cells in the upper solution remained time invariant throughout
half. the range of Ra considered. The streamfunc-
For the crosswind diﬀusion Figure 2 (d-f) tion and temperature contours for this aspect
shows the changes in the streamfunction as the ratio appear to be consistent with prior bench-
ck parameter was adjusted. In this case only mark studies.
slight changes in the streamfunction are visi- The numerical diﬀusion parameters that
ble, primarily near the upper and lower cells. give appropriate results for higher aspect ra-
As shown in Figure 3 as ck is increased oscilla- tios (A > 15) may not be appropriate for
tions in the cells are halted completely. Over this aspect ratio. It is also possible based on
the entire range of Ra examined the higher ck the summary map provided by Chenoweth and
values prevented the transition to oscillatory Paolucci (Figure 1) that the critical Rayleigh
convection. number is aﬀected by the temperature diﬀer-
Based on this analysis a low value of cross- ences examined in this work. Small diﬀerences
wind diﬀusion (ck = 0.01) and a moderate in computational techniques may also be im-
value of anisotropic diﬀusion (δsd = 0.15) were pacting this aspect ratio.
chosen for the analysis of larger aspect ratios.
5.3 High Aspect Ratios (A > 15)
5.2 Small Aspect Ratio (A=8)
In all the simulations for the higher aspect ra-
Analysis of a small aspect ratio cavity was per- tios the streamfunctions and temperature con-
formed because the benchmark of Xin  was tours are consistent with those reported in the
very complete. Using the numerical diﬀusion literature. The streamfunctions are shown in
Fig. 3: Eﬀect of variation of diﬀusion for Ra = 2.5e5, A = 15 on the air temperature located at 2/3
of the full cavity height (upper half). a) Diﬀusion where ck = 0.01. i) δsd = 0.05 ii) δsd = 0.15 iii)
δsd = 0.25. b) Diﬀusion where δsd = 0.15. i) ck = 0.01 ii) ck = 0.05 iii) ck = 0.1.
Figure 4 and the temperature contours are
shown in Figure 5. Author Racritical Wave
Most authors observe a cat-eye pattern in Number
ﬂow for high aspect ratios (A > 12) that then A=8
transitions to a traveling wave instability. In Xin and Le Quere  3.1e5 1.7
this work, the cat-eye pattern occurred as ex- Present work - -
pected and became traveling waves for A = 33.
A = 15
Analysis of the cavity for A = 15 over the Liakopoulos et al.  1.4e5 -
range of Rayleigh numbers shows that Rac is Present work 1.4e5 2.62
near 1e5. This is very consistent with the ﬁnd-
ings of other numerical work . A qualita- A = 20
tive comparison of streamfunctions also aligns Lee and Korpela  1.1e4 2.82
well with the results of Chenoweth  and Li- Liakopoulos et al.  7.1e3 -
akopoulos . Vest and Arpaci  3.7e5 3.5
Present work 3.2e4 -
For A = 20 the literature reports Rac in
the range of 1e4 for the onset of multi-cellular
A = 33
motion. A stability transition in this regime
Reeve  6.7e3 2.77
was not calculated but the secondary motion
Vest and Arpaci  6.2e3 2.74
transition was found at Ra = 3.2e4 which is
Present work 5.8e3 2.49
in good agreement with the stability map of
Chenoweth . None of the comparison work
for A = 20 was based on experimental values Table 2: Summary of Cartesian numerical
for air (Vest and Arpaci report A = 20 results convection studies.
for oil), and other authors including Liakopou-
los did not have clear computational results for For A = 33 the critical Rayleigh number
Rac at this aspect ratio. was determined to be 5.8e3 for the transition
to multi-cellular ﬂow, slightly below the exper- air layer with large horizontal tempera-
imental value determined by Vest and Arpaci. ture diﬀerences, Journal of Fluid Mechan-
This value is also consistent with other nu- ics (1996).
merical studies. As the Rayleigh number was
 M.A. Christon, P.M. Gresho, and S.B.
increased the period of the waves in the simu-
Sutton, Computational predictability of
lation decreased. The critical wavenumber for
time-dependent natural convection ﬂows
this computational work agrees well with the
in enclosures (including a benchmark so-
experimental work of Vest and Arpaci  but
lution), International Journal for Numer-
improvements may be possible.
ical Methods in Fluids (2002).
Table 2 provides a summary of all the crit-
ical Rayleigh numbers determined by prior au-  G. de Vahl Davis, Natural convection of
thors and the present work. air in a square cavity: A bench mark
numerical solution, International Journal
for Numerical Methods in Fluids (1983).
 J. Elder, Laminar free convection in a
The comparison with experimental work indi- vertical slot, Journal of Fluid Mechanics
cates that by tuning the numerical diﬀusion (1965).
it is possible to match the values the critical
 G. Facas, On the convective instability in
Rayleigh numbers and wave lengths for high
long inclined cavities, ASME Stability of
aspect ratio cavities with natural convection.
Convective Flows (1992).
The critical Rayleigh numbers agree well with
the experimental results of Vest and Arpaci  P. Haldenwang and G. Labrosse, 2-d and
 for A = 33. A further analysis of the eﬀect 3-d spectral chebysheve solutions for free
of this numerical diﬀusion may be required to convection at high rayleigh number, Sixth
understand if the low aspect ratio case (A=8) International Symposium on Finite Ele-
can give better agreement. ment Methods in Flow Problems (1986).
In general, COMSOL was found to give
appropriate qualitative stream functions and  S.A. Korpela, D. Gozum, and C.B. Baxi,
temperature contours for all aspect ratios. On the stability of the conduction regime
The critical Rayleigh numbers for the sim- of natural convection in a vertical slot, In-
ulations decrease as the aspect ratio is in- ternational Journal Heat and Mass Trans-
creased which is consistent with prior work fer (1973).
in the literature. That the values of δsd and  Y. Lee and S. A. Korpela, Multicellular
ck that were found to be optimum for A=15 natural convection in a vertical slot, Jour-
were also appropriate for the higher aspect ra- nal of Fluid Mechanics (1983).
tios provides additional assurance to users of
COMSOL Multiphysics.  A. Liakopoulos, P.A. Blythe, and P.G.
Simpkins, Convective ﬂows in tall cavi-
ties, Simulation and Numerical Methods
References in Heat Transfer (1990).
 Comsol user’s guide.  S. Paolucci and D.R. Chenoweth, Tran-
sition to chaos in a diﬀerentially heated
 R.F. Bergholz, Instability of steady nat- vertical cavity, Journal of Fluid Mechan-
ural convection in a vertical ﬂuid layer, ics (1989).
Journal of Fluid Mechanics (1978).  P. Le Quere, A note on multiple and un-
steady solutions in two-dimensional con-
 A. Chait and S.A. Korpela, The sec- vection in a tall cavity, Journal of Heat
ondary ﬂow and its stability for natu- Transfer (1990).
ral convection in a tall vertical enclosure,
Journal of Fluid Mechanics (1989).  Hayden Reeve, Eﬀect of natural convec-
tion heat transfer during polymer optical
 D.R. Chenoweth and S. Paolucci, Nat- ﬁber drawing, Ph.D. thesis, University of
ural convection in an enclosed vertical Washington, 2003.
 S.A. Suslov and S. Paolucci, Stability of chebyshev pseudo-spectral benchmark for
natural convection ﬂow in a tall verti- the 8:1 diﬀerentially heated cavity, Inter-
cal enclosure under non-boussinesq condi- national Journal for Numerical Methods
tions, International Journal of Heat and in Fluids (2002).
Mass Transfer (1995).
 C.H. Vest and V.S. Arpaci, Stability of
natural convection in a vertical slot, In-
ternational Journal of Fluid Mechanics
These results were obtained as part of the
research supported by the National Science
 S. Xin and P. Le Quere, An extended Foundation through Grant 0626533.
(a) A=8 (b) A=15 (c) A=20 (d) A=33
Fig. 4: Streamfunction for Ra = 1.6e5 at t = 100.
(a) A=8 (b) A=15 (c) A=20 (d) A=33
Fig. 5: Temperature contours for Ra = 1.6e5 at t = 100.