- Email: [email protected]
Development of a robust finite element CFD procedure for predicting indoor room air motion
Building and Environment, Vol. 29, No. 3, pp. 261 273, 1994
~
Copyright © 1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0360-1323/94 $7.00 + 0.00
Pergamon 0360-1323(94)E0008-F
Development of a Robust Finite Element CFD Procedure for Predicting Indoor Room Air Motion A.J.BAKER* P. T. WILLIAMS* R. M. KELSOt The requirement is verifiable high quality computational (CFD) solutions to Reynolds-averaged Navier-Stokes equations for buoyant and turbulent flows in genuine three-dimensional room geometries. Current-practice theory~codes stabilize intrinsic dispersive error via artificial diffusion methods, with resultant compromise of genuine physical~turbulent diffusion versus convective processes. A Taylor weak statement CFD theory, and finite element code implementation, represents a potential significant accuracy/stability improvement for three-dimensional buoyant forced-, mixed- and natural-convection room airflow prediction. This theory, documentary validations, and quantitative data comparisons from conducted CFD experiments are presented.
INTRODUCTION
significant number of assumptions, it is unreasonable to expect the results of a given CFD simulation to be accurate unless (1) the problem definition is trivial, or (2) great care is exercised in design and execution of the "CFD experiment." In truth, CFD simulation is a new "experimental protocol" that under parametric variation produces a flow field data base for interpretation. Engineers have for years conducted physical experiments in the "real laboratory," where a carefully-constructed scale model is interrogated using data-measuring devices carefully selected to control error in the acquisition process. A CFD simulation experiment is analogously conducted in the " C F D laboratory," i.e. via a math model on a computer, rather than a physical model, where careful attention is (must be) paid to sources of error that can compromise the CFD data generation process. To answer the posed questions then, one must seek out the error mechanisms existing to compromise the fidelity of a CFD room air motion simulation. The distribution of mass, velocity and temperature in a given room certainly satisfies the basic conservation principles we all studied in college. Expressed in partial differential equation form, these are called the "Navier-Stokes" equations, which mathematically describe in differential calculus form the interaction among the dependent variables density (p), pressure (p), velocity (vector) field (u) and temperature (T). However, most room air flow fields are at least locally "turbulent," and even though the NavierStokes (NS) system remains valid, a "Reynolds-averaging" manipulation (cf. [3]) is required to be applied for CFD computational tractibility. Any such resolution of velocity, into mean and fluctuating (about the mean) components, involves several assumptions, and usually leads to definition of a turbulent (Reynolds) stress-strain rate constitutive hypothesis containing a turbulent "eddy viscosity." Such an expression is similar-appearing to
COMPUTATIONAL fluid dynamics (CFD) is the maturing technical discipline of development of math/physics models amenable "to computing," that are predictive regarding the motion of fluids through and around practical geometries. Indoor room air motion is the problem class of interest, and CFD is becoming applied with increasing frequency to simulation of the buoyant and perhaps-turbulent flow of air in pertinent geometries. Nielsen [1] pioneered in applying CFD to room air motion simulation, and in the ensuing two decades, progressively more indoor air applications of CFD are becoming published. These activities have focused mainly in Europe and Japan, and for example led to formation of Annex 20 [2], under the auspices of the International Energy Agency (lEA), with the goal to organize a basis for validation of CFD methods/codes for application to indoor room air motion prediction. The fundamental question to be raised regarding CFD methodology applied to room air motion prediction is, "Just how good is it, or can it be?" The companion query must be, "What issues exist to compromise the ability of a CFD method to accurately simulate a (any) given room air motion situation?" The answers to these questions are intrinsically important to the indoor air technical community. This paper addresses them, first philosophically, and then technically with development and benchmark validation of a new CFD methodology designed to alleviate limitations in the current practice. The fundamental precept to seeking the answers is that, since a//CFD procedures constitute the union of a *Department of Engineering Science and Mechanics and tSchool of Architecture, The University of Tennessee, Knoxville, TN 37996, U.S.A. 261
A . J . Baker et al.
262
Stokes' laminar viscosity constitutive law (contained in the NS system). The distinction is that in NS the genuine fluid viscosity v depends on the fluid itself (and weakly on temperature), while in Reynolds-averaged Navier Stokes (RaNS), the turbulent eddy viscosity v' depends upon the state of the flow, i.e. v' = vt(u). Therefore, a principle fundamental assumption for C F D simulation of room air motion is the form selected for the governing Reynolds-averaged Navier-Stokes PDE system. The almost-universal selection has been to append the two-equation turbulent kinetic energy (TKE) transport PDE system to RaNS, which yields a closed, mathematically well-posed system. The drawback is that the TKE model assumes the flow is.[ully turbulent, hence predicts an excessive dose of vt for the subtly-turbulent flows existing throughout indoor rooms away from HVAC supply systems and obstructions with edges. One measure of "how turbulent" a flow (simulation) is, is the turbulence Reynolds number (Re') defined as
yl V
Re' =- -- >/- O.
(1)
For a laminar flow, Re' = 0 by definition, while Re' order (102-103 ) for fully turbulent flows. Weakly turbulent flows are characterized (roughly) by 0 < Re' ~< 30 while the transition region to full turbulence occurs (essentially) within 30 ~< Re' ~< 100. The TKE fully-turbulent closure model has been adapted with "low turbulence" adjustments to account for wall layer effects, but these constructions involve many more assumptions, and remain to be validated for general application to weakly turbulent central flows away from walls, supply diffusers, etc. Philosophically then, the answer to the premier question of "how good C F D may be" for room air motion simulation rests squarely on the validity of the assumed form of the RaNS governing PDE system, including all boundary and initial conditions. Errors committed here are not assessable until more exact PDE forms are derived and validated. This arena is a "hot topic" in turbulence research, but in practice today, it is unreasonable to expect that any given C F D simulation can be "exact." Hence, one must conduct a range o f " C F D experiments" to validate C F D simulation sensitivity to parameters embedded in the RaNS turbulence closure system. An example of such a C F D experiment is presented in this paper. One must next turn to the very fundamental questions, i.e. "How good is the C F D theory/code in emulating the genuine solution to the selected RaNS system?" All C F D theories, be they finite difference (FD), finite volume (FV), finite element (FE), or any other form, can do no more than establish an approximation to the genuine RaNS solution. One fundamental flaw in any C F D theory is that it employs a mesh to support construction of the approximate solution process. Just as Reynoldsaveraging precludes obtaining the complete NS solution itself, a C F D approximation contains no resolution of solution information of two cell-wide wavelength. Figure la illustrates a I-D mesh of measure "Ax," and a C F D solution generates information only at the nodes (labelled j, j + 1, etc.). Obviously, any data of wavelength 2Ax is
a)
" j . 1~
b)
j+2 x 1;3 j~. /"-'1
Fig. 1. Sketches of spectral content resolution on a mesh of measure Ax. never seen(!) at the nodes, hence by the C F D solution process. Therefore, the C F D mesh must be sufficiently dense to admit capture of essential solution spectral content, or be totally worthless. Figure lb confirms that "5Ax" wavelength information is quite resolvable by a mesh of "measure" Ax. In practice, the C F D code execution process transmits information about the RaNS approximate solution over the nodes of the mesh. Therefore, one measure of how well this is accomplished, by competing F D / F V / F E theories, is to compute the algorithm intrinsic phase velocity error as a function of "data" of wavelength 2 = nAx, n/> 2. Using symbolic manipulation [4], this has now been unequivocally established, and Fig. 2 summarizes these data for six competing C F D theory implementations. The current practice F D and FV methodologies are theoretically confirmed to be of quite poorer performance in the short wavelength region (in this measure) in comparison to the emerging FE weak statement methods. This predictive theoretical analysis is documented further herein. One consequence of this phase error distortion, intrinsic to all C F D methods, is that the poorly transmitted information leaves behind short wavelength oscillations in the approximate solution. These "wiggles" can easily produce non-realizable data, e.g. negative temperatures, or can cascade (the RaNS PDE system is nonlinear) and cause iterative divergence of the C F D code solution process. The almost universal "CFD-fix" for this is to add numerical dissipation to the algorithm, to diffuse the wiggles. This process unequivocally introduces a numerical "artificial viscosity" coefficient va that acts in concert with the genuine fluid viscosity v and the turbulent viscosity v'. In some instances, va can be so large as to absolutely swamp the other two physics-generated diffusion mechanisms. The control of this strictly-CFD error mechanism is probably most critical to a favorable answer to the practical question, "what issues exist to compromise . . . " posed earlier. This paper highlights the theory, and presents finn validation data, regarding development of a new, mathematically robust and versatile C F D procedure for the room air motion problem class. It is based on the Taylor weak statement (TWS) theory [5], with generalizations [6,7], which has proven capable of recovering historical dissipative C F D methods, while predicting improved, more-accurate replacement methodology. This paper
Prediction of Indoor Room Air Motion 1.0
I
I
I
263
!
Central F D -,~--QUICK-3 F V o QUICK-5 FV --+--IAaearFE -n--.
0.8
• °° *" °° *'
...."""/I//
Quadratic FE -.A-.Cubic FE --~---
...., //,,t... /'; /
•"
i "S
.."
0.6
///
,,"
/
o.°
ra ! ; / i
..
/
i '.
..'
/
/ /
;
i
.°
,/
.-"
C=1/2 0=1/2
!
."
/
..."
8
"
ii
!!
,,
/
,"/"
i ;
#
• o°
0.4
,o
• °° °°°
0.2 °, -° °° °°
°° o°
°°
/,s
°°
J
/s ,° °L'JI
s A"
i,, S
/
°°"
,
o°°° /
s s.s" ..s"
/
/ • i° ~,'
s.J
I"
•
t¢ °
I'
..[3"
°° °° °°
Olb oo
16
8
4
3
Modal Wavelength (X----~x) Fig. 2. Fourier stability phase error comparisons for various CFD algorithms.
highlights the TWS theory, applied to RaNS, and its implementation via a finite element semi-discretization of a Galerkin weak statement. Thereafter, the implicit temporal semi-discretization admits time-accurate, unsteady predictions, and sparse matrix methods are employed to iteratively solve the resulting non-linear algebraic equation system. Theoretical and computational performance validations are presented for select 2-D and 3-D benchmark problems, and for two genuine 3-D room air motion experiments.
THE PROBLEM S T A T E M E N T
A statistical manipulation, c.f. [3], produces a RaNS PDE system governing incompressible buoyant room air flows. Following nondimensionalization, and recalling the turbulent Reynolds number (Re') of the turbulence closure model, the 3-D PDE system governing mass, momentum and energy transport processes, hence requiring solution, is
80
~ (
×
= 0
(2)
au,+~( dxj u~uj+ P'~°
£e(u,) = - ~
l (l+Re')t-~xj+-~xJ)-ArOS,g=O
Re BAE 29:3-B
(3)
+ Pr'Jdxjj-s°
= 0.
(4)
The dependent variable set in (2)-(4) contains the meanflow velocity vector u ~ ui(xj,t), the kinematic pressure P=P/Po and the potential temperature ® - ~ ( T Tmin)/(Trnax-Tmin). The Boussinesq buoyancy body-force approximation is made in (3), where 6ig is the gravity unit vector. The dimensionless groups appearing in (3) and (4) exert a controlling impact on solution character. The definitions for Reynolds, Grashoff, Prandtl and Archimedes numbers are
UL
Re=--,
Gr=
3)
g fl(tma
x - - Trnin ) Z 3 V2
P r - poCpV Ar = GrRe -2. =
~e(po) = ~
1 Re
L~(O) = --~ + ~jxj ujO
k
'
(5)
'
The familiar Reynolds number (Re) is based on userselected reference velocity (U) and length (L) scales and the fluid kinematic viscosity v. The Grashoff number (Gr) contains a reference temperature spread, the gravity constant and the fluid compressibility factor (fl), roughly equal to inverse absolute temperature. The Prandtl number (Pr) is based on the reference (constant) density (P0), specific heat (cp) and thermal conductivity (k). The
A.J. Baker et al.
264
remaining parameters in (3) and (4) are turbulent Reynolds number Re', defined previously in (1), and the turbulent Prandtl number Pf, which nominally is defined identical to Pr for air. The character of solutions to (2)-(4) depends on these non-D groups, the distribution of Re' and boundary conditions. Since Re > 0 and Re' ~> 0, isolated equations (3) and (4) are elliptic boundary-value, hence the corresponding members of the mean-flow state variable q (x, t) = {u~, u2, u3, O} r, and/or their normal derivatives, are required to be known all around the problem statement boundary. The solution process for these equations is thus classic initial-boundary value, up to enforcement of the continuity constraint that us must be divergencefree, i.e. the velocity field must also satisfy eq. (2). This differential constraint is only enforceable approximately, via one of several inexact CFD theories, as a modelled action for pressure P. A thorough theoretical study [6] verifies a robust procedure appends a Poisson variable ~b satisfying the boundary value problem
where superscript tilde denotes any approximation to the true (i.e. divergence-free) vector velocity u, and n is the time-step index. The parametric variable q~ is then manipulated [6] into a kinematic model of the action of pressure via P"+ ~ _~ P " - d?.+ l/OAt, where subscript "A" denotes approximation and P" is the genuine pressure (from the previous solution at time t,). The CFD approximation PA must then be replaced, upon time-step convergence, with the genuine pressure P,+ ~, which is the solution to
t?2P - -~2 (u~u/-5fl(P) = (~x~ + Ox,axj ~ ×
+
coU(q) = ;~ + ~_ (.t - f ' i ) - s = 0
(8)
~(q,) = V2q , - s(q) = 0.
(9)
In (8),Ji = J)(q) is the kinetic" flux vector, containing convection and pressure effects, while fj' is the dissipative flux vector embodiment of the turbulence closure model (Re') plus natural diffusion (Re, Pr). In (9), the Laplacian operates on q~ = {q6, p}r, and s = s(q) is the respective source dependence on the solution to (8). Any CFD algorithm generates (only) an approximation to the true solution to (8) and (9), via a denumerable set of decisions leading to an algebraic system amenable "to computing." As discussed, a dominating dispersive error mechanism is intrinsic to this process for (8), ultimately leading to the near-universal use of "artificial diffusion" mechanisms for stabilizing computations, e.g. " Q U I C K " differencing. The emerging encompassing Taylor weak statement theory, [5,6], extended to systems (8), yields creation of the PDE system replacement for (8) in the form ~caC(q) - L,a(q) _ flAt c?~x--~( A j A k ~ ) = 0 ,
(10)
where Aj ----0fj/0 q is the RaNS kinetic flux vector Jacobian (matrix). The replacement of (8) by (10) in the CFD implementation process contributes directly to a tensor invariant, highly phase selective numerical dissipation mechanism, as will be verified. Hence, upon replacement of (8) by (10), any approximation qN(x, t) to the true solution q(x, t) to (8) and (9) in the continuum is of the functional form N
1
q(x, t) ~ q"¢ (x, t) = ~ LFj(X)Qj (t), j=J
~e(1 +Re')
)
-- A r ~ xs• 6jo. = 0.
(7)
The third (spatial) derivatives on velocity in (7) are computationally intractable to form to an acceptable discrete order of accuracy. The theoretical resolution is now achieved via a Galerkin weak statement theory [8], which lowers required differentiability order to first. Since (6) and (7) are elliptic, surrounding boundary conditions are required for well-posedness. For (6), the sole known boundary conditions are homogeneous Neumann (normal derivative), while non-homogeneous Neumann conditions are derived for (7) using (3), c.f. [8]. Finally, a closure model for turbulent Reynolds number Re'(x, t) is required, which is constrained herein to a scalar factor, as appropriate for CFD experiment comparative purposes.
(11)
where tPj(x), 1 ~
(12)
be satisfied for any (all!) given known functions ~(x). If the theory-designer picks ~i(x) and ~i(x) identical, then a "Galerkin" (finite element, FE) weak statement results. Conversely, if one chooses q)~(x) to be constants, then a "finite volume" (FV) weak statement form is produced. To actually evaluate the integrals defined in (12), for genuine problems, a spatial semi-discretization facilitates the calculus operations, whereupon, for either an FE or FV selection, (12) always produces a matrix ordinary differential equation (ODE) system of the form [9]
WEAK STATEMENT CFD ALGORITHM
d{Q} TWS h = [ M ] ~ + {R(Q)} = {0}.
A CFD algorithm for (2)-(7) generates solutions parameterized by Re, Gr, Pr and Re'. The CFD theoryderived dependent variable set (the "state variable") is q(x, t) = {ul, uz, u~, ®; 4), P} T, which at the semi-colon partition satisfies the PDE systems
In (13), the superscript h notation emphasizes that a "discretization" ~h of the domain I'~ of (8)-(10) has been invoked, in achieving this matrix expression. Further therein, [M] and {R(Q)} are global (nodal) rank square
(13)
Prediction of Indoor Room Air Motion and column matrices, respectively, and {Q} = {Q(t)} is the (unknown) state variable approximation coefficient (set) Qj at the nodes of the mesh defined by fP. Any ODE algorithm uses (13) directly for time derivative evaluation. Selecting the one-step, 0 implicit Euler family for example, and substituting (13) as needed [91, produces the terminal matrix equation system "for computing" as
{FQ} = [M] {Q.+I - Q.} + At(O{R(Q)}.+I + ( 1 - 0 ) {R(Q)}.) = {0},
(14)
where subscript "n" remains time level notation, The discrete Galerkin weak statement (12) for (9) directly produces the companion matrix systems
{FQ~} --- [D] {QA}- {S(Q) } = {0},
(15)
where square matrix [D] is the discrete Laplacian, {Q~} is the nodal array of 4? or eh, and {S(Q)} contains dependence on qh via the solution {Q(t)} from (14). The generic solution process for (I 4) and (15) involves use of an appropriate linear algebra approximation to the Newton iteration algorithm
[JA C] .~60~P+ , .+ l = - {FQ}P+~
(16)
P
p+l {Q}.+I ,+i = {Q}.+ 2 {6Q}.+~,
(17)
265
the point of destablizing the quasi-Newton iterative process for (16) and (17). Hence, the central stability issue with a CFD algorithm is control of this dispersive error mechanism, such that a convergent iteration process can produce a "smooth" solution, free of "wiggles." To date, there are no viable methods to annihilate this mode, hence the remaining (universally-applied) "fix" is to dissipate the oscillations via an added numerical diffusion process. For a linear model of equation (8), the Fourier stability analysis employed to develop Fig. 1 is extended to computation of the amplification factor for a variety of CFD algorithms. The exact solution amplification factor is unity, i.e. no (artificial) damping is present, hence departures from unity, as a function of modal wavelength ~ = nAx, characterize algorithm artificial damping. Figure 3, from [4], compares the error in (unit) amplification factor for a range of FD, FV and linear basis FE CFD algorithms, the latter for various TWS-theory E-levels, which are user specifiable. The horizontal line at zero error confirms that the non-artifically diffused FD and FE algorithms (fl --- 0) have no intrinsic dissipation error. The QUICK3 FV algorithm is most dissipative, especially in the short wavelength region, while the current FE algorithm with (user-definable) modest fl yields a more phase-selective dissipation mechanism.
i=O
where [JAC] is the Jacobian of (14) and (15). The Newton Jacobian defined in (16) is readily formed using calculus operations, cf. [10]. Thereafter, the range of quasi-Newton iterative procedures use approximations thereto. For this paper, a variable-segregated, 3-D quasi-Newton Jacobian construction was iteratively cycled, using a diagonally-preconditioned GMRES algorithm [6], in concert with an incomplete-Cholesky preconditioned conjugate gradient (PCG) solver for (15). ACCURACY, STABILITY AND DISPERSION ERROR The weak statement theory (8)-(12) confirms that
required fundamental CFD algorithm decisions focus on selection of the trial space Wi (x), and the test space ¢i(x), for defining qU and the TWS N error constraint respectively. Under a spatial semi-discretization, these decisions become commensurate with the FD concept of algorithm order of accuracy, hence asymptotic convergence under mesh refinement. Linear basis FE and traditional FV algorithms are typically second-order accurate, and coarse meshes aggravate the intrinsic semidiscrete approximation phase error mechanisms; recall Fig. 1. In practice, algorithm stability is of much greater practical significance than convergence for genuine RaNS problems, since the available mesh density is always limited by computer resources. The practical issue is thereby control of the third-order dispersive error mechanism intrinsic to all discrete approximation methodologies. As discussed, any mesh of measure Ax in one-dimension is incapable of resolving solution information of wavelength 2Ax. The non-linearity of a RaNS equation system amplifies this 2Ax oscillatory dispersive error mode, often to
RESULTS AND DISCUSSION The benchmark problem of specific pertinence to prediction of room air motion is the close-coupled step wall diffuser, cf. [111, which models the essence of a supply outlet. Figure 4 summarizes the 2-D Galerkin linear basis FE WSh steady state laminar solution obtained at Re = 390 in terms of (a) velocity vector uh distribution, with experimental location of separation zone reattachment (xl/S) noted, and (b) pressure distribution. In proceeding to Re = 650, a top-wall secondary separation occurs (Fig. 4c), in qualitative agreement with the (3-D) experiment; however, the discretized Poisson equation pressure field is noticeably polluted by a dispersive error mode (Fig. 4d). Homogeneous Poisson (Laplace) solutions are always smooth, for regular data, hence these oscillations must be due to the non-homogeniety in (7) involving derivatives of the CFD discrete approximation uh. Switching to the TWS h dissipative algorithm, using = 0.1 does not measurably alter u h, but this "filtering" totally annihilates the pressure solution oscillations (Fig. 4e). The identical trend is recorded in the fully 3-D benchmark simulation, Fig. 5. Additional details are given in [61. Interest in natural convection in geometrically-complex enclosures, with restricted communication, has been stimulated by applications involving energy-efficient passive-solar buildings, cold-storage design and HVAC room air motion characterization. The basic room ventilation model involves association of two enclosures, communicating laterally through an opening such as a doorway, window, corridor or incomplete dividing wall. At Colorado State University, Neymark [ 12] reports heat transfer and flow visualization data for a candidate benchmark two-room scale-model enclosure under laminar natural convection.
°iI
266
A . J . B a k e r et al. I
I
I
I
7
.... .....2...... ..+.. .......... -19.22.-..
... ....
"~. ....
-0.1
Linear FE,Beta=A -o--. Linear . . F E. , B e. t a = . 2 . QUICK-5 F V -~-QUICK-3 FV o .
",,.3
0
"~
".......
.......x",, ........
,
""t3 ...........
""'-'~
""
~ ",ix....,.
o -0.2 "~
-0.3
%-,.
C=I/2 0=I/2
-0.4
t. oo
16
8
4
3
Modal Wavelength (~---nAx) Fig. 3. Fourier stability amplification factor error for various CFD algorithms.
The small-scale water model is experimentally verified to retain laminar wall boundary layers for Ra ~< 1012, for a constant flux (hot) wall, a constant temperature cold wall, and an adiabatic partition with doorway of height h = 1-1/2 and width w -- 0.2 W (see Fig. 6a). Flow visualization experiments used dye injection, and for (flux) Ra ~ l0 l/, the hot wall wavy laminar boundary layer appeared to separate at wall mid-height, via the observed horizontal flow directed towards the top of the doorway. Upon entering the cold room, the laminar plume rises sharply upwards along the partition to the ceiling, where it turns and thickens before descending the cold wall as a wavy boundary layer. This experimental characterization, reported as symmetry mid-plane flow visualization data, is summarized in Fig. 6b. The CFD experiments simulating flow in this partitional room geometry employ Rayleigh number continuation on Ra >I 106. Figure 7a shows the final M = 52 x 42 x 26 mesh, containing approximately 57,000 nodes, determined minimal as the consequence of a meshing study starting with a uniform 24,000 node discretization. Based upon C F D solution character at Ra .~ 10 6, solution-adaptive meshing focused on solution stability (for fl = 0) in the region of the doorway and hot zone vertical buoyant plume. Momentum boundary conditions are no-slip on all partition surfaces, while tangency with drag was defined on all other boundaries (for mesh conservation) formulated via the Blasius solution for a developing, flat plate laminar boundary layer
(1.328
= t,/(Rex))
"~pU 2
2g< =
/['v'~U3,2
fl.328
p ",/t;))
= 0.008487 U 312 [dynes/cm2].
(18)
The energy equation boundary conditions were fixed cold wall temperature (Tc), and a fixed hot wall temperature (Tn) as an approximation to (assumed) constant heat flux. All other room surfaces are adiabatic. The base C F D solution at Ra = 106 readily achieved steady-state, and the computed velocity and temperature distributions on the symmetry plane are graphed in Fig. 7b and c. The predicted buoyant plume rising up the cold zone partition above the door agrees very well with the flow visualization data (Fig. 6b), but there is no indication of boundary layer separation mid-height on the hot wall. The CFD simulation continuation to Ra ~ 101° also does not produce a separation, but instead predicts a significant 3-D roll vortex that transports mass to the symmetry plane from the far-lateral reaches of the hot room. Velocity vector and temperature distributions, on select CFD cutting planes, clearly delineate the extent of the thermal boundary layers along the heated and cooled walls (Fig. 7d). These C F D data are smoothly continuous, and the predicted absence of hot wall separation is perhaps indicative of the difficulty of experimental dye visualization procedures applied to characterize complex 3-D flow fields. For such detailed 3-D flow fields, quality C F D data
Prediction of Indoor Room Air Motion : :~:.--_---.:-:--:.---:.-:=---~.~:..-:
~ ~-~-~
~ ~ - _ .: - .: - - .:.. :. -.
267
~ - ~ -.
:.
::.
a)
~f ! x~/S
t,)
,",~,~,~.,,,-,, ........ ~,,,.....,, ~.~. . . . .
x,,/S /
x~/S
I"l~
_
v
ql d
l
. . . .
ii'i
....
ii~iii~t777tTr77"77.Z~:,;
T
xJS
x,~
x,~S
d) X~/S
x~ I
lx,~
I
Fig. 4. TWS h solutions for benchmark, step-wall diffuser, Re = 390,/~ = 0, (a) velocity, (b) pressure; and for Re = 650, (c) velocity, (d) pressure, (e) pressure, ~ = 0.1.
can assist in visualization, via post-processing the tracks of Lagrangian particles for select release points. Figure 8a and b present such data, viewed perspectively from the hot and cold rooms, respectively, for a particle release near the top of the door jamb and offset 0.5 cm from the symmetry plane. Caught in the buoyant plume, the particle rises to the ceiling of the cold zone and then descends to the floor via a winding path. Once entrained by the main flow along the near-symmetry plane flow, the particle moves rapidly into the hot zone while closely following the wall boundary layers. The direction of motion is indicated along the track, and the time interval between symbols ("bubbles") is a uniform constant, hence the relative particle speed is discernible by local bubble separation distance. For a particle released slightly below doorway mid-
height, halfway between the symmetry plane and door jamb, Fig. 8c and d summarize engagement into the complex 3-D flow within the central core of the hot zone. It eventually proceeds horizontally towards the top of the door, re-enters the cold zone, and is then entrained by the dominant buoyant plume. Descending to near its release point, the particle then re-enters the hot zone. These Lagrangian particle tracks, in our view, clearly delineate the CFD-simulated (and complexly 3-D) motion of this two-room laminar natural convection flow field. Finally, to illustrate the quantitative diagnostic utility of conducting " C F D experiments," the full-scale forced convection experiment with significant buoyancy effects documented by Spitler [13] has been simulated for parametric Re'. Flow field magnitude (speed) and tem-
A . J . Baker et al.
268
Lmvel p F E D C
B A 9 8
0.0(~8 0.0448 0.0222 -0.0001 `0.0224 `0.0447 `0.0670 `0.O883
a)
L
\\\'r",...
~a..]~'~/~:..,~.~..~ ~ k - - ~ . ~ . . . ~ , ' , ~ C - ~ . _
2 `0=,
• I.~Mi p F
•E
0.0026 -0.0,,4
O
`0.0254
C B
`0.0394 `0.0554
A
`0.0674
9
`0.08,3
8
-o.o953
b)
ft x~ "E"l-.._"c-...I ~\\";'~"~'-s..l)~/l/~
,
-0.,5,3
Fig. 5. T W S h solutions for 3-D pressure distributions in the benchmark step-wall diffuser, Re = 650, (a) /~ = 0, (b)/~ = 0.1, from [3]•
Prediction of lndoor Room Air Motion
ty ~
H,.58.4
c~
269
symmetryp/ane
~...~.
Fig. 6. Partitioned-room 3-D benchmark geometry of [8]. (a) Illustration of problem definition. (b) Experimental flow visualization data for w/W = 0.2.
a)
b) i
t~
t~
"O O O M
O ¢-
Fig. 7a and b (caption overleaf).
A.J. Baker et
270
1
al.
f
Level. temp B
27
A
26
9
25
8
24 23
6
~
5
22 21 20
e-
hot wail-
c)
3 2
18
1
17
19
Level temp B A 9 8 7 6 5 4 3 2 1
27 26 25 24 23 22 21 20 19 18 17
d)
Fig. 7. CFD experiment data for two-room scale model natural convection experiment of [8]. (a) Halfdomain FE mesh, distributions on symmetry plane for Ra = 10 6 o f (b) velocity, (c) temperature, and (d) lateral cutting plane temperature and velocity distributions.
perature data were taken by Spitler for several settings in the University of Illinois full-scale ventilation test facility. The 15 air changes/hour (ACH) test employed wall heating panels, and the room ceiling and floor were insulated. The ventilation system supplied cold air to the room, through a non-vaned, half-elevation end-wall duct, with associated Archimedes number of 4.3 and Reynolds number of 15,000. The room exhaust was at floor level, kittycorner to the supply. The C F D experiments used a nominal 20,000 element meshing, with discretization extensions into both the supply and room exhaust ducts (for consistant boundary condition application). The C F D boundary conditions employed drag for velocity, computed from the 1/7th power law turbulent boundary layer model. The supply air temperature was Ts = 20°C(69.8°F), the side walls
were maintained at 27°C (83°F), and the bulk supply velocity of 1.3 ft/s was applied uniformly across the duct opening. For this large, vaneless supply opening, the 3-D room air flow field was anticipated to be low turbulence level, to be verified by conducting " C F D experiments." The (laminar) flow base diffusion level of the CFD simulation is characterized by Re-~, recall (3), which is then augmented by Re t according to
diffusionlevel ,,~~
[1 + Re'].
(19)
Hence, for Re' = 14, for example, the diffusion level in (3) is characterized by the premultiplier [1+14]/ 15,000 = 10 3. Continuing, Re' = 29 changes the diffusion level to 0.5 × 10 -2, while Re' - 149 produces a
Prediction of lndoor Room Air Motion
271
a)
c)
(t) ~lllllllllllllllllll~ IlUlllllllllllllll IIIIIHIlUlIIIIIII illlillllllnllllll IIIIIllllllllllllll " - ~. . . . InllllIP
!UIIIIIIIIIIIIHIII qlllllllllllllll 'lllllll~ii.Sil ~lllll.';ii~ I l l r o illlllllllu.~ll IIIIl~lllllk ~llim~Jllll|
"}
111:.-,I ":-n/~II/UI I l l l l l
._=~_ .
.
.
.
,i lllllllll
l
.
.
.
.
.
.
.
. .
. .
. .
i mlllUllllllllilrl, i I I . l ~ l l ~i ll l ll ~l l i l t lm
I~
. .
~
. .
. .
. .
. .
. .
. .
.
.
.
.
.
.
Iliilll~,Ulllillllll
I"l 'l l"l 'r :l l"
ll',
l "•"
Ilill:ll Illll-II
I"l l'l"l l'l " ?
.
.
.
.
.
- . . . . . .
r
IIII i ~ i
;llllll, ,~I
Illlllllllll IIIlllllllll ~.~*
........
.
.
. .
m
• . • m m mm m m m ~
•
Ilil I~ ]lllll Illll:; IIl~..',,g Illlll'llf~illllll~
!1
,
~-~:-=~ ~?" L.
i; ,
I l l l l U l l l l
- *u
IN
IIIIIIlillll II II l l l l l l l l l F Hi illllllllPVdll i • ~11111| ~ l l l l n i l ' ~liinmmmamnm[, illilJ i lilililll II I I l l l l i l l l M l l
""'",;lu-'J ,,,,u.
!lillimlilmmli
.......... - -----
Fig. 8. CFD experiment prediction of Lagrangian particle tracks in two-room natural convection experiment, Ra = 106, for release point near doorjamb : (a) view from cold room, (b) view from hot room ; and release point at door mid-height : (c) view from hot room, (d) view from cold room.
10 -2 relative diffusion level. Figure 9a presents a scale perspective of the room, showing supply and exhaust duct locations, while Fig. 9b graphs the Spitler data presentation of measured flow field speed distribution on the vertical plane bisecting (laterally) the supply duct. These "rastor graphic" data clearly confirm that the supply jet descends rapidly to the floor, and hence proceeds across it while decelerating near the exhaust plane wall. Figure 9c presents the C F D simulation results for Re' = 14, as a 3-D perspective of select velocity vectors on cutting planes, and as speed contours on the supply mid-plane (Fig. 9f). Figure 9d and e present the comparative C F D simulation velocity vector data for R e ' = 29 and R e ' = 149, respectively, while Fig. 9g and h present planar speed profiles. Clearly, the C F D simulation result agreeing best with the (available) data is for Re' = 14, BAE 29:3-C
which confirms the room flow is low turbulence level. Since p = 0 for these tests, no artificial diffusion effects were present to confuse these conclusions. Having such control is critical to conducting quality C F D experiments. SUMMARY AND CONCLUSIONS A finite element spatial semi-discretization of a Taylor weak statement forms the theoretical basis for a robust and phase-accurate C F D algorithm, for application to prediction of natural, mixed and forced convection 3D room air motion. The key formulation ingredient of stability, without excess artificial diffusion, is validated for 2-D and 3-D benchmark problems. Three-dimensional room air motion C F D simulation data, with corn-
A.J. Baker et al.
272
a)
g~ m~ mlm
b)
E~
ou rne 1
- - - 4 ft/s
d)
I
"J
,t
T
F E D G B A 9 8 ? 6 5 4 3 2 1 tJVa
143 134 124 115 105 9IS M 76 07 57 48 38 29 19 10 vlW
F
100
D C 6 A 9 e ? 6 5 4 3 2 1
87 SO 73 67 6O 53 47 4O 33 27 20 13 7
0
g)
Level vel
e)
F
85
E I) C B A 9 6 7 6 5 4 3 2 1
79 "74 68 62 S7 51 45 4O 34 26 23 17 11 e
Fig. 9. Full scale UI r o o m , Re = 15,000, Ar = 4.3, A C H = 15, (a) perspective of r o o m geometry, (b) speed raster graphic from UI experiment [13] ; select C F D experiment velocity vectors, (c) Re' = 14, (d) Re t = 29, (e) Re' = 149 ; C F D planar isovel distributions, (f) Re' = 14, (g) Re' = 29, (h) Re' = 149.
h)
Prediction o f Indoor R o o m Air M o t i o n parative experimental data, have quantized cause-effect relationships for a forced convection, empty 3-D full scale room, and a natural convection, partitioned 3-D model r o o m geometry. Acknowledgements--This work was partially supported by grant
273
MSS-9015912 from the National Science Foundation. The theory formative-phase research support provided by the American Society of Heating, Refrigerating and Air-Conditioning Engineers (ASHRAE), under project RP-464, is also gratefully acknowledged, as is the corporate sustaining member financial support of the UT CFD Laboratory, where these experiments were conducted on Unix workstations.
REFERENCES 1. P.V. Nielson, Prediction of air flow and comfort in air conditioned spaces. ASHRAE Transactions 81, II, 247 (1975). 2. Annex 20, Air flow patterns within buildings. Int. Energy Agency, A. Moser, Opr. Agent, Swiss Federal Inst. of Tech., Zurich ( 1991). 3. H. Tennekes and J. L. Lumby, A First Course in Turbulence, MIT Press, Cambridge, MA (1972). 4. D.J. Chaffin and A. J. Baker, On Taylor weak statement finite element methods for computational fluid dynamics. International Journal Numerical Methods in Fluids, in review. 5. A.J. Baker and J. W. Kim, A Taylor weak statement for hyperbolic conservation laws. International Journal Numerical Methods Fluids 7, 489-520 (1987). 6. P. T. Williams, A three dimensional, time-accurate, incompressible, Navier Stokes, finite element CFD algorithm. Ph.D. dissertation, University of Tennessee, TN (1993). 7. A.J. Baker, P. T. Williams and R. M. Kelso, Numerical calculation of room air motion, I. Math, physics & CFD modeling. ASHRAE Transactions 100, 1,514 (1994). 8. P.T. Williams and A. J. Baker, A well-posed weak statement algorithm for the pressure Poisson equations for incompressible flows. AIAA-93-3340 (1993). 9. A.J. Baker, Finite Element Computational Fluid Mechanics, Hemisphere, Washington, DC (1983). 10. A.J. Baker and D. W. Pepper, Finite Elements 1-2-3, McGraw-Hill, New York (1991). 11. B.F. Armaly, F. Durst, J. Pereira and B. Schonung, Experimental and theoretical investigation of backward-facing step flow. Journal of Fluid Mechanics 127, 473-496 (1983). 12. J. Neymark, Aperture and Prandtl number effect on interzonal natural convection. M.Sc. thesis, Colorado State Univ., CO (1988). 13. J . D . Spitler, An experimental investigation of air flow and convective heat transfer in enclosures having large ventilation flow rates. Ph.D. dissertation, University of Illinois, IL (1990).
~
Copyright © 1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0360-1323/94 $7.00 + 0.00
Pergamon 0360-1323(94)E0008-F
Development of a Robust Finite Element CFD Procedure for Predicting Indoor Room Air Motion A.J.BAKER* P. T. WILLIAMS* R. M. KELSOt The requirement is verifiable high quality computational (CFD) solutions to Reynolds-averaged Navier-Stokes equations for buoyant and turbulent flows in genuine three-dimensional room geometries. Current-practice theory~codes stabilize intrinsic dispersive error via artificial diffusion methods, with resultant compromise of genuine physical~turbulent diffusion versus convective processes. A Taylor weak statement CFD theory, and finite element code implementation, represents a potential significant accuracy/stability improvement for three-dimensional buoyant forced-, mixed- and natural-convection room airflow prediction. This theory, documentary validations, and quantitative data comparisons from conducted CFD experiments are presented.
INTRODUCTION
significant number of assumptions, it is unreasonable to expect the results of a given CFD simulation to be accurate unless (1) the problem definition is trivial, or (2) great care is exercised in design and execution of the "CFD experiment." In truth, CFD simulation is a new "experimental protocol" that under parametric variation produces a flow field data base for interpretation. Engineers have for years conducted physical experiments in the "real laboratory," where a carefully-constructed scale model is interrogated using data-measuring devices carefully selected to control error in the acquisition process. A CFD simulation experiment is analogously conducted in the " C F D laboratory," i.e. via a math model on a computer, rather than a physical model, where careful attention is (must be) paid to sources of error that can compromise the CFD data generation process. To answer the posed questions then, one must seek out the error mechanisms existing to compromise the fidelity of a CFD room air motion simulation. The distribution of mass, velocity and temperature in a given room certainly satisfies the basic conservation principles we all studied in college. Expressed in partial differential equation form, these are called the "Navier-Stokes" equations, which mathematically describe in differential calculus form the interaction among the dependent variables density (p), pressure (p), velocity (vector) field (u) and temperature (T). However, most room air flow fields are at least locally "turbulent," and even though the NavierStokes (NS) system remains valid, a "Reynolds-averaging" manipulation (cf. [3]) is required to be applied for CFD computational tractibility. Any such resolution of velocity, into mean and fluctuating (about the mean) components, involves several assumptions, and usually leads to definition of a turbulent (Reynolds) stress-strain rate constitutive hypothesis containing a turbulent "eddy viscosity." Such an expression is similar-appearing to
COMPUTATIONAL fluid dynamics (CFD) is the maturing technical discipline of development of math/physics models amenable "to computing," that are predictive regarding the motion of fluids through and around practical geometries. Indoor room air motion is the problem class of interest, and CFD is becoming applied with increasing frequency to simulation of the buoyant and perhaps-turbulent flow of air in pertinent geometries. Nielsen [1] pioneered in applying CFD to room air motion simulation, and in the ensuing two decades, progressively more indoor air applications of CFD are becoming published. These activities have focused mainly in Europe and Japan, and for example led to formation of Annex 20 [2], under the auspices of the International Energy Agency (lEA), with the goal to organize a basis for validation of CFD methods/codes for application to indoor room air motion prediction. The fundamental question to be raised regarding CFD methodology applied to room air motion prediction is, "Just how good is it, or can it be?" The companion query must be, "What issues exist to compromise the ability of a CFD method to accurately simulate a (any) given room air motion situation?" The answers to these questions are intrinsically important to the indoor air technical community. This paper addresses them, first philosophically, and then technically with development and benchmark validation of a new CFD methodology designed to alleviate limitations in the current practice. The fundamental precept to seeking the answers is that, since a//CFD procedures constitute the union of a *Department of Engineering Science and Mechanics and tSchool of Architecture, The University of Tennessee, Knoxville, TN 37996, U.S.A. 261
A . J . Baker et al.
262
Stokes' laminar viscosity constitutive law (contained in the NS system). The distinction is that in NS the genuine fluid viscosity v depends on the fluid itself (and weakly on temperature), while in Reynolds-averaged Navier Stokes (RaNS), the turbulent eddy viscosity v' depends upon the state of the flow, i.e. v' = vt(u). Therefore, a principle fundamental assumption for C F D simulation of room air motion is the form selected for the governing Reynolds-averaged Navier-Stokes PDE system. The almost-universal selection has been to append the two-equation turbulent kinetic energy (TKE) transport PDE system to RaNS, which yields a closed, mathematically well-posed system. The drawback is that the TKE model assumes the flow is.[ully turbulent, hence predicts an excessive dose of vt for the subtly-turbulent flows existing throughout indoor rooms away from HVAC supply systems and obstructions with edges. One measure of "how turbulent" a flow (simulation) is, is the turbulence Reynolds number (Re') defined as
yl V
Re' =- -- >/- O.
(1)
For a laminar flow, Re' = 0 by definition, while Re' order (102-103 ) for fully turbulent flows. Weakly turbulent flows are characterized (roughly) by 0 < Re' ~< 30 while the transition region to full turbulence occurs (essentially) within 30 ~< Re' ~< 100. The TKE fully-turbulent closure model has been adapted with "low turbulence" adjustments to account for wall layer effects, but these constructions involve many more assumptions, and remain to be validated for general application to weakly turbulent central flows away from walls, supply diffusers, etc. Philosophically then, the answer to the premier question of "how good C F D may be" for room air motion simulation rests squarely on the validity of the assumed form of the RaNS governing PDE system, including all boundary and initial conditions. Errors committed here are not assessable until more exact PDE forms are derived and validated. This arena is a "hot topic" in turbulence research, but in practice today, it is unreasonable to expect that any given C F D simulation can be "exact." Hence, one must conduct a range o f " C F D experiments" to validate C F D simulation sensitivity to parameters embedded in the RaNS turbulence closure system. An example of such a C F D experiment is presented in this paper. One must next turn to the very fundamental questions, i.e. "How good is the C F D theory/code in emulating the genuine solution to the selected RaNS system?" All C F D theories, be they finite difference (FD), finite volume (FV), finite element (FE), or any other form, can do no more than establish an approximation to the genuine RaNS solution. One fundamental flaw in any C F D theory is that it employs a mesh to support construction of the approximate solution process. Just as Reynoldsaveraging precludes obtaining the complete NS solution itself, a C F D approximation contains no resolution of solution information of two cell-wide wavelength. Figure la illustrates a I-D mesh of measure "Ax," and a C F D solution generates information only at the nodes (labelled j, j + 1, etc.). Obviously, any data of wavelength 2Ax is
a)
" j . 1~
b)
j+2 x 1;3 j~. /"-'1
Fig. 1. Sketches of spectral content resolution on a mesh of measure Ax. never seen(!) at the nodes, hence by the C F D solution process. Therefore, the C F D mesh must be sufficiently dense to admit capture of essential solution spectral content, or be totally worthless. Figure lb confirms that "5Ax" wavelength information is quite resolvable by a mesh of "measure" Ax. In practice, the C F D code execution process transmits information about the RaNS approximate solution over the nodes of the mesh. Therefore, one measure of how well this is accomplished, by competing F D / F V / F E theories, is to compute the algorithm intrinsic phase velocity error as a function of "data" of wavelength 2 = nAx, n/> 2. Using symbolic manipulation [4], this has now been unequivocally established, and Fig. 2 summarizes these data for six competing C F D theory implementations. The current practice F D and FV methodologies are theoretically confirmed to be of quite poorer performance in the short wavelength region (in this measure) in comparison to the emerging FE weak statement methods. This predictive theoretical analysis is documented further herein. One consequence of this phase error distortion, intrinsic to all C F D methods, is that the poorly transmitted information leaves behind short wavelength oscillations in the approximate solution. These "wiggles" can easily produce non-realizable data, e.g. negative temperatures, or can cascade (the RaNS PDE system is nonlinear) and cause iterative divergence of the C F D code solution process. The almost universal "CFD-fix" for this is to add numerical dissipation to the algorithm, to diffuse the wiggles. This process unequivocally introduces a numerical "artificial viscosity" coefficient va that acts in concert with the genuine fluid viscosity v and the turbulent viscosity v'. In some instances, va can be so large as to absolutely swamp the other two physics-generated diffusion mechanisms. The control of this strictly-CFD error mechanism is probably most critical to a favorable answer to the practical question, "what issues exist to compromise . . . " posed earlier. This paper highlights the theory, and presents finn validation data, regarding development of a new, mathematically robust and versatile C F D procedure for the room air motion problem class. It is based on the Taylor weak statement (TWS) theory [5], with generalizations [6,7], which has proven capable of recovering historical dissipative C F D methods, while predicting improved, more-accurate replacement methodology. This paper
Prediction of Indoor Room Air Motion 1.0
I
I
I
263
!
Central F D -,~--QUICK-3 F V o QUICK-5 FV --+--IAaearFE -n--.
0.8
• °° *" °° *'
...."""/I//
Quadratic FE -.A-.Cubic FE --~---
...., //,,t... /'; /
•"
i "S
.."
0.6
///
,,"
/
o.°
ra ! ; / i
..
/
i '.
..'
/
/ /
;
i
.°
,/
.-"
C=1/2 0=1/2
!
."
/
..."
8
"
ii
!!
,,
/
,"/"
i ;
#
• o°
0.4
,o
• °° °°°
0.2 °, -° °° °°
°° o°
°°
/,s
°°
J
/s ,° °L'JI
s A"
i,, S
/
°°"
,
o°°° /
s s.s" ..s"
/
/ • i° ~,'
s.J
I"
•
t¢ °
I'
..[3"
°° °° °°
Olb oo
16
8
4
3
Modal Wavelength (X----~x) Fig. 2. Fourier stability phase error comparisons for various CFD algorithms.
highlights the TWS theory, applied to RaNS, and its implementation via a finite element semi-discretization of a Galerkin weak statement. Thereafter, the implicit temporal semi-discretization admits time-accurate, unsteady predictions, and sparse matrix methods are employed to iteratively solve the resulting non-linear algebraic equation system. Theoretical and computational performance validations are presented for select 2-D and 3-D benchmark problems, and for two genuine 3-D room air motion experiments.
THE PROBLEM S T A T E M E N T
A statistical manipulation, c.f. [3], produces a RaNS PDE system governing incompressible buoyant room air flows. Following nondimensionalization, and recalling the turbulent Reynolds number (Re') of the turbulence closure model, the 3-D PDE system governing mass, momentum and energy transport processes, hence requiring solution, is
80
~ (
×
= 0
(2)
au,+~( dxj u~uj+ P'~°
£e(u,) = - ~
l (l+Re')t-~xj+-~xJ)-ArOS,g=O
Re BAE 29:3-B
(3)
+ Pr'Jdxjj-s°
= 0.
(4)
The dependent variable set in (2)-(4) contains the meanflow velocity vector u ~ ui(xj,t), the kinematic pressure P=P/Po and the potential temperature ® - ~ ( T Tmin)/(Trnax-Tmin). The Boussinesq buoyancy body-force approximation is made in (3), where 6ig is the gravity unit vector. The dimensionless groups appearing in (3) and (4) exert a controlling impact on solution character. The definitions for Reynolds, Grashoff, Prandtl and Archimedes numbers are
UL
Re=--,
Gr=
3)
g fl(tma
x - - Trnin ) Z 3 V2
P r - poCpV Ar = GrRe -2. =
~e(po) = ~
1 Re
L~(O) = --~ + ~jxj ujO
k
'
(5)
'
The familiar Reynolds number (Re) is based on userselected reference velocity (U) and length (L) scales and the fluid kinematic viscosity v. The Grashoff number (Gr) contains a reference temperature spread, the gravity constant and the fluid compressibility factor (fl), roughly equal to inverse absolute temperature. The Prandtl number (Pr) is based on the reference (constant) density (P0), specific heat (cp) and thermal conductivity (k). The
A.J. Baker et al.
264
remaining parameters in (3) and (4) are turbulent Reynolds number Re', defined previously in (1), and the turbulent Prandtl number Pf, which nominally is defined identical to Pr for air. The character of solutions to (2)-(4) depends on these non-D groups, the distribution of Re' and boundary conditions. Since Re > 0 and Re' ~> 0, isolated equations (3) and (4) are elliptic boundary-value, hence the corresponding members of the mean-flow state variable q (x, t) = {u~, u2, u3, O} r, and/or their normal derivatives, are required to be known all around the problem statement boundary. The solution process for these equations is thus classic initial-boundary value, up to enforcement of the continuity constraint that us must be divergencefree, i.e. the velocity field must also satisfy eq. (2). This differential constraint is only enforceable approximately, via one of several inexact CFD theories, as a modelled action for pressure P. A thorough theoretical study [6] verifies a robust procedure appends a Poisson variable ~b satisfying the boundary value problem
where superscript tilde denotes any approximation to the true (i.e. divergence-free) vector velocity u, and n is the time-step index. The parametric variable q~ is then manipulated [6] into a kinematic model of the action of pressure via P"+ ~ _~ P " - d?.+ l/OAt, where subscript "A" denotes approximation and P" is the genuine pressure (from the previous solution at time t,). The CFD approximation PA must then be replaced, upon time-step convergence, with the genuine pressure P,+ ~, which is the solution to
t?2P - -~2 (u~u/-5fl(P) = (~x~ + Ox,axj ~ ×
+
coU(q) = ;~ + ~_ (.t - f ' i ) - s = 0
(8)
~(q,) = V2q , - s(q) = 0.
(9)
In (8),Ji = J)(q) is the kinetic" flux vector, containing convection and pressure effects, while fj' is the dissipative flux vector embodiment of the turbulence closure model (Re') plus natural diffusion (Re, Pr). In (9), the Laplacian operates on q~ = {q6, p}r, and s = s(q) is the respective source dependence on the solution to (8). Any CFD algorithm generates (only) an approximation to the true solution to (8) and (9), via a denumerable set of decisions leading to an algebraic system amenable "to computing." As discussed, a dominating dispersive error mechanism is intrinsic to this process for (8), ultimately leading to the near-universal use of "artificial diffusion" mechanisms for stabilizing computations, e.g. " Q U I C K " differencing. The emerging encompassing Taylor weak statement theory, [5,6], extended to systems (8), yields creation of the PDE system replacement for (8) in the form ~caC(q) - L,a(q) _ flAt c?~x--~( A j A k ~ ) = 0 ,
(10)
where Aj ----0fj/0 q is the RaNS kinetic flux vector Jacobian (matrix). The replacement of (8) by (10) in the CFD implementation process contributes directly to a tensor invariant, highly phase selective numerical dissipation mechanism, as will be verified. Hence, upon replacement of (8) by (10), any approximation qN(x, t) to the true solution q(x, t) to (8) and (9) in the continuum is of the functional form N
1
q(x, t) ~ q"¢ (x, t) = ~ LFj(X)Qj (t), j=J
~e(1 +Re')
)
-- A r ~ xs• 6jo. = 0.
(7)
The third (spatial) derivatives on velocity in (7) are computationally intractable to form to an acceptable discrete order of accuracy. The theoretical resolution is now achieved via a Galerkin weak statement theory [8], which lowers required differentiability order to first. Since (6) and (7) are elliptic, surrounding boundary conditions are required for well-posedness. For (6), the sole known boundary conditions are homogeneous Neumann (normal derivative), while non-homogeneous Neumann conditions are derived for (7) using (3), c.f. [8]. Finally, a closure model for turbulent Reynolds number Re'(x, t) is required, which is constrained herein to a scalar factor, as appropriate for CFD experiment comparative purposes.
(11)
where tPj(x), 1 ~
(12)
be satisfied for any (all!) given known functions ~(x). If the theory-designer picks ~i(x) and ~i(x) identical, then a "Galerkin" (finite element, FE) weak statement results. Conversely, if one chooses q)~(x) to be constants, then a "finite volume" (FV) weak statement form is produced. To actually evaluate the integrals defined in (12), for genuine problems, a spatial semi-discretization facilitates the calculus operations, whereupon, for either an FE or FV selection, (12) always produces a matrix ordinary differential equation (ODE) system of the form [9]
WEAK STATEMENT CFD ALGORITHM
d{Q} TWS h = [ M ] ~ + {R(Q)} = {0}.
A CFD algorithm for (2)-(7) generates solutions parameterized by Re, Gr, Pr and Re'. The CFD theoryderived dependent variable set (the "state variable") is q(x, t) = {ul, uz, u~, ®; 4), P} T, which at the semi-colon partition satisfies the PDE systems
In (13), the superscript h notation emphasizes that a "discretization" ~h of the domain I'~ of (8)-(10) has been invoked, in achieving this matrix expression. Further therein, [M] and {R(Q)} are global (nodal) rank square
(13)
Prediction of Indoor Room Air Motion and column matrices, respectively, and {Q} = {Q(t)} is the (unknown) state variable approximation coefficient (set) Qj at the nodes of the mesh defined by fP. Any ODE algorithm uses (13) directly for time derivative evaluation. Selecting the one-step, 0 implicit Euler family for example, and substituting (13) as needed [91, produces the terminal matrix equation system "for computing" as
{FQ} = [M] {Q.+I - Q.} + At(O{R(Q)}.+I + ( 1 - 0 ) {R(Q)}.) = {0},
(14)
where subscript "n" remains time level notation, The discrete Galerkin weak statement (12) for (9) directly produces the companion matrix systems
{FQ~} --- [D] {QA}- {S(Q) } = {0},
(15)
where square matrix [D] is the discrete Laplacian, {Q~} is the nodal array of 4? or eh, and {S(Q)} contains dependence on qh via the solution {Q(t)} from (14). The generic solution process for (I 4) and (15) involves use of an appropriate linear algebra approximation to the Newton iteration algorithm
[JA C] .~60~P+ , .+ l = - {FQ}P+~
(16)
P
p+l {Q}.+I ,+i = {Q}.+ 2 {6Q}.+~,
(17)
265
the point of destablizing the quasi-Newton iterative process for (16) and (17). Hence, the central stability issue with a CFD algorithm is control of this dispersive error mechanism, such that a convergent iteration process can produce a "smooth" solution, free of "wiggles." To date, there are no viable methods to annihilate this mode, hence the remaining (universally-applied) "fix" is to dissipate the oscillations via an added numerical diffusion process. For a linear model of equation (8), the Fourier stability analysis employed to develop Fig. 1 is extended to computation of the amplification factor for a variety of CFD algorithms. The exact solution amplification factor is unity, i.e. no (artificial) damping is present, hence departures from unity, as a function of modal wavelength ~ = nAx, characterize algorithm artificial damping. Figure 3, from [4], compares the error in (unit) amplification factor for a range of FD, FV and linear basis FE CFD algorithms, the latter for various TWS-theory E-levels, which are user specifiable. The horizontal line at zero error confirms that the non-artifically diffused FD and FE algorithms (fl --- 0) have no intrinsic dissipation error. The QUICK3 FV algorithm is most dissipative, especially in the short wavelength region, while the current FE algorithm with (user-definable) modest fl yields a more phase-selective dissipation mechanism.
i=O
where [JAC] is the Jacobian of (14) and (15). The Newton Jacobian defined in (16) is readily formed using calculus operations, cf. [10]. Thereafter, the range of quasi-Newton iterative procedures use approximations thereto. For this paper, a variable-segregated, 3-D quasi-Newton Jacobian construction was iteratively cycled, using a diagonally-preconditioned GMRES algorithm [6], in concert with an incomplete-Cholesky preconditioned conjugate gradient (PCG) solver for (15). ACCURACY, STABILITY AND DISPERSION ERROR The weak statement theory (8)-(12) confirms that
required fundamental CFD algorithm decisions focus on selection of the trial space Wi (x), and the test space ¢i(x), for defining qU and the TWS N error constraint respectively. Under a spatial semi-discretization, these decisions become commensurate with the FD concept of algorithm order of accuracy, hence asymptotic convergence under mesh refinement. Linear basis FE and traditional FV algorithms are typically second-order accurate, and coarse meshes aggravate the intrinsic semidiscrete approximation phase error mechanisms; recall Fig. 1. In practice, algorithm stability is of much greater practical significance than convergence for genuine RaNS problems, since the available mesh density is always limited by computer resources. The practical issue is thereby control of the third-order dispersive error mechanism intrinsic to all discrete approximation methodologies. As discussed, any mesh of measure Ax in one-dimension is incapable of resolving solution information of wavelength 2Ax. The non-linearity of a RaNS equation system amplifies this 2Ax oscillatory dispersive error mode, often to
RESULTS AND DISCUSSION The benchmark problem of specific pertinence to prediction of room air motion is the close-coupled step wall diffuser, cf. [111, which models the essence of a supply outlet. Figure 4 summarizes the 2-D Galerkin linear basis FE WSh steady state laminar solution obtained at Re = 390 in terms of (a) velocity vector uh distribution, with experimental location of separation zone reattachment (xl/S) noted, and (b) pressure distribution. In proceeding to Re = 650, a top-wall secondary separation occurs (Fig. 4c), in qualitative agreement with the (3-D) experiment; however, the discretized Poisson equation pressure field is noticeably polluted by a dispersive error mode (Fig. 4d). Homogeneous Poisson (Laplace) solutions are always smooth, for regular data, hence these oscillations must be due to the non-homogeniety in (7) involving derivatives of the CFD discrete approximation uh. Switching to the TWS h dissipative algorithm, using = 0.1 does not measurably alter u h, but this "filtering" totally annihilates the pressure solution oscillations (Fig. 4e). The identical trend is recorded in the fully 3-D benchmark simulation, Fig. 5. Additional details are given in [61. Interest in natural convection in geometrically-complex enclosures, with restricted communication, has been stimulated by applications involving energy-efficient passive-solar buildings, cold-storage design and HVAC room air motion characterization. The basic room ventilation model involves association of two enclosures, communicating laterally through an opening such as a doorway, window, corridor or incomplete dividing wall. At Colorado State University, Neymark [ 12] reports heat transfer and flow visualization data for a candidate benchmark two-room scale-model enclosure under laminar natural convection.
°iI
266
A . J . B a k e r et al. I
I
I
I
7
.... .....2...... ..+.. .......... -19.22.-..
... ....
"~. ....
-0.1
Linear FE,Beta=A -o--. Linear . . F E. , B e. t a = . 2 . QUICK-5 F V -~-QUICK-3 FV o .
",,.3
0
"~
".......
.......x",, ........
,
""t3 ...........
""'-'~
""
~ ",ix....,.
o -0.2 "~
-0.3
%-,.
C=I/2 0=I/2
-0.4
t. oo
16
8
4
3
Modal Wavelength (~---nAx) Fig. 3. Fourier stability amplification factor error for various CFD algorithms.
The small-scale water model is experimentally verified to retain laminar wall boundary layers for Ra ~< 1012, for a constant flux (hot) wall, a constant temperature cold wall, and an adiabatic partition with doorway of height h = 1-1/2 and width w -- 0.2 W (see Fig. 6a). Flow visualization experiments used dye injection, and for (flux) Ra ~ l0 l/, the hot wall wavy laminar boundary layer appeared to separate at wall mid-height, via the observed horizontal flow directed towards the top of the doorway. Upon entering the cold room, the laminar plume rises sharply upwards along the partition to the ceiling, where it turns and thickens before descending the cold wall as a wavy boundary layer. This experimental characterization, reported as symmetry mid-plane flow visualization data, is summarized in Fig. 6b. The CFD experiments simulating flow in this partitional room geometry employ Rayleigh number continuation on Ra >I 106. Figure 7a shows the final M = 52 x 42 x 26 mesh, containing approximately 57,000 nodes, determined minimal as the consequence of a meshing study starting with a uniform 24,000 node discretization. Based upon C F D solution character at Ra .~ 10 6, solution-adaptive meshing focused on solution stability (for fl = 0) in the region of the doorway and hot zone vertical buoyant plume. Momentum boundary conditions are no-slip on all partition surfaces, while tangency with drag was defined on all other boundaries (for mesh conservation) formulated via the Blasius solution for a developing, flat plate laminar boundary layer
(1.328
= t,/(Rex))
"~pU 2
2g< =
/['v'~U3,2
fl.328
p ",/t;))
= 0.008487 U 312 [dynes/cm2].
(18)
The energy equation boundary conditions were fixed cold wall temperature (Tc), and a fixed hot wall temperature (Tn) as an approximation to (assumed) constant heat flux. All other room surfaces are adiabatic. The base C F D solution at Ra = 106 readily achieved steady-state, and the computed velocity and temperature distributions on the symmetry plane are graphed in Fig. 7b and c. The predicted buoyant plume rising up the cold zone partition above the door agrees very well with the flow visualization data (Fig. 6b), but there is no indication of boundary layer separation mid-height on the hot wall. The CFD simulation continuation to Ra ~ 101° also does not produce a separation, but instead predicts a significant 3-D roll vortex that transports mass to the symmetry plane from the far-lateral reaches of the hot room. Velocity vector and temperature distributions, on select CFD cutting planes, clearly delineate the extent of the thermal boundary layers along the heated and cooled walls (Fig. 7d). These C F D data are smoothly continuous, and the predicted absence of hot wall separation is perhaps indicative of the difficulty of experimental dye visualization procedures applied to characterize complex 3-D flow fields. For such detailed 3-D flow fields, quality C F D data
Prediction of Indoor Room Air Motion : :~:.--_---.:-:--:.---:.-:=---~.~:..-:
~ ~-~-~
~ ~ - _ .: - .: - - .:.. :. -.
267
~ - ~ -.
:.
::.
a)
~f ! x~/S
t,)
,",~,~,~.,,,-,, ........ ~,,,.....,, ~.~. . . . .
x,,/S /
x~/S
I"l~
_
v
ql d
l
. . . .
ii'i
....
ii~iii~t777tTr77"77.Z~:,;
T
xJS
x,~
x,~S
d) X~/S
x~ I
lx,~
I
Fig. 4. TWS h solutions for benchmark, step-wall diffuser, Re = 390,/~ = 0, (a) velocity, (b) pressure; and for Re = 650, (c) velocity, (d) pressure, (e) pressure, ~ = 0.1.
can assist in visualization, via post-processing the tracks of Lagrangian particles for select release points. Figure 8a and b present such data, viewed perspectively from the hot and cold rooms, respectively, for a particle release near the top of the door jamb and offset 0.5 cm from the symmetry plane. Caught in the buoyant plume, the particle rises to the ceiling of the cold zone and then descends to the floor via a winding path. Once entrained by the main flow along the near-symmetry plane flow, the particle moves rapidly into the hot zone while closely following the wall boundary layers. The direction of motion is indicated along the track, and the time interval between symbols ("bubbles") is a uniform constant, hence the relative particle speed is discernible by local bubble separation distance. For a particle released slightly below doorway mid-
height, halfway between the symmetry plane and door jamb, Fig. 8c and d summarize engagement into the complex 3-D flow within the central core of the hot zone. It eventually proceeds horizontally towards the top of the door, re-enters the cold zone, and is then entrained by the dominant buoyant plume. Descending to near its release point, the particle then re-enters the hot zone. These Lagrangian particle tracks, in our view, clearly delineate the CFD-simulated (and complexly 3-D) motion of this two-room laminar natural convection flow field. Finally, to illustrate the quantitative diagnostic utility of conducting " C F D experiments," the full-scale forced convection experiment with significant buoyancy effects documented by Spitler [13] has been simulated for parametric Re'. Flow field magnitude (speed) and tem-
A . J . Baker et al.
268
Lmvel p F E D C
B A 9 8
0.0(~8 0.0448 0.0222 -0.0001 `0.0224 `0.0447 `0.0670 `0.O883
a)
L
\\\'r",...
~a..]~'~/~:..,~.~..~ ~ k - - ~ . ~ . . . ~ , ' , ~ C - ~ . _
2 `0=,
• I.~Mi p F
•E
0.0026 -0.0,,4
O
`0.0254
C B
`0.0394 `0.0554
A
`0.0674
9
`0.08,3
8
-o.o953
b)
ft x~ "E"l-.._"c-...I ~\\";'~"~'-s..l)~/l/~
,
-0.,5,3
Fig. 5. T W S h solutions for 3-D pressure distributions in the benchmark step-wall diffuser, Re = 650, (a) /~ = 0, (b)/~ = 0.1, from [3]•
Prediction of lndoor Room Air Motion
ty ~
H,.58.4
c~
269
symmetryp/ane
~...~.
Fig. 6. Partitioned-room 3-D benchmark geometry of [8]. (a) Illustration of problem definition. (b) Experimental flow visualization data for w/W = 0.2.
a)
b) i
t~
t~
"O O O M
O ¢-
Fig. 7a and b (caption overleaf).
A.J. Baker et
270
1
al.
f
Level. temp B
27
A
26
9
25
8
24 23
6
~
5
22 21 20
e-
hot wail-
c)
3 2
18
1
17
19
Level temp B A 9 8 7 6 5 4 3 2 1
27 26 25 24 23 22 21 20 19 18 17
d)
Fig. 7. CFD experiment data for two-room scale model natural convection experiment of [8]. (a) Halfdomain FE mesh, distributions on symmetry plane for Ra = 10 6 o f (b) velocity, (c) temperature, and (d) lateral cutting plane temperature and velocity distributions.
perature data were taken by Spitler for several settings in the University of Illinois full-scale ventilation test facility. The 15 air changes/hour (ACH) test employed wall heating panels, and the room ceiling and floor were insulated. The ventilation system supplied cold air to the room, through a non-vaned, half-elevation end-wall duct, with associated Archimedes number of 4.3 and Reynolds number of 15,000. The room exhaust was at floor level, kittycorner to the supply. The C F D experiments used a nominal 20,000 element meshing, with discretization extensions into both the supply and room exhaust ducts (for consistant boundary condition application). The C F D boundary conditions employed drag for velocity, computed from the 1/7th power law turbulent boundary layer model. The supply air temperature was Ts = 20°C(69.8°F), the side walls
were maintained at 27°C (83°F), and the bulk supply velocity of 1.3 ft/s was applied uniformly across the duct opening. For this large, vaneless supply opening, the 3-D room air flow field was anticipated to be low turbulence level, to be verified by conducting " C F D experiments." The (laminar) flow base diffusion level of the CFD simulation is characterized by Re-~, recall (3), which is then augmented by Re t according to
diffusionlevel ,,~~
[1 + Re'].
(19)
Hence, for Re' = 14, for example, the diffusion level in (3) is characterized by the premultiplier [1+14]/ 15,000 = 10 3. Continuing, Re' = 29 changes the diffusion level to 0.5 × 10 -2, while Re' - 149 produces a
Prediction of lndoor Room Air Motion
271
a)
c)
(t) ~lllllllllllllllllll~ IlUlllllllllllllll IIIIIHIlUlIIIIIII illlillllllnllllll IIIIIllllllllllllll " - ~. . . . InllllIP
!UIIIIIIIIIIIIHIII qlllllllllllllll 'lllllll~ii.Sil ~lllll.';ii~ I l l r o illlllllllu.~ll IIIIl~lllllk ~llim~Jllll|
"}
111:.-,I ":-n/~II/UI I l l l l l
._=~_ .
.
.
.
,i lllllllll
l
.
.
.
.
.
.
.
. .
. .
. .
i mlllUllllllllilrl, i I I . l ~ l l ~i ll l ll ~l l i l t lm
I~
. .
~
. .
. .
. .
. .
. .
. .
.
.
.
.
.
.
Iliilll~,Ulllillllll
I"l 'l l"l 'r :l l"
ll',
l "•"
Ilill:ll Illll-II
I"l l'l"l l'l " ?
.
.
.
.
.
- . . . . . .
r
IIII i ~ i
;llllll, ,~I
Illlllllllll IIIlllllllll ~.~*
........
.
.
. .
m
• . • m m mm m m m ~
•
Ilil I~ ]lllll Illll:; IIl~..',,g Illlll'llf~illllll~
!1
,
~-~:-=~ ~?" L.
i; ,
I l l l l U l l l l
- *u
IN
IIIIIIlillll II II l l l l l l l l l F Hi illllllllPVdll i • ~11111| ~ l l l l n i l ' ~liinmmmamnm[, illilJ i lilililll II I I l l l l i l l l M l l
""'",;lu-'J ,,,,u.
!lillimlilmmli
.......... - -----
Fig. 8. CFD experiment prediction of Lagrangian particle tracks in two-room natural convection experiment, Ra = 106, for release point near doorjamb : (a) view from cold room, (b) view from hot room ; and release point at door mid-height : (c) view from hot room, (d) view from cold room.
10 -2 relative diffusion level. Figure 9a presents a scale perspective of the room, showing supply and exhaust duct locations, while Fig. 9b graphs the Spitler data presentation of measured flow field speed distribution on the vertical plane bisecting (laterally) the supply duct. These "rastor graphic" data clearly confirm that the supply jet descends rapidly to the floor, and hence proceeds across it while decelerating near the exhaust plane wall. Figure 9c presents the C F D simulation results for Re' = 14, as a 3-D perspective of select velocity vectors on cutting planes, and as speed contours on the supply mid-plane (Fig. 9f). Figure 9d and e present the comparative C F D simulation velocity vector data for R e ' = 29 and R e ' = 149, respectively, while Fig. 9g and h present planar speed profiles. Clearly, the C F D simulation result agreeing best with the (available) data is for Re' = 14, BAE 29:3-C
which confirms the room flow is low turbulence level. Since p = 0 for these tests, no artificial diffusion effects were present to confuse these conclusions. Having such control is critical to conducting quality C F D experiments. SUMMARY AND CONCLUSIONS A finite element spatial semi-discretization of a Taylor weak statement forms the theoretical basis for a robust and phase-accurate C F D algorithm, for application to prediction of natural, mixed and forced convection 3D room air motion. The key formulation ingredient of stability, without excess artificial diffusion, is validated for 2-D and 3-D benchmark problems. Three-dimensional room air motion C F D simulation data, with corn-
A.J. Baker et al.
272
a)
g~ m~ mlm
b)
E~
ou rne 1
- - - 4 ft/s
d)
I
"J
,t
T
F E D G B A 9 8 ? 6 5 4 3 2 1 tJVa
143 134 124 115 105 9IS M 76 07 57 48 38 29 19 10 vlW
F
100
D C 6 A 9 e ? 6 5 4 3 2 1
87 SO 73 67 6O 53 47 4O 33 27 20 13 7
0
g)
Level vel
e)
F
85
E I) C B A 9 6 7 6 5 4 3 2 1
79 "74 68 62 S7 51 45 4O 34 26 23 17 11 e
Fig. 9. Full scale UI r o o m , Re = 15,000, Ar = 4.3, A C H = 15, (a) perspective of r o o m geometry, (b) speed raster graphic from UI experiment [13] ; select C F D experiment velocity vectors, (c) Re' = 14, (d) Re t = 29, (e) Re' = 149 ; C F D planar isovel distributions, (f) Re' = 14, (g) Re' = 29, (h) Re' = 149.
h)
Prediction o f Indoor R o o m Air M o t i o n parative experimental data, have quantized cause-effect relationships for a forced convection, empty 3-D full scale room, and a natural convection, partitioned 3-D model r o o m geometry. Acknowledgements--This work was partially supported by grant
273
MSS-9015912 from the National Science Foundation. The theory formative-phase research support provided by the American Society of Heating, Refrigerating and Air-Conditioning Engineers (ASHRAE), under project RP-464, is also gratefully acknowledged, as is the corporate sustaining member financial support of the UT CFD Laboratory, where these experiments were conducted on Unix workstations.
REFERENCES 1. P.V. Nielson, Prediction of air flow and comfort in air conditioned spaces. ASHRAE Transactions 81, II, 247 (1975). 2. Annex 20, Air flow patterns within buildings. Int. Energy Agency, A. Moser, Opr. Agent, Swiss Federal Inst. of Tech., Zurich ( 1991). 3. H. Tennekes and J. L. Lumby, A First Course in Turbulence, MIT Press, Cambridge, MA (1972). 4. D.J. Chaffin and A. J. Baker, On Taylor weak statement finite element methods for computational fluid dynamics. International Journal Numerical Methods in Fluids, in review. 5. A.J. Baker and J. W. Kim, A Taylor weak statement for hyperbolic conservation laws. International Journal Numerical Methods Fluids 7, 489-520 (1987). 6. P. T. Williams, A three dimensional, time-accurate, incompressible, Navier Stokes, finite element CFD algorithm. Ph.D. dissertation, University of Tennessee, TN (1993). 7. A.J. Baker, P. T. Williams and R. M. Kelso, Numerical calculation of room air motion, I. Math, physics & CFD modeling. ASHRAE Transactions 100, 1,514 (1994). 8. P.T. Williams and A. J. Baker, A well-posed weak statement algorithm for the pressure Poisson equations for incompressible flows. AIAA-93-3340 (1993). 9. A.J. Baker, Finite Element Computational Fluid Mechanics, Hemisphere, Washington, DC (1983). 10. A.J. Baker and D. W. Pepper, Finite Elements 1-2-3, McGraw-Hill, New York (1991). 11. B.F. Armaly, F. Durst, J. Pereira and B. Schonung, Experimental and theoretical investigation of backward-facing step flow. Journal of Fluid Mechanics 127, 473-496 (1983). 12. J. Neymark, Aperture and Prandtl number effect on interzonal natural convection. M.Sc. thesis, Colorado State Univ., CO (1988). 13. J . D . Spitler, An experimental investigation of air flow and convective heat transfer in enclosures having large ventilation flow rates. Ph.D. dissertation, University of Illinois, IL (1990).