Special Session

Computational Predictability Of
Natural Convection Flows in Enclosures

First M.I.T. Conference on
Computational Fluid and Solid Mechanics

Massachusetts Institute of Technology Cambridge, Massachusetts U.S.A.

June 12-14, 2001

proposed by:

Mark A. Christon¹

Livermore Software Technology Corporation

University of New Mexico

Philip M. Gresho

Steven B. Sutton

Lawrence Livermore National Laboratory

August 14, 2000

Download this document: pdf, or gzip'd postscript.

Introduction

Modern thermal design practices often rely on a `predictive' simulation capability - although predictability is rarely quantified and often difficult to confidently achieve in practice. The computational predictability of natural convection in enclosures is a significant issue for many industrial thermal design problems. One example of this is the design for mitigation of optical distortion due to buoyancy-driven flow in large-scale laser systems.

In many instances the sensitivity of buoyancy-driven enclosure flows can be linked to the presence of multiple bifurcation points that yield laminar thermal convective processes that transition from steady to various modes of unsteady flow [6]. This behavior is brought to light by a problem as `simple' as a differentially-heated tall rectangular cavity (8:1 height/width aspect ratio) filled with a Boussinesq fluid with $Pr=0.71$ - which defines, at least partially, the focus of this special session. For our purposes, the differentially-heated cavity provides a virtual fluid dynamics laboratory as pointed out by Le Quéré [4]:

``In conclusion let us emphasize that the differentially-heated cavity, in addition to its relevance as a model of convective heat transfer, turns out to be a real fluid mechanics laboratory in itself. The spatial structure of the flow is made of vertical and horizontal boundary layers, of corner structures, of a stratified core ... which depend very sensitively on the aspect ratio, Prandtl number and thermal boundary conditions (even a fly-wheel structure can be found at low $Pr$). All these features cooperate to give rise to very complex time behaviors resulting from several instability mechanisms, traveling waves in the vertical boundary layers, thermal instabilities along the horizontal walls in particular, which can interact strongly with internal wave dynamics.''

In the 8:1 cavity, the spectrum of the Jacobian of the Navier-Stokes equations about a steady-state solution is characterized by an infinite number of eigenvalues, either real or complex conjugates. For increasing Rayleigh number, some of the eigenvalues can cross the imaginary axis indicating bifurcation points (i.e., steady flow becomes unsteady, á la Hopf). Preliminary computations in the air filled 8:1 cavity (Xin and Le Quéré[7]) have indicated that two pairs of complex conjugate eigenvalues cross the imaginary axis in the vicinity of $Ra = 3.1 \times 10^5$. One of the corresponding eigenmodes has the skew-symmetry property of the base flow, while the other, the first unstable mode, is not skew-symmetric. This suggests that there will be significant sensitivity to the choice of boundary and initial conditions. That is, the choice of skew-symmetric conditions can promote the saturation of the second unstable mode which is skew-symmetric. In contrast, a random perturbation of the temperature field around the mean can promote the growth of the first unstable mode which is not skew-symmetric - at least for a finite period of time. Due, in part, to the presence of multiple unstable modes with a relatively small separation in $Ra$, the apparently simple differentially-heated cavity problem is not as simple as one might initially believe.

Additionally, the simulation of this buoyancy-driven flow is remarkably susceptible to the deleterious effects of numerical damping and/or dispersion introduced by commonly used `tricks-of-the-trade', thus making it surprisingly challenging. For example, numerical tests have demonstrated that the damping/dispersion artifacts from the simplest time-marching advection treatment with balancing tensor diffusivity can destroy the delicate thermal convective processes present in this enclosure when close to the critical Rayleigh number. In fact, coarse-grid computations may exhibit steady-state solutions even though the true solution is unsteady - requiring higher resolution grids than may be initially thought. Even when dissipation and dispersion have seemingly been minimized, computational experiments have shown that the amplitude of the periodic temperature oscillations can vary by as much as an order of magnitude depending on the specifics of the spatial discretization, grid resolution, stopping criteria for iterative solvers, and even the use of advective vs. conservative forms of the governing equations.

Ultimately the sensitivity of this class of flow problem to initial and boundary conditions, formulation details, and numerical procedures raises at least the following questions: What is the critical Rayleigh number, above which the flow will be unsteady, for the 8:1 enclosure? What is the behavior of the flow field at Rayleigh numbers slightly above critical? What is the role of linear stability analyses in predicting unstable modes? Can non-linear dynamics provide any insight into the behavior of the 8:1 cavity? What can be said about the relationship between (unstable) steady-state solutions and time-averaged periodic solutions? What is the best formulation and associated numerical procedure to use in order to ameliorate the sensitivities observed in practice and raise the level of accurate predictability? Which `other' numerical methods, i.e., discretization, time integrator, stabilization, preconditioned iterative technique, and `tricks-of-the-trade' are at least viable - and which are not?

