v'/EPA
United States
Environmental Protection
Agency
Environmental Sciences Research EPA-600/4-78-050
Laboratory August 1978
Research Triangle Park NC 27711
JW . ..' .
Research and Development
Turbulence Modeling
Applied to Buoyant
Plumes
PROPERTY OF
DIVISION
OF
METEOROLOGY
-------
RESEARCH REPORTING SERIES
Research reports of the Office of Research and Development, U S Environmental
Protection Agency, have been grouped into nine series. These nine broad cate-
gories were established to facilitate further development and application of en-
vironmental technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in related fields.
The nine series are'
1 Environmental Health Effects Research
2 Environmental Protection Technology
3. Ecological Research
4. Environmental Monitoring
5. Socioeconomic Environmental Studies
6. Scientific and Technical Assessment Reports (STAR)
7 Interagency Energy-Environment Research and Development
8. "Special" Reports
9 Miscellaneous Reports
This report has been assigned to the ENVIRONMENTAL MONITORING series.
This series describes research conducted to develop new or improved methods
and instrumentation for the identification and quantification of environmental
pollutants at the lowest conceivably significant concentrations. It also includes
studies to determine the ambient concentrations of pollutants in the environment
and/or the variance of pollutants as a function of time or meteorological factors.
This document is available to the public through the National Technical Informa-
tion Service, Springfield, Virginia 22161.
-------
EPA-600/4-78-050
August 1978
TURBULENCE MODELING APPLIED TO BUOYANT PLUMES
by
M. E. Teske, W. S. Lewellen, and H. S. Segur
Aeronautical Research Associates of Princeton, Inc
Princeton, New Jersey 085^0
Contract No. 68-02-2285
Project Officer
Francis S. Binkowski
Meteorology and Assessment Division
Environmental Sciences Research Laboratory
Research Triangle Park, North Carolina 27711
ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NORTH CAROLINA 27711
-------
DISCLAIMER
This report has been reviewed by the Environmental Sciences
Research Laboratory, U.S. Environmental Protection Agency, and
approved for publication. Approval does not signify that the
contents necessarily reflect the views and policies of the U.S.
Environmental Protection Agency, nor does mention of trade names
or commercial products constitute endorsement or recommendation
for use.
ii
-------
ABSTRACT
This report details the progress made at A.R.A.P. between
March 1976 and April 1977 towards the goal of developing a viable
computer model based on second-order closure of the turbulent cor-
relation equations for predicting the fate of nonchemically re-
acting contaminants released in the atmospheric boundary layer.
The invariant turbulence model discussed in previous reports has
been extended to compute the development of buoyant plumes. Nu-
merical program capability has been extended by improving the
speed and accuracy of the two-dimensional unsteady turbulent flow
calculation. Plume calculations are made for buoyant plumes
rising into a stable quiescent atmosphere and stable, neutral,
and unstable moving atmospheres. An examination of the application
of an integral approach to our turbulent boundary layer model is
also included.
This report was submitted in partial fulfillment of Contract
No. EPA-68-02-2285 by Aeronautical Research Associates of
Princeton, Inc. under the sponsorship of the Environmental Protec-
tion Agency. This report covers a period from March 1976 to
April 1977, and work was completed as of April, 1977.
iii
-------
CONTENTS
Abstract ill
Figures vi
Abbreviations and Symbols ix
1. Introduction 1
2. Conclusions 3
3. Recommendations 4
4. Computer Code Improvements and Modifications. . . 5
5. Model Comparison With Two-Dimensional Jets and
Plumes Into a Quiescent Neutral Fluid 9
6. Two-Dimensional Plume In a Moving Stream 22
7. Buoyant Plume Rise Into a Stable Quiescent
Atmosphere 26
8. Plume Rise into Moving Atmospheres 32
9. An Examination of an Integral Approach to Use in
Regional Air Quality Models 37
References 44
v
-------
FIGURES
Number . Page
1 Comparison of the predicted cross-sectional
mean velocity profile in a two-dimensional
self-similar jet with the experimental
observations of Kotsovinos (19) 13
2 Comparison of the predicted cross-sectional
mean temperature profile in a two-dimensional
self-similar jet with the observations of
Kotsovinos (19) 13
3 Comparison of the predicted cross-sectional
profile of the streaming velocity variance
in a two-dimensional self-similar jet. The
observations are from Kotsovinos (19) and
Davies, et al. (20) 14
Comparison of the predicted cross-sectional
profile of the temperature variance in a
two-dimensional self-similar jet. The
observations are from Kotsovinos (19), Davies,
et al. (20), and Uberoi and Singh (21)
Comparison of the predicted cross-sectional
mean velocity profile in a two-dimensional
self-similar plume with the experimental
observations of Kotsovinos (19) .......... 15
Comparison of the predicted cross-sectional
mean temperature profile in a two-dimensional
self-similar plume with the observations of
Kotsovinos (19) .................. 15
Comparison of the predicted cross-sectional
profile of the streaming velocity variance
in a two-dimensional self-similar plume
with the observations of Kotsovinos (19) ..... 17
Comparison of the predicted cross-sectional
profile of the temperature variance in a two-
dimensional self-similar plume with the
observations of Kotsovinos (19) .......... 17
vi
-------
Number
9 Decay of the normalized buoyancy flux as a
function of downstream distance for a jet
transition to a plume (data from Kotsovinos
and List (22)) 18
10 Growth of the local Richardson number as
a function of downstream distance for a
jet transition to a plume (data from
Kotsovinos and List (22)) 18
11 Predicted variation of the entrainment
coefficient in two-dimensional jets and
plumes compared with the average observed
values of Kotsovinos (19) 21
12 The experimental result of run 83 of Ceder-
wall (24) for a vertical buoyant slot jet
exiting into a cross stream. The box size
is our simulation region 24
13 An intensity plot of the variation of jet
temperature within the box defined in
Figure 12 for the Cederwall experiment.
Maximum intensity corresponded to maximum
temperature 25
14 Prediction of the development of a rising
plume into a quiescent stable atmosphere.
Shown here are the contour lines at several
times after flow initialization for the
q2 = 0.7 m2/s2 value of turbulent
kinetic energy 29
15 Prediction of the development of a rising
plume into a quiescent stable atmosphere.
Shown here are various stream-function
contours at several times after flow
initialization 30
16 Prediction of the downstream development
of a moving buoyant plume into a stable
atmosphere. Shown here are various pas-
sive tracer (smoke) contours at several
distances downwind of the initialization
plane 35
vii
-------
Number
17 Buoyant plume centerllne rise height as a
function of distance downstream of release
point. Curves show the effects of unstable,
neutral, and stable background temperature
gradients. A0 is the maximum difference
iricix
between the initial plume temperature and
the local ambient atmospheric temperature 36
18 Surface shear velocity as a function of
Rossby number in a steady neutral planetary
boundary layer. The exact solution is
found from a solution of the partial dif-
ferential equations of motion; the two
dashed curves represent the results of our
simplified integral theory. The decay
parameter (1 or 0.2) multiplies the assumed
exponential behavior of q in Equation 54 ^3
19 Surface windshear angle as a function of
Rossby number for the same conditions as
Figure 18 43
viii
-------
LIST OP ABBREVIATIONS AND SYMBOLS
A,bjS,v ,a invariant model constants (0 .75, 0.125,1.8 ,0. 3,0. 65',
\*/
a., ,a? ,a~ ,aj, ,X , B integral model constants
c phase speed
C mean smoke (passive tracer) concentration level
D diameter, or jet width
f Coriolis parameter (2n sin <|>)
P buoyancy flux
Pr Proude number (Equation 31)
g gravity
k wave number
a distance between cars
c
m momentum flux (Equation 18)
N Brunt-Vaisala frequency (Equation 2)
q square root of twice the turbulent kinetic energy
Q driving heat flux
r radius
Ri Richardson number (Equation 19)
Ro Rossby number (U /fz )
t time
u* surface shear stress velocity
U geostrophic velocity
o
U. mean fluid velocity
U freestream velocity
u.u. velocity correlation ensemble average
u.0 heat flux correlation ensemble average
x. coordinate (x,y,z)
z hydrodynamic roughness
z Ł distance to where U falls to one-half its
centerline value
IX
-------
n
0
eP
A
V
P
fi
entrainment coefficient
buoyancy flux (Equation 17)
vorticity
mean temperature
temperature correlation ensemble average
scale length
mass flux (Equation 20)
density
latitude
streamfunction
angular velocity
Subscripts
B
max
o
r
background
maximum
initial or surface
reference height
-------
SECTION 1
INTRODUCTION
A.R.A.P. has been developing a program under EPA sponsorship
since 1971 for practical predictions of the dispersal of pollu-
tants in the atmosphere based on a second-order closure model of
turbulence. This work is detailed in a number of reports and
publications (1-10) and together with the work of others (11-15)
has led to the general recognition that second-order closure
approaches to turbulent modeling have a distinct advantage in
terms of general validity over first-order eddy viscosity ap-
proaches. "Turbulence Modeling and Its Application to Atmos-
pheric Diffusion, Part II" (4) contains a critical review of the
current status of our A.R.A.P. invariant model. It contains a
detailed discussion of the modeled terms, the empirical determi-
nation of the model coefficients, validation tests, and sample
applications. A number of different plume disperal calculations
have been made (1-5).
Section 4 contains a discussion of the modifications imple-
mented within our computer codes to permit computation of buoyant
plume rise. Sections 5 through 8 contain the results of the
major effort in this contract, that of examining the behavior of
buoyant plumes using these codes. Section 5 contains a compari-
son with one-dimensional plumes and jets, showing how well our
predictions compare with results from a careful laboratory ex-
periment. Section 6 shows the results of a calculation of a one-
dimensional buoyant plume released into an ambient stream. In
Section 7, we study the characteristics of a plume rising into a
stable quiescent atmosphere. This computation demonstrates the
-------
interaction of the heated plume with the surrounding stable
atmosphere, and shows the production of Brunt-Vaisala waves from
the plume and the development of the turbulence within the plume
itself.
In Section 8, we discuss plumes rising into stable, neutral,
and unstable moving atmospheres. Here we show the cross-sectional
appearance of the plume, and comment upon plume bifurcation, rise
height, and entrainment.
Section 9 deals with the other major task of our contract:
an examination of the feasibility of second-order closure in
three-dimensional applications. We propose an integral model
that holds the promise of permitting a two-dimensional, unsteady
computer code to compute a three-dimensional, unsteady regional
air quality problem reasonably well.
-------
SECTION 2
CONCLUSIONS
The A.R.A.P. computer model for predicting the fate of non-
chemically reacting pollutants released into the atmospheric
boundary layer has been extended to compute the development of
buoyant plumes.
The sample calculations confirm the laboratory observations
that positive buoyancy increases the entrainment rate into a
plume.
When a plume is released into a stable atmosphere, sample
calculations suggest that approximately 10 percent of the plume's
energy is converted into radiating internal waves.
Under certain conditions, a rising plume from a point source
can divide into two distinct plumes due to the vortex pair gene-
rated by the interaction of the wind and the plume.
Finally, in a separate integral analysis, we demonstrate a
promising means of incorporating some of the capability of
second-order closure modeling in the fully three-dimensional,
unsteady case required for urban air quality models without going
to a costly, fully three-dimensional code.
-------
SECTION 3
RECOMMENDATIONS
The ultimate goal of this research is to permit the full power
of the second-order closure approach to turbulent modeling to be
applied to practical regional air quality models. One approach
is to use individual plume predictions for a number of different
combinations of atmospheric and release conditions to parameterize
required modifications in the Gaussian plume distributions pri-
marily used in air quality models. We believe such individual
plume calculations as we have detailed in this report and in
references 1-5 should be continued to utilize our model's current
capabilities. However, computational efficiency may be achieved
by vertically integrating the basic equations across the boundary
layer. We then reduce the general problem to a two-dimensional,
unsteady calculation for engineering estimates of the boundary-
layer-averaged wind velocities and accompanying pollutant dis-
persal. The potential advantage in both accuracy and computer
time over fully three-dimensional simulations based on first-order
eddy viscosity approaches justifies pursuing the approach further.
-------
SECTION 4
COMPUTER CODE IMPROVEMENTS AND MODIFICATIONS
A major effort of this past year was directed toward en-
hancing the two-dimensional, unsteady code developed over the
last few years for the EPA, the Navy, and NASA to enable us to
accurately compute plumes with significant vertical momentum.
The most significant improvement came after our attempts at
studying two-dimensional planetary boundary layers, particularly
in a coastal environment (16). With a primitive equation code,
using the hydrostatic pressure approximation and solving for the
U and V velocities using their partial differential equations,
we used a quadrature of the continuity equation to determine W ,
and therefore were forced to relax one boundary condition on W .
Ideally, in this problem, we could anticipate that W approaches
zero at the surface and at the top of the computational domain.
Only one of these conditions could be used, and we had to permit
flow through the top of the domain, since we needed W approach-
ing zero at the surface. If we relax the hydrostatic approxi-
mation and include a dynamic equation for W , we are permitted
both boundary conditions on W . However, we also had to include
a loop to solve for the local pressure value. Fortunately, the
solution of the Poisson equation in two dimensions has been
explored carefully, and several reliable direct solution tech-
niques are available for obtaining rapid results (17,18). How-
ever, the satisfaction of continuity (not an a priori assumption)
had to be maintained by the addition of a correction term to the
pressure forcing function. This technique was not always success-
ful.
-------
This problem Is intensified whenever the hydrostatic pres-
sure forms a dominant share of the pressure field. Then two
terms (the vertical pressure gradient and temperature terms) tend
to balance in the vertical momentum equation. The smallness of
the remaining terms which determine W leads to excessive de-
mands on the numerical accuracy.
In order to make the same program operational whether the
hydrostatic approximation is valid or not, we recase our equa-
tions with streamfunction and vorticity as dependent variables.
By defining V = -|| and W = +||- , so that,
W2, f3V 3WA (Fn -n
V \i> — —I-*— — -K—I = —ft ^^q• -Li
\3z 3y/
we obtain an equation for ri that does not contain the pressure.
By discarding the imbalance problem, we do however acquire the
difficulty of having bo specify boundary conditions on i|» and
n rather than V and W .
While we were modifying the code, we decided to gain some
numerical accuracy in our implicit scheme by going to a Crank-
Nicolson form of the tridiagonal matrix elements. This improve-
ment, coupled with the streamfunction-vorticity formulation and
the use of a generalized direct solver (now for i|/ ), have made
significant improvements to our speed of computation (on the
order of a 40 percent increase on our in-house computer) and the
accuracy of the results. The computations discussed in Sections
7 through 9 were made with the enhanced version of the two-dimen-
sional, unsteady code.
One of the major areas of concern, buoyant plume rise in a
stable layer, forced us to study and handle the internal wave
radiating boundary problem. Here a rising buoyant plume will
radiate waves at frequencies that do not exceed the Brunt-Vaisala
frequency
N = CEq. 2)
-------
These waves are constantly produced and tend to travel along the
inversion, eventually reaching the boundary of our computational
domain. In addition, the later stages of plume development see
the propagation of waves of smaller and smaller wavelength. As
the simulation proceeds, the computational mesh spacing becomes
too large to resolve the wave structure. Either one or both of
these problems may cause premature program termination by an
inadequacy in our matching outflow boundary conditions or our
mesh descretization.
At large horizontal distances from the source of the distur-
bances, the streamfunction, vorticity, and temperature each obey a
wave radiation condition of the form
= ° ' (Eq* 3)
where c is the linearized horizontal phase speed of the locally
dominant wave,
f \jj n _ f0 Cartesian
axisymmetric
and we use the local values of ty and n • A stable finite
difference form of Equation 3 is
*n,t - *n,t-l + C M- (*n,t ~ V-l.t-l' = ° (Eq. 5)
Similar conditions also hold for the vorticity and temperature.
Since the dominant waves are inviscid, a simple zero slope at
large x is sufficient for all turbulent correlations. These
conditions appear to permit the smooth passage of Brunt-Vaisala
waves through the outer boundary.
The mesh spacing near the plume centerline, however, is a
more difficult problem. In this regard, we are constrained by a
maximum possible number of mesh points, plus an empirical rule
regarding the ratio of the maximum to mimimum spacing difference.
Our algorithm for determining the streamfunction appears to set
restrictions on the neighbor-to-neighbor spacing, Axn//Axn_i > or
7
-------
on the maximum to minimum spacing, Ax _ /Ax . . The precise
max mm
limits have not yet been determined.
Documented versions of our one-dimensional boundary layer
code and our two-dimensional pollutant dispersal code were deli-
vered early in the contract period to our EPA technical monitor,
A version of the full atmospheric code, containing the improve-
ments discussed above, is currently being tested on the EPA Uni-
vac 1110.
-------
SECTION 5
MODEL COMPARISON WITH TWO-DIMENSIONAL JETS AND PLUMES INTO
A QUIESCENT NEUTRAL FLUID
Reliable atmospheric plume data taken under carefully con-
trolled conditions are difficult to obtain. In this section, we
compare model results with data on two-dimensional jets and
plumes from controlled laboratory flows.
As long as the jet or plume has a parabolic character, it is
possible to compute it using our one-dimensional code modified to
permit a gravity component in the direction counter to the flow.
As was done with the jet flows, we recast the equations into
stream-line coordinates, avoiding the problem of computing the
normal velocity W by defining the independent variable \\i as
|i = U (Eq. 6)
d 2
and write the equations in their two-dimensional, steady version
(marching downstream in x while solving for the one-dimensional
variations in z ) . In these equations we assume that gravity
now acts against the flow, so that in Donaldson's formulation the
gravity vector is written as
K± = (g, 0, 0)
The resulting equations for the mean streaming velocity and tem-
perature, in the large Reynolds number limit, become
__ (Eq. 7)
0 U
-------
1®.
3x
!*§.
9^
(Eq. 8)
The appropriate Reynolds stress equations then become
IT
(Eq-10>
8x
3uu
Ix
H
— 9U
If
) - ^
(Eq. 13)
(Eq. 14)
(Eq. 15)
with the scale length given dynamically by the equation
• 0.35
0.6
0.375U
(Eq. 16)
When g = 0 these equations may be used to compute a two-dimen-
sional jet and the transition to a plume obtained as g is
increased.
10
-------
Our prediction of the flow in a two-dimensional wake has
been presented previously (4) and shows how well the model does
(without the permission of model constant change) at predicting
the deficit velocity, temperature, and turbulent correlation pro-
files in the self-similar region. Similar results for the two-
dimensional jet are shown in Figures 1-4 and compared against
available laboratory data (19-21). We begin our calculation with
Gaussian distributions of velocity, temperature, and turbulence,
and run until similarity exists. In his experiment, Kotsovinos
exhausts heated water from a slot into a large quiescent water
tank. Thermistors and a laser Doppler velocimeter permit him to
determine the flowfield. In all cases, we show his furthest
downstream results only; one of his experimental curves (his
Figure 4.2.4) clearly shows that he has not yet reached simi-
larity. Our numerical predictions show that his last station may
be taken as a reasonable estimate of similarity, however.
Figures 1 and 2 show the normalized cross-sectional profiles
of streaming velocity and temperature. All curves are plotted
normalized by the half-width of the velocity profile (as appro-
priate from the experiment or the prediction). These figures
show that our model does a reasonable job of predicting the mean
profiles, although one would like better agreement at the larger
z distances. These curves will be discussed in greater detail
later.
Figures 3 and 4 give the similarity profiles of the variance
of the velocity in the streaming direction, vuu , and the tempera-
fsasf
ture ve"7 . Both predictions fall within the scatter of the data
and tend to tail out a little high near the edges of the jet.
Our results for these same model constants in the two-dimensional
wake overpredicted the behavior of ef7 . A modification of the
constants to improve the prediction in one case would probably
worsen the prediction in the other. The predicted mean profiles
11
-------
do not fall off as rapidly as the data suggest. The reasons for
this behavior are unclear, but may involve the numerical approach
used to solve for the expanding jet, a technique in streamfunc-
tion coordinates that requires U , > 0 .
H edge
Kotsovinos also studied the two-dimensional plume by heating
his exiting water to a high enough level to permit gravity to
influence the development of the flow. In this way, he was able
to arrive at the entrainment properties of the simple plume; our
comparisons with his results are given in Figures 5-8.
Figure 5 shows the good agreement we obtain with the mean
velocity profile. In fact, in the normalized units plotted here,
our predictions for U for both the jet and the plume are nearly
identical. The one difference not seen from this plot is that
U_._ in both experiments and predictions in the jet are decaying
m.cix -i /Q
as (x/D) ~ , while the plume value reaches a fixed, constant
value of IL,Q /U „ =1.1 . The spread rate Az c/Ax grows at
lUclX UlcLX • _2
a slope of 0.13 for the plume and 0.12 for the jet, both consis-
tent with the experiments.
In Figure 6, we show a disappointing comparison with the
mean temperature measurements. Here we underpredict the plume
profiles, while we tended to overpredict the jet profiles.
Before we begin questioning the model, however, it is best to
note that we are using the .same model constants to predict three
types of two-dimensional flows: jets, wakes, and plumes. Gravity
definitely affects the cross-sectional temperature profile, yet
the similarity between the plume and jet predictions still exists.
A plot of the behavior of the maximum temperature as a function
of downstream distance shows that in the experiment the jet
—1/2
temperature falls asymptotically as 0.28(x/D)~ while the
_i
plume decays as 0.28(x/D) . Our simulation gives 0.3^(x/D)
for the jet, and 0.32(x/D)~ for the plume.
12
-------
Ir
Figure 1.
-2
Comparison of the predicted cross-sectional mean velo-
city profile in a two-dimensional self-semilar jet
with the experimental observations of Kotsovinos (19)
.8
'max
.6
.4
0
-2 -1.6 -1.2 -.8 -.4
0 -.4
.8
1.2 1.6
Figure 2. Comparison of the predicted cross-sectional mean tem-
perature profile in a two-dimensional self-similar jet
with the observations of Kotsovinos (19).
13
-------
Vuu
.4r
.3
U
max
.2
o Kotsovinos
x Davies, et. al.
o o
-2 -1.6 -1.2 -.8 -.4
0
2/1
A .8 1.2 1.6
.5
Figure 3.
Comparison of the predicted cross-sectional profile of
the streaming velocity variance in a two-dimensional
self-similar jet. The observations are from Kotsovinos
(19) and Davies, et al. (20).
.4
.3
®
max
.2
o Kotsovinos
x Davies, et. al.
a Uberoi & Singh
o o
9 x
x a
-2 -1.6 -1.2 -.8 -.4
.4
1.2 1.6
Z/Z
.5
Figure 4.
Comparison of the predicted cross-sectional profile
of the temperature variance in a two-dimensional self-
similar jet. The observations are from Kotsovinos (19)>
Davies, et al. (20), and Uberoi and Singh (21).
14
-------
Ir
-2 -1.6 -1.2 -.8 -.4
.4
1.2 1.6
Figure 5- Comparison of the predicted cross-sectional mean velo-
city profile in a two-dimensional self-similar plume
with the experimental observations of Kotsovinos (19)
.6
"mox
.4
.2
(9
-2 -1.6 -1.2 -.8 -.4
0
Z/Z
.4
.8
1.2 1.6
.5
Figure 6.
Comparison of the predicted cross-sectional mean tempera-
ture profile in a two-dimensional self-similar plume
with the observations of Kotsovinos (19).
15
-------
Figure 7 shows a comparison of the velocity correlation u
uu . It is surprisingly good in light of the mean profiles. The
temperature variance is plotted in Figure 8 for the plume. Now
we see that we overpredict ez for the plume, but underpredict
it for a jet.
An important point of interest is to observe that while the
normalized velocity and temperature profiles from the jet and the
plume appear similar, the profiles of the turbulent correlations
are not. Gravity modifies the jet-like profile behavior into a
shape that is distinctly different.
In a later publication, Kotsovinos and List (22) reworked
some of the original results into nondimensional variables. Two
of these variables are plotted as a function of their nondimen-
sional downstream distance in Figures 9 and 10. In Figure 9 we
show the variation of the buoyancy flux
00
/
g = |_ I U0dz (Eq. 17)
"o
— 00
as a function of distance. (In these figures,
•/
_00
m = I U2dz (Eq. 18)
is a momentum flux of the jet or plume and m denotes the ini-
tial value.) The numerical run was made by beginning a jet cal-
culation with a small amount of gravity. The initial decay of
the normalized buoyancy flux is at the -1/2, power, both for the
experimental data and for the simulation. When the gravity
effect increases, the jet transitions into a plume and the normal-
ized buoyancy flux becomes constant.
They quote an asympototic average value of 0.42 , while we
obtain 0.33 •
16
-------
u
max
-2
Figure 7 •
Comparison of the predicted cross-sectional profile of
the streaming velocity variance in a two-dimensional
self-similar plume with the observations of Kotsovinos
(19).
.5
.4
.3
8'
'max
-0-J-
-2 -1.6 -1.2 -.8 -.4
.8 1.2 1.6
Z/Z
.5
Figure 8.
Comparison of the predicted cross-sectional profile
of the temperature variance in a two-dimensional self-
similar plume with the observations of Kotsovinos (19)
17
-------
10
IO
X.
-------
In Figure 10, we compare the predicted and observed behavior
of the local Richardson number,
19)
_
Ri
Where
J
_00
(Eq> 20)
is the mass flux. The predicted initial growth of Ri is at the
3/2 power, consistent with the experimental observations. The
predicted values typically exceed the experimental data at inter-
mediate distances i but eventually oscillate about a mean quite
close to the observed value of 0.63 . Both of these curves show
a favorable comparison between numerical predictions and experi- .
mental observations .
The biggest, discrepancy between the data of Kotsovinos and
List, and our predictions, is in their plot of the mass-to-
moment um ratio
(Eq. 21)
They demonstrate that C may be taken at a constant value of
0.5^ throughout the jet-to-plume transition. Our simulation is
consistent with a constant, also, but at a higher value of C =
0.70.
An interesting feature of the predictions is the variation •
in the entrainment coefficient ae , for both the jet and the
plume, where ae is defined as
' 22)
with Q the mass flux in the x direction, measured along the
centerline. In Figure 11, we plot ae for both the jet and the
19
-------
plume. The predicted asympototic behavior for the jet is
a = 0.042 , while for the plume we obtain a = 0.115 • These
e e
values may be compared with the observed averaged values quoted
by Kotsovinos, ag = 0.055 for the jet and 0.11 for the plume
over the approximate duration of his experiments. The predicted
similarity value for plume entrainment is within 5 percent of his
average values, but is only within 13 percent for jet entrainment
at x/D = 200 . In their later publication, Kotsovinos and List
(22) note that a constant entrainment coefficient does not de-
scribe their experimental measurements as well as the constant C.
Further comparison may be made with Briggs (23), who cites three-
dimensional values of 0.090 for the jet and 0.125 for the
plume.
The largest uncertainty in the model itself is the edge
boundary condition on the dynamic scale. A fine tuning of this
condition from a value which gives good agreement in the two-
dimensional wake (10) to a 25 percent higher value would raise
the asymptotic value of the jet entrainment to 0.055 • However,
we have not considered such fine tuning as justified at the pre-
sent time.
20
-------
.12
.10
.08
.06
.04
.02
Plume
Jet
40
80
x/D
120
160
200
Figure 11.
Predicted variation of the entrainment coefficient
two-dimensional jets and plumes compared with the
average observed values of Kotsovinos (19).
in
21
-------
SECTION 6
TWO-DIMENSIONAL PLUME IN A MOVING STREAM
The previous section dealt with our prediction (and com-
parison with experiment) of a slot plume exhausted into a quie- •
scent, surrounding fluid. The next level of sophistication,
both in our code simulation and in the laboratory, is to exhaust
a buoyant plume into a moving fluid. Qualitative experimental
results of this effect were performed by Cederwall (24). The
case of particular interest here is the perpendicular exit of a
plume into a moving stream.
This case runs the full spectrum of the plume rise problem.
In its very initial stages the plume exit velocity and turbulence
dominate the problem, and the plume continued rising from its
source without influence from the ambient crosswlnd. Eventually,
however, the flowing stream forces the plume to begin bending
over; it naturally continues to spread during this process, and
would continue to rise (however slightly) unless its vertical
movement is capped by a temperature inversion. In the laboratory
configuration, the plume spreads until it penetrates the entire
streaming depth.
If we confine our attention to a single experimental run -
one shown visually as exiting normal to the freestream Chis run
83) - then all of the necessary data are available to simulate
the experimental setup. In this particular case, the jet was
flowing at 17 times the speed of the freestream cross-flow; the
source Froude number was set to 0.57 (this, with an assumption
on the shape of the exiting jet and the freestream speed, permits
the evaluation of the temperature difference between jet and
22
-------
ambient fluid); and the normalized volume flux was 23 (giving
the needed fluid depth to balance the exiting jet flow volume).
Figure 12 shows the experimental run of Cederwall for this
particular case.
To shorten the computer time to reach steady state, we
confine our attention to the superimposed box region as shown in
Figure 12. Within this region, we assume that the incoming
uniform stream (from the left) enters at V = 4 m/sec , to
interact with an exiting slot jet of the form
W1et = Wmax ^-(W^ (Eq' ^
j " v nidx
with VT = 17 and D = 1 . We assume that the exit tempera-
in 3.x
ture may be represented by a Gaussian distribution with
0 = 1.8°C (consistent with the experimental source Froude
number).
For this simulation we need to assume a jet turbulence
level. From Kotsovino's work in the last section (and Figure
3), we see that the streaming correlation may reach a value of
ww/W = 0.3 . Assuming that this constitutes half of the
TUclX Q Q
turbulence level, we set the maximum q = 50 m2/sec , also in
a Gaussian distribution with D = 1 .
At the jet exit (the lower boundary), the equation for W. ,
may be integrated and differentiated to give the lower boundary
conditions on \i> and n , respectively. The freestream inflow
vorticity is assumed equal to zero, while the inflow ^ is con-
sistent with the uniform inflow. Along the top f is held at
its maximum inflow value, and n = 0 . A zero slope condition
is imposed on y along the outflow boundary. The temperature
and turbulence profiles are likewise assumed known along the
lower and inflow boundaries, and are permitted zero slopes at
the top and downstream boundaries.
23
-------
Our simulation is shown in Figure 13. Qualitative compari-
son is quite good, with the general bending-over character of
the exit slot jet clearly evident.
Figure 12.
The experimental result of Cederwall's run 83 (24)
for a vertical buoyant slot jet exiting into a
cross-stream. The box size is our simulation region,
-------
===++)
..... —===+ + ) )1
- --------- ====++))11
= = = = ------ = = = + +)) 11 ZZ
= = = = = = = --- ===++))1lZZX
============++)11ZZXXA
=+++=======++) 1 1 ZZXAA^
+++++++===+) )11 ZXXAi-1 *W
=+++))))) )))))) 1ZZXAMMMM66
-n-)))1111111111ZZXAMMM€€€€
)11lzl11111lZXXA*M€€€e€€
++++) ))111ZZXXXXXXXXXAMMM€€€€€€€€
))111ZZZXXXAAAAAA A
€€€€
ttt*
e€e€ €€€€<«€<
-=)X!4€€€€t<€€ttt€tftfff€f€ftffftfeffttCtCtfCtCttCttg
-)< €€«ttttttttfttuttttf tttttteeftttttttt
t f««f«t«t €€€€€€€€
-+X6€« I liflftftftlitfttf II tlflt «««€«€€«««€€€ 6e6€€€MMM
--)€<€tIIfHHtffft«««€t««)+ + =. = =
=A«mif «€€M«AXXXZZZZZZZZZZZZZZZZZ11 111))) + + + = = =
-=X€li§lt€Z=-
Figure 13. An intensity plot of the variation of jet temperature
within the box defined in Figure 12 for the Cederwall
experiment. Maximum intensity corresponds to maximum
temperature.
25
-------
SECTION 7
BUOYANT PLUME RISE INTO A STABLE QUIESCENT ATMOSPHERE
We next examine a two-dimensional, unsteady calculation of
buoyant plume rise from a heated surface into a quiescent atmos-
phere. Here, we wish to study not only the time-development of
the plume and the periodic character of the wave structure out-
side the plume, but also the effective interaction of the plume
with the surface layer.
For this calculation, we use an axisymmetric version of the
two-dimensional code developed under contract to the Nuclear
Regulatory Commission to study strongly swirling flows. In the
present case we are not interested in the swirl, but only in the
radial coordinate system. We represent the stable atmosphere by
a constant background temperature gradient 90B/9z of 0.003°C/m ,
Along the lower surface we assume that the flow is sufficiently
unstable so that the surface will be free-convection-like (8) at
a reference height of z = 1 m above a hydrodynamic roughness
height of z = 0.1 m . In this case, we may set fy = 0 at
z . The vorticity carries a compatibility relationship to the
streamfunction through the Poisson equation, which reduces to the
form 2^
_ +
n F 2z in I
where iK is the streamfunction value at the point above z
The presence of a heated surface is felt by assuming that
the temperature at z is a simple Gaussian distribution of 3°C
26
-------
maximum value and 100 m spread. The wide spread is suggestive of
cooling tower dimensions; a surface temperature assumed no larger
than the stable temperature at the top of the computational
domain (at 1000 m) ensures that the developing plume will be
capped within our field of interest. The temperature at z is
then given by formula
(Z
0 =0 + (0 -0 )
°r °s <°+ °a> (z_
where 9+ is the temperature at z+ > z , and e is the sur-
face temperature.
We anticipate that as the plume develops, it will entrain
fluid from along the surface and accelerate it into the unstable
layer near the origin. The presence of the ambient stable tem-
perature gradient permits the development of Brunt-Vaisala waves
above the plume. These waves will interact with the far-field
boundaries, and will also lead to the radiation of waves of
smaller and smaller wavelengths. Our radiating boundary condi-
tion, discussed in Section 5, will handle the far-field condi-
tions; however, the short wavelength waves will eventually become
too small to be defined by the mesh spacing near r = 0 .
An illustration of the plume buildup and oscillation may be
observed in Figure 14 where we have plotted a near-edge coutour
value of the turbulent kinetic energy q2 , its 0.1 m2/s2
value. Starting from nominally zero values, the plume had reach-
ed an altitude of nearly 200 meters and a radial spread of 100
meters by a time of 300 seconds. The continued upward movement
of the plume is arrested near t = 600 seconds, after which time
the top of it drops. At this time, the first internal gravity
wave is generated. Further time slots show the movement of the
top of the plume and the regular (almost steady) profile to the
plume below 100 meters. The top of the plume is located near
27
-------
330 meters. The rise height in a calm, stable atmosphere may be
estimated from a formula given by Briggs (23)
Az = 5.0 F1/4 N~3/i| (Eq. 26)
where the buoyancy flux is given by
/ T T/"\ I A \ J A
.. 27)
Evaluating these formulas at the reference height, we find that
the rise height of the plume should be about 300 meters. Thus,
it would appear as though the early time behavior of the simula-
tion is quite consistent with available observations.
The nature of the radiating wave patterns may be seen by
presenting the streamline patterns at several times (Figure 15).
Figure 15a gives the pattern at t - 300 seconds, when the simple
structure is a vortex ring (i.e., at one-half of a Brunt-Vaisala
period, the first internal gravity wave has been generated). By
t = 600, Figure 15b, the ring is flattened slightly by a second
ring (or wave) developing near z = 300 meters (where the top of
the plume exists). By t = 900, the ring in Figure 15b has envel-
oped a large portion of the domain of observation, and a second
ring of opposite sign has appeared above the plume. A new cell
is just beginning near r = 0. The structure near z = 300 meters
is quite complicated and exibits two vertical stagnation points
at 300 and 350 meters. Subsequent frames show the continued
pattern development. The streamline cell structure is now clearly
visible, and we observe that the wave structure does not appear
to have any problems in passing through the large radial boundary.
This is a sample of a plume problem where the application of
a simple turbulent entrainment coefficient would be inappropriate.
The low-level convergence due to the buoyant acceleration of the
plume, and the upper-level divergence due to the negative buoyancy,
28
-------
600t
2100
r, m
r, m
r, m
r, m
0 100 0 100
r, m r, m
Figure 14.
Prediction of the development of a rising plume
into a quiescent stable atmosphere. Shown here
are the contour lines at several times after flow
initialization for the q2 = 0.1 m2/s2 value of
turbulent kinetic energy.
29
-------
. t-300 sec
t"600 sec
C. t-900 sec
t-1200 sec
g m t-1500 sec
f. t-1800 sec
Figure 15-
Prediction of the development of a rising plume into
a quiescent stable atmosphere. Shown here are various
streamfunction contours at several times after flow
initialization.
30
-------
are intimately connected with the turbulence, so that a simple
turbulent entrainment cannot be separated from the rest of the
flow mechanisms operating in the problem.
The rising plume establishes a complicated cross-sectional
flow pattern into which the stable background temperature gradient
permits the generation of Brunt-Vaisala waves. Surface heat flux
drives the plume development, while the Brunt-Vaisala waves carry
some energy from the plume (turbulent dissipation handles the rest;
To estimate the significance of the waves as energy carriers, we
compare the ratio of energy lost along the outflow boundary to
the energy inflow at the surface. It may be shown that by linear-
izing the governing equations in the large r region, the average
energy flux through the boundary may be written as
t
out o max
evaluated at r = r , where 1 is the wave number of the out-
ma x
going waves, {Nz2/kr2} is the group velocity, and <•> denotes
the average over a wave period. This quantity is compared with
the total turbulent heat flux coming through the surface
rmax
/
0
"in ' <*>•
evaluated at z = z . Using the results for t > 1200 seconds,
we find that the energy carried away from the vicinity of the
plume by the waves may amount to 10 percent of the entering
energy. Over 70 percent of the energy is dissipated by turbulent
kinetic energy processes within this same region. At larger
radii, energy is lost from the waves to thermal dissipation of the
potential energy.
31
-------
SECTION 8
PLUME RISE INTO MOVING ATMOSPHERES
We now turn our attention to the corresponding plume rise
problem in a moving atmosphere. We use the two-dimensional, un-
steady code, marching in x downwind of the plume release point.
This technique forces us to begin computation after the buoyant
plume has become sufficiently bent-over that (W/U )2 « 1 . All
CO
of our calculations begin with a Gaussian distribution of the
heated air. Our simulation tracks its behavior downstream, as it
forces the formation of a vortex pair as it rises, entraining
surrounding air, and producing turbulence.
A straightforward application of our work with aircraft vor-
tex wakes leads to the solution for plumes rising in a neutral
atmosphere. These results were presented at the AMS Third Sym-
posium on Atmospheric Turbulence, Diffusion, and Air Quality (2).
Our discussion will center primarily upon the behavior of the
plumes emitted into moving stable atmospheres.
We again may begin with a Gaussian distribution of tempera-
ture at the assumed point of plume emission. In the initial
development of the plume, we do not anticipate any appreciable inter-
actions with the surface; thus, we use a simple no-slip boundary
condition there. At large y distances, we use the two-dimenional
analogue to Equation 5 for the radiating wave condition. We
again begin with a background temperature gradient of
= 0 . 003°C/meter . For this simulation we are able to
^
B
calculate the transport of a passive species, as well as the
temperature and velocity components, and will use the tracer pro-
32
-------
file behavior to illustrate the problem simulation.
Beginning with a simple Gaussian distribution for the tracer,
we show in Figure 16 its later downstream development. ^igure l6a
shows the contour shapes x = 500 meters after initialization.
The simple Gaussian structure with maximum value located at z =
200 meters and a spread of 65 meters has risen above z = 400
meters by this frame. The top of the plume is flattening out as
the stable layer is penetrated. At x = 1000 meters. Figure l6b,
the plume has continued diffusing upward and outward into the
stable layer. Figure l6c shows that the plume is in a downward
oscillation by x = 1500 meters, while the final frame, l6d, shows
that the plume at x = 2000 meters is spreading laterally along
the stable inversion layer. Cross-sectional profiles of the
streamlines are very similar to those produced by the quiescent
plume rise.
A formula from Briggs (23) for computing plume rise height
for bent-over plumes in an atmosphere of constant stratification
is
A—}
\U N2/
00
Az = 2.6 -S ) (Eq. 30)
N:
where F has been defined in Equation 27. Substituting into
Equation 30 with the data from this simulation, we find that
Az ^ 260 meters. Inspecting the 2 percent contour on the
tracer, we find that our simulation gives Az ^ 300 meters.
We may compare the rise height of the center of the plume
rising into stable layers with the rise heights for neutral and
unstable layers as shown in Figure 17. Normalization by the
Froude number
~1/2 (Eq. 3D
33
-------
permits the collapse of all buoyant plumes released Into neutral
surroundings to the single solid curve shown, with the asymptote
as given. Unstable layers accelerate the plume rise, while a
stable layer will stop the plume rise. The oscillations about
its inversion level are the Brunt-Vaisala waves.
-------
800
700
Cma«
= 049
800,
100
a. x = 500 meters
= 0.25
too Soo boo '400 boo
y, m
c. x = 1500 meters
800 T
700.
— =0 30
800,
700
600
500. .
z, m
4OO
3OO.
y, m
b. x = 1000 meters
— =0.21
d. x = 2000 meters
Figure 16. Prediction of the downstream development of a moving
buoyant plume Into a stable atmosphere. Shown here
are various passive tracer (smoke) contours at several
distances downwind of the Initialization plane.
35
-------
dz
3
D
lOOr—
z/D
— =0 , all Fr numbers
) asympto
- = -l, Fr = l
2/3
l.8(x/Fr-D) asymptote
= .3,Fr = .84
Figure 17.
100
Buoyant plume centerline rise height as a function
of distance downstream of release point. Curves
show the effects of unstable, neutral, and stable
background temperature gradients. A0max is the
maximum difference between the initial plume
temperature and the local ambient atmospheric
temperature.
36
-------
SECTION 9
AN EXAMINATION OP AN INTEGRAL APPROACH TO USE IN REGIONAL
AIR QUALITY MODELS
One of the tasks in the current contract is to identify how
second-order closure may best be utilized in practical urban or
regional air quality models involving complex terrain and numerous
plumes. One way to utilize a sophisticated model is to use in-
dividual plume predictions for a number of different combinations
of atmospheric and chimney exit conditions to parametrize modifi-
cations in the Gaussian plume distributions presently used in air
quality models. We have made such calculations (1) and plan to
continue this type of calculation in the future. However, it
seems probable that in many situations the full influence of com-
plicated terrain or unsteadiness cannot be parameterized on the
basis of the individual steady plume calculations. It is pos-
sible to consider using the full model for direct three-dimen-
sional, unsteady numerical calculations, but the enormous com-
puting load this entails makes it highly desirable to seek some
practical approximations.
We have experimented with vertically integrating the basic
equations across the boundary layer to remove one dimension from
the problem. This will reduce the general problem to a two-
dimensional, unsteady calculation for the boundary-layer-averaged
variables. The integration process is exact, but approximations
must be introduced to parameterize the relationships between the
different integrals which occur in order to close the set of
equations. Results of a simple test of this approach indicate
that it may be the best approach to follow to obtain a practical
model for air quality studies.
37
-------
To see how this process may work, we consider the limiting
case of a neutral planetary boundary layer. The equations for
the mean velocities and second-layer correlations as given in
reference 4 may be integrated to yield
(wu) . _r . ^ (Eq. 32)
o o
u> + (Eq. 33,
"v" sin 4> - uĄ cos <|>) - 2uĄ -
= <-4fluv sin
c
['•
D _ ,,,^ cos ^ + ,Aq gw, _ ^^ . u, + Ł^_ > (Eq. 36)
Vc
38
(Eq. 34)
Dt
„,_ „ , ,>,„ (Eq' 35)
+ v, ' rtv" '
Dt c L --jo
gy
= <2fl(vv sin t() - vw cos - uu sin 4>) - uw •*—
Dt ,_ -,„ dZ (Eq. 37)
— STT ^
l_ —JU
— .. — 3U
sin d> - ww cos 4> + uu cos 4>) - ww •%— >
(Eq. 38)
sin ij> + uv cos ) - vm ->
(Eq. 39)
-------
where
Dt -- 3t~ 9x ~ i^ --
To make this a closed system, we must relate such integrals
_ g TJ _
as to , to , etc. This may
a Z
be done by parameterizing the profile shape of the variable with
respect to z . These profile shapes may contain only one in-
dependent parameter for each independent variable. Of course,
the profile distributions may be chosen to be consistent with as
many boundary condition constraints as desired.
To be specific, we will work through a limiting case for a
sample choice of profile distributions. To achieve the desired
asymptotic approach to geostrophic conditions, we use an exponen-
tial z dependence. At the surface we force the profiles to be
compatible with the surface layer conditions obtained by requiring
that Equations 32-39 also be satisfied across the thin constant
flux layer next to the surface. For neutral flow these surface
conditions reduce to the familiar law-of-the-wall conditions
(U2 + Vp = ln(zr/z0)
r0 40)
V vw
(Eq.
r uwr
with
and the turbulent correlations related by their super-equilibrium
values (7)
qj = /32u*2 (Eq. 43)
wĄ = q/4 (Eq. 44)
39
-------
To these constraints we add matching conditions at z^ , the
upper edge of the surface layer
9U_ I =
Hr "
u! E _ (Eq. 45)
kz (Tj2 + v2)l/2
1 ^ r 'r'
=0 ' (Eq. 46)
Jr
(Eq.
- v
For purposes of this sample calculation, we simplify the
correlation equations, Equations 34-39 , by assuming ww = q2/4
and carrying only two integrated equations, one for q2" and Equa-
tion 38 for uw . With the exception of the Coriolis terms,
Equation 39 for vw is similar to Equation 38 for uw . In the
present approximation we neglect the Coriolis terms in the shear
stress equations so it is consistent to carry only one integral
equation for the shear stress. The equation for q2 is obtained
by adding Equations 34 through 36
> (Eq.
Dt
Profile shapes with a total of 10 independent parameters may
now be used to approximate the five variables (U , V , uw , vw ,
q). We approximate the profiles as:
j _ u = (U - U )e~A^ (Eq. 50)
g g r
V = Vr(l + alZ)e~Xz (Eq. 5D
uw = uw (1 + a9z)e~Xz (Eq. 52)
40
-------
vw = vwr(l + a.,z)e~Az (Łq- 53)
q = qr(l + ai|z)e-Az (Eq< 54)
where z is measured from z
r '
Equations 50-54 yield a closed system when the model scale
A is determined. In a general scheme we would want to utilize
an integral of the scale equation, but for this simple test we
choose to use an expression related to the distance from the
ground and the spread of q ; thus we choose
A = a(ze"3z + zr) (Eq. 55)
The parameter g is chosen so that the relationship between A
max
and the height at which q falls to one-half its maximum value
remains constant at 0.21, approximately the same value which
would be maintained for the dynamic scale equation in the case of
a neutral, steady state boundary layer (25).
Equations 32, 33, 38, and 40-55 form a complete set that may
be used to compute the approximate variables as a function of
x, y , and t . In the limit of steady, spatially homogeneous
flow, the system reduces to a set of algebraic equations to
determine the profile shapes as a function of Rossby number,
Ro = U /fz . Although the profile shapes are only rough approxi-
mations to their true shapes, the integrated variables should be
reasonably represented as functions of Ro . A good way to view
the results is as plots of u*/U and tan" (V/U) versus Ro .
o
The approximate integral results compare with the complete
numerical integration of our equations as shown in Figures 18 and
19. The u* variation is quite reasonable but the surface angle
is a little rougher than we would like. This can be adjusted by
observing from the exact solutions that the height where q falls
-------
to 10 percent of its maximum value is significantly higher than
the height where (U -U) falls to 10 percent of its maximum value.
CJ
If the adjustment is made in the assumed profile shape for z so
that it is proportional to e~ * Xz (in Equation 55), then as
seen in Figure 19, tan" (V/U) agrees much better with the exact
solution.
We could, of course, continue this process of modifying the
profile choices to increas'e the accuracy of the approximation,
but predicting the integrated variables to within about 10 percent
of their exact values (as done by these approximate profile
choices) would be adequate for engineering calculations. Note
from Equations 32 and 33 that predicting u* and surface wind
shear angle correspond to predicting the boundary layer averaged
wind velocities directly for this limiting case.
These sample calculations suggest that a method based on
generalizing the above approach to include stability variations
and pollution dispersal should be a viable procedure for incor-
porating second-order closure modeling in practical air quality
models.
The model by Busch, Chang, and Anthes (26) is representative
of currently available boundary layer models for use with regional
air quality models. They suggest a grid with 20 vertical points.
Thus, in spite of the fact that our second-order closure model
deals with two to three times as many basic equations, our bound-
ary layer averaged model should run computationally several times
faster than theirs. Thus, it is quite conceivable that our more
sophisticated approach may yield an advantage in both accuracy
and computer time. We believe the potential payoff, at least,
justifies pursuing the approach further.
-------
24
20. .
16
z, m
12
4.
30
-20
-10
y, m
Figure 18. Predicted temperature levels for continous automobile
traffic in a surface wind normal to the roadway.
Normalization is by local maximum value.
24,
y, m
Figure 19. Predicted temperature levels for continous automobile
traffic in a surface wind parallel to the roadway.
Normalization is by local maximum value.
-------
REFERENCES
1. Lewellen, W.S., and M.E. Teske. Second-Order Closure Model-
ing of Diffusion in the Atmospheric Boundary Layer. Boundary-
Layer Meteor., 10:69-90, 1976.
2. Teske, M.E., and W.S. Lewellen. Example Calculations of At-
Atmospheric Dispersion Using Second-Order Closure Modeling.
In: Proceedings of the Third Symposium on Atmospheric
Turbulence, Diffusion, and Air Quality, Amer. Meteor. Soc.,
Raleigh, North Carolina, Oct. 19-22, 1976. pp. 149-154.
3- Lewellen, W.S., and M.E. Teske. Atmospheric Pollutant Dis-
persion Using Second-Order Closure Modeling of the Turbulence.
In: Proceedings of the EPA Conference on Environmental
Modeling and Simulation, Research Triangle Park, North
Carolina, April 19-22, 1976. EPA-600/9-76-106, U.S.
Environmental Protection Agency, Research Triangle Park,
North Carolina, pp. 714-718.
4. Lewellen, W.S., and M.E. Teske. Turbulence Modeling and Its
Application to Atmospheric Diffusion, Parts I and II. EPA-
600/4-75-0166, U.S. Environmental Protection Agency, Research
Triangle Park, North Carolina, 1975.
5. Lewellen, W.S., M.E. Teske, R.M. Contiliano, G.R. Hilst, and
C.dup. Donaldson. Invariant Modeling of Turbulent Diffusion
in the Planetary Boundary Layer. EPA-650/4-7.4-035, U.S.
Environmental Protection Agency, Research Triangle Park,
North Carolina, 1974.
6. Donaldson, C. duP. Calculation of Turbulent Shear Flows for
Atmospheric and Vortex Motions. AIAA Journal, 6:4-12, 1972.
7. Donaldson, C. duP. Atmospheric Turbulence and the Dispersal
of Atmospheric Pollutants. In: Proceedings of the Workshop
on Micrometeorology, D.A. Haugen, ed., Amer. Meteor. Soc.,
Boston, Massachusetts, August 1972. pp. 313-390.
8. Lewellen, W.S., and M.E. Teske. Prediction of the Monin-
Obukhov Similarity Functions from an Invariant Model of
Turbulence. J. Atmos. Sci., 30:1340-1345, 1973-
44
-------
9. Lewellen, W.S, M.E. Teske, and C.duP. Donaldson. Turbulence
Model of Diurnal Variations in the Planetary Boundary Layer.
In: Proceedings of the 24th Heat Transfer and Fluid Mechanics
Institute, L.R. Davis and R.E. Wilson, eds., Corvallis,
Oregon, June 1974. Stanford University Press, Stanford,
California, 1974. pp. 301-319.
10. Lewellen, W.S., M.E. Teske, and C.duP. Donaldson. Variable
Density Plows Computed by a Second-Order Closure Description
of Turbulence. AIAA Journal, 14:382-387, 1976.
11. Launder, B.E., and D.B. Spalding. Mathematical Models of
Turbulence. Academic Press, New York, New York, 1972.
12. Lumley, J.L., and B. Khajeh-Nouri. Computational Modeling
of Turbulent Transport. Advances in Geophysics, Volume 18A:
Academic Press, New York, New York, 1974.
13. Mellor, G.L., and T. Yamada. A Hierarchy of Turbulence
Closure Models for Planetary Boundary Layers. J. Atmos.
Sci., 31:1791-1806, 1974.
14. Wyngaard, J.C., and O.K. Cote. The Evolution of a Convective
Planetary Boundary-Layer - A Higher-Order Closure Model
Study. Boundary-Layer Meteor., 7:289-308, 1974.
15. Zeman, 0., and J.L. Lumley. Modeling Buoyancy Driven Mixed
Layers. J. Atmos. Sci, 33:1974-1988, 1976.
16. Lewellen, W.S., and M.E. Teske. A Second-Order Closure
Model of Turbulent Transport in the Coastal Planetary
Boundary Layer. In: Proceedings of the Conference on
Coastal Meteorology, Amer. Meteor. Soc., Virginia Beach,
Virginia, Sept. 1976. pp. 118-123.
17. Buneman, 0. A Compact Non-Iterative Poisson Solver. SUIPR
Kept. No. 294. Inst. for Plasma Research, Stanford Univer-
sity, Stanford, California, 1969.
18. Swarztrauber, P., and R. Sweet. Efficient Fortran Sub-
programs for the Solution of Elliptic Partial Differential
Equations. Rept. No. NCAR-TN/IA-109. NTIS, Springfield,
Virginia, 1975-
19 Kotsovinos, N.E. A Study of the Entrainment and Turbulence
in a Plane Buoyant Jet. Rept. No. KH-R-32. Calif. Inst. of
Tech., Pasadena, California, 1975.
-------
20. Davies, A.E., J.P. Keffer, and W.D. Baines. Spread of a
Heated Plane Turbulent Jet. Physics of Fluids, 18:770-775,
1975-
21. Uberoi, N.S., and P.I. Singh. Turbulent Mixing in a Two-
Dimensional Jet. Physics of Fluids, 18:764-769, 1975.
22. Kotsovinos, N.E., and E.J. List. Plane Turbulent Buoyant
Jets. Part 1. Integral Properties. J. of Fluid Mech., 81:
25-44, 1977.
23. Briggs, G.A. Plume Rise Predictions. In: Lectures in Air
Pollution and Environmental Impact Analysis, D.A. Haugen,
ed., Amer. Meteor. Soc., Boston, Massachussetts, Sept. 1975-
pp. 59-111.
24. Cederwall, K. Buoyant Slot Jets into Stagnant or Flowing
Environments. Kept. No. KH-R-25- Calif. Inst. of Tech.,
Pasadena, California, 1971.
25. Lewellen, W.S., and G.G. Williamson. Wind Shear and Turbu-
lence Around Airports. Kept. No. NASA-CR-2752. NTIS,
Springfield, Virginia, 1976.
26. Busch, N.E., S.W. Chang, and R.A. Anthes. A Multilevel
Model of the Planetary Boundary Layer Suitable for Use with
Mesoscale Dynamic Models. J. of Appl. Meteor., 15:909-919,
1976.
46
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
1 RE
2.
3. RECIPIENT'S ACCESSION-NO.
4. TITLE AND SUBTITLE
TURBULENCE MODELING APPLIED TO BUOYANT PLUMES
5. REPORT DATE
August 1978
6. PERFORMING ORGANIZATION CODE
7. AUTHOR(S)
M.L. Teske, W.S. Lewellen, and H.S. Segur
8. PERFORMING ORGANIZATION REPORT NO,
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Aeronautical Research Associates of Princeton, Inc.
50 Washington Road
Princeton, New Jersey 08540
10. PROGRAM ELEMENT NO.
1AA601 CA-31 (FY-78)
11. CONTRACT/GRANT NO.
68-02-2285
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Sciences Research Laboratory - RTF, NC
Office of Research and Development
U.S. Environmental Protection Agency
Research Triangle Park, North Carolina 27711
13. TYPE OF REPORT AND PERIOD COVERED
Interim 6/76 - 6/77
14. SPONSORING AGENCY CODE
EPA/600/09
15. SUPPLEMENTARY NOTES
16. ABSTRACT
A viable computer model was developed that is based on second-order closure
of the turbulent correlation equations for predicting the fate of nonchemically
reacting contaminants released in the atmospheric boundary layer. The invariant
turbulence model discussed in previous reports has been extended to compute the
development of buoyant plumes. Numerical program capability has been extended by
improving the speed and accuracy of the two-dimensional unsteady turbulent flow
calculation. Plume calculations are made for buoyant plumes rising into a stable
quiescent atmosphere and stable, neutral,,and unstable moving atmospheres. An
examination of the application of an integral approach to our turbulent boundary
layer model is also included.
7.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lDENTIFIERS/OPEN ENDED TERMS
c. COSATI Field/Group
*Air pollution
*Meteorology
*Turbulence
*Turbulent flow
*Atmospheric diffusion
^Mathematical Models
13B
04B
20D
04A
12A
8. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19. SECURITY CLASS (ThisReport)
UNCLASSIFIED
21. NO. OF PAGES
20. SECURITY CLASS (Thispage)
UNCLASSIFIED
22. PRICE
_52_
EPA Form 2220-1 (9-73)
47
-------
"C3 -< d) XT
o o ^ o
^ C O C
w
a
o
CO
i
o
m
o
a: S ? 5
rr, ;; 3 „
^ - 3
»i r* o o
i—1° ~ o
o
rt>
CD
Cr •
g, g
§: *
5 S-
% §
r
O
TJ
-n
U
0
~%
C
z
H
<
m
j*
s»
T)
1 —
1
0
-<
m
~o
m
7
f^
>
r~
H
~n
0
~n
A)
-D
Zl
<
>
H
m
C
0)
m
to
\jf
u
o
o
o
~n
-n
C>
>
r
C
w
z
m
(Si
(A
a
D
a
3
D
m
-
O
IT
6
.&.
01
tV)
0)
00
m
13
<
O
n
3
CD
•3
CO
-TJ
~M
CD
C/)
CD
Q)
— ^
O
i
^
— h
O
— t
co
*"*"
0
O
CD
3
g
o
0)
o
3J
CD
CD
CD
o
zr
Ł11
13
d
D
(D
<
CD
o"
13
3
CD
n
m
Z
<
33
O
^
m
H
>
r~
Tl
JD
O
_!
m
o
—I
O
'Z.
>
O
C
U)
rn
Z
o o
z w
>
z
H
CO
ffi
m
J>
Z
O
O m
H m
O M
a
m
Z
n
------- |