In order to answer these questions and more, a special session is being organized for the First MIT Conference on Computational Fluid and Solid Mechanics. The session organizers are soliciting the contribution of solutions to the 8:1 differentially-heated cavity problem for near-critical Rayleigh numbers, and hope to see finite difference, finite volume, finite element, and spectral methods applied to this seemingly-simple 2-D problem. The application of commercial CFD codes is also highly encouraged.

Problem Definition

In this section, we present the 8:1 differentially-heated cavity as an initial boundary value problem with one set of initial conditions that may be used for a purely transient simulation. Here the emphasis is on a time-accurate computation with the fluid initially at rest. However, we recognize that a variety of analysis techniques ranging from steady-state to direct time integration and computation of unstable eigenmodes may be applied to this problem. We encourage the use of multiple analysis techniques, but leave the details of the steady-state and stability analyses to the contributors.

The 8:1 buoyancy driven enclosure flow problem is based upon the geometrical configuration shown in Figure 1 where $W$ is the width and $H$ the height of the enclosure. The enclosure aspect ratio is $A=H/W$ and takes on the value $A=8$. The gravity vector is directed in the negative $y$-coordinate direction, and the Boussinesq approximation for the buoyancy forces is assumed to be valid; i.e., only small temperature excursions from the mean temperature are admitted.

Figure 1: Differentially heated enclosure with 8:1 aspect ratio, insulated horizontal walls and constant temperature vertical walls.

The non-dimensional governing equations for the time-dependent thermal convection problem are the incompressible Navier-Stokes equations, conservation of mass, and the energy equation written in terms of temperature:

$$\frac{\partial \vec{u} }{\partial t} + \vec{u} \cdot \nabla ... ... P + \sqrt{\frac{Pr}{Ra}} \nabla^{2} \vec{u} + \hat{j} \theta,$$ (1)

$$\nabla \cdot \vec{u} = 0,$$ (2)

and

$$\frac{\partial \theta }{\partial t} + \vec{u} \cdot \nabla \theta = \frac{1}{\sqrt{Ra Pr}} \nabla^{2} \theta,$$ (3)

where $\vec{u} = (u,v)$, $P$ and $\theta$ are the velocity, the deviation from hydrostatic pressure, and temperature respectively, and $\hat{j}$ the unit vector in the $y$-direction. These non-dimensional equations were obtained using the characteristic length $W$, velocity $U=\sqrt{g \beta W \Delta T}$, time scale $\tau=W/U$, and pressure $\tilde{P}=\rho U^2$. Here, $\rho$ is the mass density, $g$ the gravitational acceleration, and $\beta$ the coefficient of thermal expansion. The non-dimensional temperature is defined in terms of the wall temperature difference and a reference temperature as

$$\theta = \frac{T - T_{r}}{T_h - T_c},$$ (4)

where

$$T_{r} = \frac{T_h + T_c}{2}.$$ (5)

$T_h$ is the prescribed temperature of the hot wall, and $T_c$ is the temperature of the cold wall.

The Rayleigh number is

$$Ra = \frac{g \beta \Delta T W^3}{\nu \alpha},$$ (6)

where $\alpha$ is the thermal diffusivity, $\nu$ the kinematic viscosity, and $\Delta T = T_h - T_c$ the temperature difference between the hot and cold walls. The Prandtl number is $Pr=\nu/\alpha$ and fixed at $Pr=0.71$. Additional information relating to this problem may be found in de Vahl Davis [1], Le Quéré and deRoquefort [3], Le Quéré[4] and Le Quéré and Behnia[5].

Boundary Conditions

The enclosure boundary conditions are simple and consist of no-slip walls, insulated (zero heat flux) horizontal walls, and constant-temperature vertical walls. The no-slip and no-penetration conditions are prescribed as $u=v=0$ on all walls, the left-wall is held at a constant `hot' temperature, while the right-wall is held at a constant `cold' temperature. The non-dimensional boundary conditions are summarized in Table 1.

Table 1: Boundary conditions for the differentially-heated enclosure.

	Left Wall	Right Wall	Bottom Wall	Top Wall
	$x=0$	$x=W$	$y=0$	$y=H$
Velocity BC's	$u=v=0$	$u=v=0$	$u=v=0$	$u=v=0$
Thermal BC's	$\theta=+1/2$	$\theta=-1/2$	${\partial \theta}/{\partial y}=0$	${\partial \theta}/{\partial y}=0$

Initial Conditions

In this section, we describe one set of initial conditions that may be used for a transient simulation. Here, the fluid is isothermal and initially at rest:

$$\vec{u}(\vec{x},0) = 0,$$ (7)

and

$$\theta(\vec{x},0) = 0.$$ (8)

We encourage the use of alternative initial conditions as a means to test the sensitivity of this problem and to further probe the space of solutions. For example, a random perturbation of the initial constant temperature is also an acceptable initial condition - as are any that are not skew-symmetric.

Paper Submissions

Contributions to the special session will be handled as follows. First, all submissions to the special session should follow the format described on the conference web page - see http://www.firstmitconference.org/Call.html.

Due to limited space for the publications, the authors should not repeat the introductory material or problem definition in their submissions. Instead, authors are asked to simply cite this document for the problem definition and present only their compulsory results as discussed below.

Submissions for the special session are due by October 15, 2000 and should be sent directly to:

           Mark A. Christon
           LSTC
           7209 Aztec Rd., NE
           Albuquerque, NM 87110

Authors will be informed prior to December 1, 2000 regarding acceptance of their contribution. In addition to the conference proceedings, we anticipate that the final results presented at the special session will be published in summary form in the International Journal for Numerical Methods in Fluids.

We realize that complete results for this problem may not be available by October 15th, and note that submissions with only preliminary results - particularly all or most of the tabulated results outlined below - along with a description of anticipated additional results are acceptable.

The following sections detail the parameter space and components of an `ideal' submission and include both compulsory and optional components.

Compulsory Results

This section outlines the quantities of primary interest that are to be reported, and suggests some additional quantities that contributors may wish to report. The compulsory data should be viewed as highly desirable for the sake of performing meaningful comparisons rather than an absolute requirement for paper submission.

The compulsory results for the special session are to be presented in 3 tables and two plots as outlined below. In the ensuing description, the compulsory data for the transient and steady-state computations are categorized according to the date type, i.e., point, wall and global data. We encourage and hope to see results from steady-state, transient, and stability computations.

Parameter Space

In this problem, the primary physical parameters consist of the Prandtl number $Pr=0.71$ and a super-critical Rayleigh number which is fixed at $Ra=3.4 \times 10^5$. Due to the amount of data that is required for this problem and the fact that computations on multiple grids will be necessary, only one Rayleigh number is being specified. There is, of course, no restriction on submission of results for additional Rayleigh numbers and we encourage participants to include their best prediction of the critical Rayleigh number, $Ra_{crit}$: the transition from steady to time-periodic behavior.

Grid Resolution Guidelines

There are no strict requirements on grid resolution, however, all participants are asked to produce results that they believe are sufficiently accurate. In order to provide some guidance on grid resolution, we report here some preliminary results. Using a second-order method, a smoothly graded `coarse-grid' resolution consisting of $27 \times 121$ grid points has been used to successfully compute time-dependent results for $Ra=3.4 \times 10^5$. Doubling the grid resolution, e.g., $53 \times 241$ grid points, has shown that there is still significant sensitivity in (at least) the amplitude of temperature oscillations. Interestingly, the amplitude of the temperature oscillations is also dependent on the specifics of the Navier-Stokes formulation.

Some guidance on the vertical-wall boundary layer thickness may be found in Gill[2]. Note that graded meshes are suggested since the resolution for a uniform grid may make the computations prohibitively time consuming. We suggest using grids with approximately a 1:5 x-to-y ratio of grid points starting with a coarse grid of $21 \times 101$. Increasing grid resolution should be obtained by grid doubling, e.g., a medium grid of $41 \times 201$ and a fine grid of $81 \times 401$.

Point Data

A series of physical data should be recorded at the compulsory time-history points from Table 2. The coordinates that are identified in Table 2 are non-dimensional and shown in Figure 1.

Table 2: Non-dimensional coordinates of time-history points.

point	x-coordinate	y-coordinate
1	$0.1810$	$7.3700$
2	$0.8190$	$0.6300$
3	$0.1810$	$0.6300$
4	$0.8190$	$7.3700$
5	$0.1810$	$4.0000$

The time history data for point 1 should be tabulated as shown in Table 3. In addition, a plot of the oscillatory variation in temperature at point 1, $\theta_1$, should accompany the tabulated results. A second plot showing the skewness, $\epsilon_{12}$, should be included only if the skewness is found to be non-zero during the periodic phase. For all time-dependent computations, the average value and oscillation amplitude are to be reported, along with the period of oscillation. The average should be taken over one or more complete periods where the amplitude and period should be constant. The use of an FFT to obtain the amplitude and period is encouraged, but not mandatory. Where steady-state computations are concerned, only the point-values are required and may be tabulated in a similar format.

The computation of the cycle average is based upon achieving a statistically stationary state where the period and amplitude are constant. For a generic variable, $\phi=u,v,\theta,\epsilon,...$, the average may be computed as

$$\overline{\phi} = \frac{1}{T} \int_t^{t+T} \phi(\vec{x},t) dt$$ (9)

where $T$ represents the period of time for which the average is computed. This computation implies that any starting transients have completed thereby permitting the average to be based on a stable oscillatory solution.

Using the average value, perturbation fields may be computed as

$$\delta \phi(x,t) = \phi(x,t) - \overline{\phi}.$$ (10)

The vorticity is defined as

$$\omega = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y},$$ (11)

and the stream function as

$$u=\frac{\partial \psi}{\partial y},$$ (12)

and

$$v=-\frac{\partial \psi}{\partial x},$$ (13)

with $\psi=0$ on the walls.

In Table 3, the skewness provides a measure of the skew-symmetry of the temperature field. Using time-history points 1 and 2, the skewness is computed as

$$\epsilon_{12} = \theta_1 + \theta_2.$$ (14)

This parameter should be zero for a skew-symmetric temperature field, i.e., $\theta_1 = - \theta_2$. If the skewness is non-zero, please provide a time-history plot with the tabulated results.

The pressure differences in Table 3 are defined as

$$\Delta P_{ij} = P_i - P_j$$ (15)

where $i$ and $j$ indicate the time-history points used to compute the pressure difference.

Table 3: Tabulated results should identify the grid resolution ($Nx \times Ny$), time history point, time duration used to obtain the average ($T$), amplitude and period, and the number of time steps per period, $Nper$, in the computation.

Increasing Grid Resolution $\longrightarrow$
	Grid Resolution: $Nx_1 \times Ny_1$			Grid Resolution: $Nx_2 \times Ny_2$
	Time Duration: $T_1$			Time Duration: $T_2$
	Steps per Period: $N_{per_1}$			Steps per Period: $N_{per_2}$
Quantity	Average	Amplitude	Period	Average	Amplitude	Period
X-Velocity	$\overline{u}$	$\cal{U}$	${\cal T}_u$	$\overline{u}$	...
Y-Velocity	$\overline{v}$	$\cal{V}$	${\cal T}_v$	$\vdots$
Temperature	$\overline{\theta}$	$\Theta$	${\cal T}_{\theta}$
Skewness	$\overline{\epsilon}_{12}$	$\Xi_{12}$	${\cal T}_{\epsilon}$
Stream Function	$\overline{\psi}$	$\Psi$	${\cal T}_{\psi}$
Vorticity	$\overline{\omega}$	$\Omega$	${\cal T}_{\omega}$
$\Delta P_{14}$	$\overline{\Delta P}_{14}$	$\Delta P_{14}$	${\cal T}_P$
$\Delta P_{51}$	$\overline{\Delta P}_{51}$	$\Delta P_{51}$	${\cal T}_P$
$\Delta P_{35}$	$\overline{\Delta P}_{35}$	$\Delta P_{35}$	${\cal T}_P$

Wall Data

The wall Nusselt numbers are also a component of the compulsory results and are defined as

$$Nu(t) = \frac{1}{H} \int_0^H \left . \frac{\partial \theta}{\partial x} \right \vert _{x=0,W} dy.$$ (16)

The evaluation of the Nusselt number should be performed for each wall, i.e., $x=0$ and $x=W$. This data should be tabulated as shown in Table 4. For all time-dependent computations, the average, and the oscillation amplitude are to be reported along with the period of oscillation. Here, the average is taken to be over one or more complete periods of oscillation where the amplitude and period are constant. Again, the use of an FFT to obtain the amplitude and period is encouraged, but not mandatory. Where steady-state computations are concerned, the Nusselt number for each wall should be reported.

Table 4: Tabulated results should identify the grid resolution ($Nx \times Ny$), time history point, time duration used to obtain the average ($T$), amplitude and period, and the number of time steps per period, $Nper$, in the computation.

Increasing Grid Resolution $\longrightarrow$
	Grid Resolution: $Nx_1 \times Ny_1$			Grid Resolution: $Nx_2 \times Ny_2$
	Time Duration: $T_1$			Time Duration: $T_2$
	Steps per Period: $N_{per}$			Steps per Period: $N_{per}$
Quantity	Average	Amplitude	Period	Average	Amplitude	Period
Nusselt ($x=0$)	$\overline{Nu}$	$Nu$	${\cal T}_{Nu}$	$\overline{Nu}$	...
Nusselt ($x=W$)	$\overline{Nu}$	$Nu$	${\cal T}_{Nu}$	$\vdots$

Global Data

In addition to the point and wall data, a number of relevant global quantities are required as delineated in Table 5. For our purposes, the kinetic energy and enstrophy provide useful metrics; e.g., the square-root of kinetic energy provides a measure of the average velocity in the enclosure. That is, our velocity metric is

$${\hat u}(t)=\sqrt{\frac{1}{2A} \int_A \vec{u} \cdot \vec{u} \ dA},$$ (17)

where $A$ is the area of the enclosure, $A=W \times H$. Similarly, the measure of the average vorticity is based on the enstrophy,

$${\hat \omega}(t)=\sqrt{\frac{1}{2A} \int_A \omega^2 \ dA}.$$ (18)

Table 5: Tabulated results should identify the grid resolution ($Nx \times Ny$), time history point, time duration used to obtain the average ($T$), amplitude and period, and the number of time steps per period, $Nper$, in the computation.

Increasing Grid Resolution $\longrightarrow$
	Grid Resolution: $Nx_1 \times Ny_1$			Grid Resolution: $Nx_2 \times Ny_2$
	Time Duration: $T_1$			Time Duration: $T_2$
	Steps per Period: $N_{per}$			Steps per Period: $N_{per}$
Quantity	Average	Amplitude	Period	Average	Amplitude	Period
${\hat u}$	$\overline{\hat u}$	${\hat {\cal U}}$	${\cal T}_{\hat u}$	$\overline{\hat u}$	...
${\hat \omega}$	$\overline{\hat \omega}$	${\hat\Omega}$	${\cal T}_{\hat \omega}$	$\overline{\hat \omega}$	...

Methodology

A summary of the method used to solve the problem is required and should provide a concise description of the following items:

• Spatial discretization method and a description of the grids used.
• Solution procedure, e.g., time integration procedure, linear solvers, etc.
• Stopping criteria used for all iterative procedures.
• Description of the advective treatment - including the use of artificial viscosity, limiters and blended first/second order methods.
• Description of any `stabilizing' terms in the formulation or solution methodology, e.g., Galerkin Least-Squares (GLS) stabilization.
• An estimate of the order of accuracy of the method and, if possible, the accuracy of the submitted results.
• Any new or unique features of the method used to solve the differentially-heated cavity problem.

Computational Resources

A summary of the computational resources is required and should provide a concise description of the following items:

• Machine, e.g., Compaq, DEC, SGI, etc.,
• Clock rate in MegaHertz [MHz],
• Total memory MegaBytes [MBytes],
• peak FLOP (floating point operation) rate [MFLOPs] and/or the specfp95 rating - see: http://www.spec.org/osg/cpu95/ ,
• Number of processors used and the number of grid points per processor (degree-of-granularity) for parallel computations.
• CPU (or wall clock time for parallel) used per grid point per time step [micro-seconds/point/step], and
• Memory used (per processor in parallel) [MBytes].

Optional Results

Extensions of the results beyond the compulsory data is encouraged. However, due to space limitations, the optional results should not be reported in the paper submission. Several possibilities for optional results include:

• Compute the the period and amplitude of the velocity and temperature oscillation for supercritical Rayleigh numbers in the range $3.1 \times 10^5 \leq Ra \leq 6.2 \times 10^5$, e.g., for $Ra=3.2 \times 10^5$.

• Estimate the the critical Rayleigh number for the 8:1 enclosure based on your computed results.

• The point data may be reported at any/all of the time history points shown in Figure 1.

• As a check report the skewness based on points 3 and 4, $\epsilon_{34}$. If the skewness is non-zero, plot it as a function of time.

• Plot the skewness metrics, $\epsilon_{12}$ and $\epsilon_{34}$, for the entire computation.

• Report the maximum stream function value and location over one oscillation cycle.

• Report minimum and maximum vorticity values and locations over one oscillation cycle.

• Plot the vertical force for each constant-temperature wall over one cycle. Report the average over the cycle.

• Provide contour plots (color not required) of the instantaneous perturbation velocity, temperature, and pressure fields per Eq.(10). Plots should be provided at 4 equally-spaced points in time during an oscillation cycle.

Important Dates

October 15, 2000	Deadline for submissions
December 15, 2001	Notification of acceptance
June 12-14, 2001	First M.I.T. Conference

Organizers & E-mail Addresses

Questions regarding the special session may be sent via e-mail to any of the organizers.

Mark A. Christon	Philip M. Gresho	Steven B. Sutton
LSTC/Univ. of New Mexico	LLNL	LLNL
e-mail: christon@lstc.com	e-mail: pgresho@llnl.gov	e-mail: sutton4@llnl.gov
Ph: (505) 875-0746	Ph: (925) 422-1812	Ph: (925) 422-0322
Fx: (505) 875-0430	Fx: (925) 423-5167	Fx: (925) 423-6195

Bibliography

1

G. DE VAHL DAVIS, Finite difference methods for natural and mixed convection in enclosures, in 8th International Heat Transfer Conference, ASME, August 1986, pp. 101-109.

2

A. E. GILL, The boundary-layer regime for convection in a rectangular cavity, Journal of Fluid Mechanics, 26 (1966), pp. 515-536.

3

P. LEQU´ER´E AND T. A. DEROQUEFORT, Transition to unsteady natural convection of air in vertical differentially heated cavities: Influence of thermal boundary conditions on the horizontal walls, in 8th International Heat Transfer Conference, ASME, August 1986.

4

P. L. QU´ER´E, Onset of unsteadiness, routes to chaos and simulations of chaotic flows in cavities heated from the side: a review of present status, in Proceedings of the Tenth International Heat Transfer Conference, G. F. Hewitt, ed., Brighton, UK, 1994, pp. 281-296.

5

P. L. QU´ER´E AND M. BEHNIA, From onset of unsteadiness to chaos in a differentially heated cavity, Journal of Fluid Mechanics, 359 (1998), pp. 81-107.

6

K. H. WINTERS, Hopf bifurcation in the double-glazing problem with conducting boundaryes, Journal of Heat Transfer, 109 (1988).

7

S. XIN AND P. L. QU´ER´E, personal communication, July 2000.

Web Counter Page hits:

About this document ...

This document was generated using the LaTeX2HTML translator Version 99.1 release (March 30, 1999)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -split 0 mit_convection.tex

The translation was initiated by Mark A. Christon on 2000-08-14

Footnotes

... Christon¹

Organizers are listed in alphabetical order.