EPA-600/4-75-016a
December 1975
Environmental Monitoring Series
MODELING AND
APPLICATION TO ATMOSPHERIC DIFFUSION
Parti
Environmental Sciences Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Research Triangle Park, N.C. 27711
-------
RESEARCH REPORTING SERIES
Research reports of the Office of Research and Development,
U.S. Environmental Protection Agency, have been grouped into
five series. These five broad categories were established to
facilitate further development and application of environmental
technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in
related fields. The five series are:
4
1. Environmental Health Effects Research
2. Environmental Protection Technology
3. Ecological Research
4. Environmental Monitoring
5. Socioeconomic Environmental Studies
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 Information Service, Springfield, Virginia 22161.
-------
EPA-600/4-75-016a
December 1975
TURBULENCE MODELING AND ITS APPLICATION TO ATMOSPHERIC DIFFUSION
PART I: RECENT PROGRAM DEVELOPMENT, VERIFICATION, AND APPLICATION
by
W. S. Lewellen and M. Teske
Aeronautical Research Associates of Princeton, Inc,
Princeton, New Jersey 08540
Contract No. 68-02-1310
Project Officer
Kenneth L. Calder
Meteorology and Assessment Division
Environmental Sciences Research Laboratory
Research Triangle Park, North Carolina 27711
U.S. ENVIRONMENTAL PROTECTION AGENCY
OFFICE OF RESEARCH AND DEVELOPMENT
ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
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.
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
1. REPORT NO.
EPA-600/4-75-016a
3. RECIPIENT'S ACCESSIOWNO.
4. TITLE ANDSUBTITLE
TURBULENCE MODELING AND ITS APPLICATION TO ATMOSPHERIC
DIFFUSION. PART I: RECENT PROGRAM DEVELOPMENT,
VERIFICATION, AND APPLICATION
5. REPORT DATE
December 1975
6. PERFORMING ORGANIZATION CODE
7. AUTHOR(S)
W. S. Lewellen and M. Teske
8. PERFORMING ORGANIZATION REPORT NO,
A.R.A.P. Report 254, Part I
9. PERFORMING ORG \NIZATION NAME AND ADDRESS
Aeronautical Research Associates of Princeton,
50 Washington Road
Princeton, New Jersey 08540
Inc.
10. PROGRAM ELEMENT NO.
1AA009
11. CONTRACT/GRANT NO.
EPA 68-02-1310
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Sciences Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Research Triangle Park. North Carolina 27711
13. TYPE OF REPORT AND PERIOD COVERED
Interim
14. SPONSORING AGENCY CODE
EPA-ORD
15. SUPPLEMENTARY NOTES
Issued as Part I of 2 Parts
16. ABSTRACT
This report details the progress made at A.R.A.P. during fiscal year 1975
towards the goal of developing a viable computer model based on second-order
closures of the turbulent correlation equations for predicting the fate of
nonchemically reacting contaminants released in the atmospheric boundary layer.
The invariant turbulent model discussed in previous reports has been further
verified by new comparisons between model predictions and experimental
observations. Model capability has been extended by increasing the dimensions
of the program to permit the calculation of two-dimensional, unsteady turbulent
flow. A number of practical plume calculations were made and compared with
standard Gaussian plume assumptions. Calculations are also made for comparison
with the complementary turbulent models being developed at The Pennsylvania
State University and Flow Research, Inc.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lDENTIFIERS/OPEN ENDED TERMS
c. COSATI Field/Group
* Turbulence
* Turbulent flow
* Turbulent diffusion
* Atmospheric diffusion
* Boundary layer
* Air pollution
* Mathematical models
12A
20D
04A
13B
14G
8. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19. SECURITY CLASS (ThisReport)
UNCLASSIFIED
21. NO. OF PAGES
86
20. SECURITY CLASS
URITY CLASS (Thispage)
UNCLASSIFIED
22. PRICE
EPA Form 2220-1 (9-73)
-------
ABSTRACT
This report details the progress made at A.R.A.P. during
fiscal year 1975 towards the goal of developing a viable
computer model 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 turbulent model discussed in previous
reports has been further verified by new comparisons between
model predictions and experimental observations. Model capa-
bility has been extended by increasing the dimensions of the
program to permit the calculation of two-dimensional unsteady
turbulent flow. A number of practical plume calculations were
made and compared with standard Gaussian plume assumptions.
Calculations were also made for comparison with the comple-
mentary turbulent models being developed at The Pennsylvania
State University and at Plow Research, Inc. of Boston, Mass.
Part II of this report contains a critical review of the
current status of our A.R.A.P. invariant model, including a
detailed discussion of the modeled terms, the empirical
determination of the model coefficients, the validation tests,
and sample applications.
This report was submitted in fulfillment of Project
Number 21AUS-12, Contract Number EPA-68-02-1310 by Aeronauti-
cal Research Associates of Princeton, Inc. under the sponsor-
ship of the Environmental Protection Agency. Work was
completed as of August 1975.
ii
-------
TABLE OF CONTENTS
Page
Abstract ii
List of Figures iii
I. Conclusions 1
II Recommendations 3
III. Introduction 5
IV. Program Development 8
Introduction 8
Transfer of Program to EPA UNIVAC 1110 8
Development of the Two-Dimensional
Unsteady Model 8
Numerical Test Computations 13
V. New Verification of the Basic Turbulence
Model 17
VI. Dispersion Estimates Compared with Gaussian
Plume Models 23
Introduction 23
Ambient Turbulence Field 23
Summary of Diffusion Model 28
Line Source Release at Ground Level 30
Point Source Ground Release 38
Elevated Line Release 44
VII. Comparison with Turbulence Models of PSU
and FRI 51
VIII. Comments on "Some Views of Modeling Dispersion
and Vertical Flux" by Pasquill and Smith 58
IX. References 60
Nomenclature 64
Appendix A 66
Appendix B 71
iii
-------
LIST OF FIGURES
Figure No. Page No,
1 Definition sketch for finite difference
scheme 11
2 Comparison of the turbulent velocity fluctua-
tions in the vicinity of a change of roughness
as computed by the steady-state CORE program
( - ) and by the new PEL program ( --- ) 15
3 Comparison of the concentration isopleths for
a line source located at x = 0, z = 0 for
the turbulent fields shown in Figure 2 16
4 Concentration isopleths at x = 0.44zj_U/w*
downwind of a continuous point source release
into a free convective layer as simulated in
a laboratory test by Deardorff and Willis28
( - ) and as predicted by our model ( --- ).
Numbers on the contours are in nondimensional
units of z2UC/Q . 19
5 Centerline concentration isopleths downwind
of a continuous point source release into a
free convective layer as predicted by our
model; notation as in Figure 4 21
6 Centerline concentration isopleths downwind
of a continuous point source release into a
free convective layer as simulated in a
laboratory test by Deardorff and Willis;28
notation as in Figure 4 22
7 Velocity and turbulence distribution pre-
dicted for a neutral, steady-state planetary
boundary layer for a Rossby number
Ug/z0f =10* 24
8 Ratio of the vertical velocity variance to
10% of the wind at 10 m altitude for three
different values of Rossby number 26
>s
9 Correspondence between Ri -and Pasquill's
stability parameter P or stability classes
A-F 2?
10 Ratio of the vertical velocity variance to
u for different values of Ri and mixed-
layer depths for Ro = 106 29
11 The vertical spread az as a function of
distance x downwind of a line source ground
release in a neutral atmosphere (UK = 10 m/sec,
, f = lO^sec"1)
z0 = 10 cm, f = lO^sec"1) 31
iv
-------
Figure No. Page No.
12 Comparison of species distribution with a
Gaussian distribution for the plume of
Figure 11. Normalized profiles are shown
for several distances downstream of the
release point. The numbers give the
exponent in tens of meters; i.e., 0 is the
Gaussian distribution assumed for x = 0,
2 is the profile at x = 100, etc. 32
13 Normalized surface species concentration as
a function of downwind distance ( )
compared with the Gaussian solution ( )
for a wind assumed constant at its values
at a height of 10 meters 34
14 Comparison of vertical spread in a neutral
atmosphere from a line source ground release
for changes in the Rossby number based on
surface roughness zo 35
15 Comparison of vertical spread from a line
source ground release for different values
of Ri for Ro = 106 (z±f/u = 0.451 for
the unstable case) 36
16 Repeat of Figure 15 in unnormalized form
for Ug/f = 100 km for comparison with
Pasquill's Figure 6.1330 ( —) 37
17 Normalized profiles of vertical spread from
a line source ground release into a stable
atmosphere (Ri = 0.134, Ro = 106) as a
function of distance downwind, denoted in
powers of 10 meters (0 is Gaussian distri-
bution assumed for x = 0) 39
18 Normalized profiles of vertical spread from
a line source ground release into an
unstable atmosphere (Ri = -0.222, Ro = 106)
as a function of distance downwind, denoted
in powers of 10 meters (0 is Gaussian
distribution assumed for x = 0) 40
19 Concentration isopleths downwind of a point
source ground release into a neutral atmos-
phere (Ro =50000, u/f = 5000 m); notation
refers to tenths of maximum value, i.e.,
9 is 0.9cmax, etc. 42
20 Concentration isopleths downwind of a point
source ground release into a stable atmos-
phere (Ri = 1.1, fto = 340000, u/f = 2740 m);
notation as in Figure 19 43
-------
Figure No. Page No,
21 Concentration isopleths downwind of a point
source ground, release into an unstable
atmosphere (Ri = -1.1, Ro = 4750, u/f =
4750m); notation as in Figure 19 45
22 Comparison of the vertical spread as a
function of distance downwind of a line
source for ground and elevated release
into a neutral atmosphere with Ro = 106 46
23 Concentration isopleths downwind of a line
source elevated release into a neutral
atmosphere. Concentration levels are
consistent with the numerical experiment
of Lamb et al.31*: our predictions 47
24 Concentration isopleths downwind of a line
source elevated release into a neutral
atmosphere. Concentration levels are
consistent with the numerical experiment
of Lamb et al.31*: their results 49
25 Predictions for the mean velocity and
shear stress distribution in an equilibrium,
self-similar, neutral shear layer: data as
observed by Wygnanski and Fiedler 53
26 Predictions for the distributions of the
second-order correlations of the normal
velocity fluctuations in an equilibrium,
self-similar, neutral shear layer; data
as observed by Wygnanski and Fiedler1*0 54
27 Maximum species concentration downwind of
a point source release at the center of a
200-m thick shear layer with the wind going
from 6 m/sec at .4 km to 10 m/sec at .6 km 56
28 Vertical and horizontal dispersion as a
function of distance downwind for the same
plume as in Figure 27 57
VI
-------
I. CONCLUSIONS
The A.R.A.P. computer model for predicting the fate of
nonchemically reacting pollutants released in the atmospheric
boundary layer has been extended and its utility demonstrated
by detailed calculations. The program extension to provide
the capability of computing unsteady, two-dimensional turbu-
lent flows opens up the possibility of making a wide variety
of interesting calculations. A number of new comparisons
between model predictions and experimental observations
provides added verification of the model's fundamental
validity.
Specific conclusions from the detailed model calculations
are:
1) Modeling of the temperature fluctuations is given
strong support by the close agreement between model predic-
tions and experimental observations in the wake of a two-
dimensional, slightly heated cylinder.
2) Comparisons between laboratory results and model
predictions of turbulent fluctuations in the limiting case of
free convection when the turbulence is completely produced by
buoyant forces rather than mean velocity shear verify the
model's basic validity.
3) Comparison of the model predictions with laboratory
simulation of diffusion in a free convection, mixed layer is
good, with the model able to predict the maximum concentration
rising from the ground as observed in an instantaneous line
release.
4) Comparison with the Gaussian plume dispersion estimates
given by Pasquill shows close agreement for a line source
released at ground level under neutral atmospheric conditions.
But the calculations indicate that Pasquill's variation with
zo should be interpreted as a variation with Rossby number.
5) The specification of one stability parameter, a
stability class or Richardson number, does not adequately
specify the influence of atmospheric stability on plume
dispersion. Some information on the time history of the
turbulence must be provided, such as is contained in the
mixing layer height.
-------
6) For a point source release in a stable atmosphere, a
strong departure from a Gaussian plume distribution is caused
by variation in mean wind direction with altitude. The
horizontal dispersion due to the wind veering with altitude
appears much greater than that caused by horizontal wind
fluctuations.
7) Comparison with the PSU model predictions for an
incompressible, two-dimensional wake shows little to be gained
by the incorporation of their more complex second-order
closure model in its present form.
-------
II. RECOMMENDATIONS
The process of extending the model capability concurrent
with making detailed calculations with the existing model
should be continued. The current model considers the pollutant
as passive with no feedback on the turbulence field. We
recommend that this restriction be removed by incorporating a
separate momentum equation for the species. This will permit
buoyant plumes and plumes with significant initial velocity
to be calculated properly. This increased capability is
highly desirable because of the relatively high rate of dilu-
tion which often occurs during the rising phase of real plumes.
This report gives the results of a number of detailed
calculations for various atmospheric conditions. We believe
such calculations should be continued, both for the purposes
of further detailed verification by comparison with available
measurements and to exemplify the various possible phenomena
to be expected in unmeasured situations. Such calculations,
using the existing A.R.A.P. programs, should also be used to
provide direct technical support to EPA personnel, including
the provision for utilization of A.R.A.P. programs on the EPA
UNIVAC 1110 by EPA personnel.
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 urban or regional air quality
problems. At the present time, it is not clear whether this
can be most fruitfully accomplished by solving the complete
three-dimensional unsteady problem numerically or by approxi-
mating the problem in some fashion. A standard technique for
reducing the complexity is to assume a succession of steady-
state, three-dimensional, Gaussian plumes. It may be possible
to use our model to parameterize non-Gaussian plume distribu-
tions for incorporation into such a model. Another, perhaps
more promising, possibility is that of using vertical integrals
of various moments of the desired flow properties to reduce the
inherently'three-dimensional, unsteady problem to a two-
dimensional, unsteady problem. We believe it is time to do
some preliminary investigation of these and possibly other
alternative approaches to determine which appears most feasible
for the urban model.
-------
It is also desirable to continue refinements to the
basic turbulence model as better data become available from
experiments, either physical or numerical (such as that being
conducted by Flow Research Inc. of Boston, Mass, for EPA), or
as the desirability for more complicated models for various
terms in the equation is established by other modeling efforts
at A.R.A.P. or by the efforts of others such as Lumley or
Mellor.
-------
III. INTRODUCTION
The development of a program for practical predictions
of the dispersal of pollutants in the atmosphere based on a
second-order closure scheme under development at A.R.A.P. was
first funded by EPA in February 1971. To our knowledge, this
represented the first attempt at applying higher-order turbu-
lence modeling to the pollutant dispersal problem. Previously,
in a paper presented in January 1970* on the application of
second-order closure modeling to the problem of clear air
turbulence, it had been mentioned that there would be no real
difficulty in extending such a model to permit the calculation
of the dispersal of passive atmospheric pollutants in the
atmospheric boundary layer under arbitrary stability conditions,
Through the work at A.R.A.P.2"9 and the work of others10'19,
the use of the explicit Reynolds stress equations to compute
correlations of turbulent fluctuations with modeling techniques
applied to close the system of equations has now become a
viable means of solving problems involving turbulent transport.
This report details the progress made at A.R.A.P. in fiscal
1975 towards utilizing a model based on such a technique to
solve pollutant dispersal problems.
The computer program simulating unsteady one-dimensional
(or steady two-dimensional flow that is parabolic in one
direction) flow has been continually refined since its earliest
version, the Clear Air Turbulence (CAT) program.3 The CAT
version included the invariantly modeled turbulence equations
with known velocity and temperature gradients (functions of z
only). Concurrent with this development was the construction
of the Line Pollutant Dispersal (LPD) program which computed
one-dimensional unsteady pollutant dispersal given all the
necessary turbulent fluctuations and mean flow gradients. A
typical pollutant problem was handled by assuming distributions
for velocity U and temperature 0 and executing CAT to
compute the steady turbulent correlations, obtained by running
the CAT program to large times. These results were then used
as input into the LPD program to predict the dispersal of a
line source.
The first step in the refinement of CAT was Free CAT, in
which the turbulent equations of motion were regrouped in the
solution scheme to include the dynamic behavior of the U and
-------
6 mean equations. The excitement these solutions generated
stimulated the development of a more complete model for atmos-
pheric flows. That model became CORE, a program for solving
the complete set of one-dimensional unsteady equations,
including the coriolis body force. Here the V mean velocity
is included in the solution scheme, as is a pollutant species
C . All of the requisite turbulent correlations are also
present within the computer program. CORE was written to
permit the solution of the turbulence either by their modeled
differential equations or by a refinement of the superequili-
brium limit3 into the quasi-equilibrium equations6 where the
q2 turbulent correlation is the only one solved in its complete
differential form. CORE was also extended to permit solution
to many of the types of laboratory flows for which A.R.A.P.
verifies its turbulent modeling. One of CORE'S improvements
over CAT and Free CAT is the matching of the lower boundary
condition into the law-of-the-wall and Monin-Obukhov relations.
All of these refinements were an effort on our part to permit
one-dimensional solutions in as short a computer time as
possible while retaining the accuracy needed to yield a reli-
ably predictive result.
Concurrent with the refinements in the turbulent program-
ming and the invariant model itself came the further development
of the pollutant programming. Even though the line dispersal
is a valuable step towards predicting and understanding pollu-
tant dispersal, it was realized that a really worthwhile
modeling effort would have to include a point-source or
"spatial" pollutant dispersal (SPD) program. The development
of the initial program was funded by EPA. The numerical
procedure used in the solution of the partial differential
equations was an extension of the same technique used quite
successfully in the one-dimensional program. SPD computes
the steady three-dimensional solution of the modeled turbulent
diffusion equations when the ambient velocity, temperature,
and turbulent field are prescribed or previously calculated
from a program such as CORE. As a part of the current effort,
SPD was programmed on the EPA UNIVAC 1110. The experience
gained in the transfer of this program from our META-4 compu-
ter to the UNIVAC 1110 is discussed in the next section.
In mid-1972, A.R.A.P. began a two-year effort with the
Navy for the development of a submarine wake program. The
rudimentary features of SPD were taken over into a model for
turbulence in the two-dimensional marching framework (WAKE).20'21
At that time, computer availability was such that the turbu-
lence was entirely the quasi-equilibrium form of the invariant
equations. Only the q turbulent energy correlation was
retained in its fully differential form. Mean flow variables
consisted of a density perturbation and three velocities.
Due to the interaction between potential and kinetic energy
-------
leading to horizontal and vertical velocities and internal
gravity waves, it is essential to predict the perturbation
pressure accurately. One of the main results of the Navy
effort was the development of a solution technique to solve
the Poisson equation in two dimensions.
Section IV deals with the development of a two-dimen-
sional marching program incorporating the turbulence features
of CORE, the pollutant dispersal experience of SPD, and the
pressure solution of WAKE. Section V deals with new compar-
isons between model results and experimental observations
which provide additional verification of the basic turbulence
model. The results of detailed diffusion calculations for a
number of different atmospheric conditions are presented in
Section VI. Results for both line and point sources are
compared with dispersion estimates for Gaussian plumes.
The results of comparison calculations for the comple-
mentary turbulence models currently supported by EPA at Flow
Research Inc. and at The Pennsylvania State University are
described in Section VII. Section VIII contains some comments
on the questions related to second-order closure raised by
Drs. P.B. Smith and F. Pasquill in their recent paper, "Some
Views of Modeling Dispersion and Vertical Flux."
Two appendices are included. Appendix A details the one-
dimensional unsteady equations, together with the boundary
conditions appropriate for diffusion in the planetary boundary
layer. This provides a convenient reference for the complete
mathematical specification of the problems solved in Section
VI. Appendix B is a reprint of an AIAA paper, to be published
in the AIAA Journal, giving some comparisons between model
predictions and experimental observations for some variable
density flows.
Part II 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 determination
of the model coefficients, the validation tests, and sample
applications.
-------
IV. PROGRAM DEVELOPMENT
INTRODUCTION
One task of the current EPA contract period was the
development of computer code incorporating the invariantly
modeled turbulence equations for the solution of unsteady two-
dimensional flows. By the incorporation of an additional
advection term, this same code can be used to solve steady
three-dimensional flow problems that are parabolic in one
direction and thus provide one direction for numerical
marching. Our thrust was directed towards a generalization
of the one-dimensional unsteady CORE turbulence program, with
the structure of the spatial pollutant dispersal (SPD) program
as a guide.
TRANSFER OP PROGRAM TO EPA UNIVAC 1110
Since we anticipated the need for access to the EPA
computer as the only substantial way to resolve the future
computational demands of this large program, our first step
in its development was to make the current SPD computer code
operational on the Raleigh facility. In mid-November 197^,
we spent a pleasant and most instructive week in Raleigh.
Primed with many of the necessary manuals giving us the infor-
mation we needed about conversion to the UNIVAC 1110, our task
was accomplished rather painlessly. Special thanks must go to
Bob Jergens and Ginger Smiley for their cooperation and help-
fulness during this time. Since November we have made a number
of production runs using the SPD version on the 1110. The
speed ratio between our META-4 and the UNIVAC machine was
30/1, enabling us to complete a simulation in less than an
hour. Our only displeasure with the set-up has been the
expected one of having the machine down for maintenance quite
often and of having to wait until after the dinner hour to
obtain a communication line. Some of the results presented in
Section VI were obtained on the UNIVAC 1110.
DEVELOPMENT OP THE TWO-DIMENSIONAL UNSTEADY MODEL
The SPD program3 solves for the sp_ecies_distributipn C
and its turbulence correlations cv , cw , c6 , and cc by
-------
solving numerically the modeled partial differential equations
DC
, .
9x2
- cuj HT
Aqcu.
--i (2)
H: - ^ Hr
J J
Dec. = _2—- 3C_ + 9 LA 9cc \ 2bsqcc
Dt ^CUj 3x. Vc 9x qA 9x A
J J \ J
The effort for this year has been directed towards using
the same computer code structure, to be described below, to
include the solution for the other turbulence variables in
the y-z coordinate system, with marching direction of either
x_ or t . The£3e_ variables_are u ,_v , w , 0 , p , q2 , A ,
uv , uw , vw , u0 , v6 , w0 , 92 , vv , and ww . The
modeled equations governing these variables are given in
Part II of this report. They consist of 4 scalar equations,
(2), (3), (26) and (29); 2 vector equations, (1) and (25);
along with the tensor equation (24) representing 6 scalar
equations. For convenient reference, the equations and bound-
ary conditions are detailed in Appendix A for the one-
dimensional planetary boundary layer code. The resultant
computer code is called the Planetary Boundary Layer (PEL)
code.
A great deal of sophisticated coding is required to
transform the differential equations into finite difference
form and place them within an auxiliary code that performs
the many tests and numerical controls needed to maintain the
solution accuracy. For our present purposes, we will attempt
to give the basic information essential to the reproduction
of our work. If the inclination of the reader leans else-
where, we suggest that he now skip to the last few paragraphs
of this section for a discussion of our simulation comparison.
-------
The necessary nonlinear parabolic partial differential
equations are written in finite difference form using a two-
point forward analog for the marching step At = tn - tn
and a three-point centered analog for the two spatial
dimensions. In order to make efficient use of a relatively
modest number of mesh points and to solve accurately, not
only near the lower surface but also out to large vertical z
heights, we determined to continue using our unequally spaced
spatial grid mesh. Thus, Ay and Az need not be the same
within the solution space and, in fact, Az may vary by as
much as two orders of magnitude from the surface to the upper
level geostrophic wind. We permit the program to adjust this
mesh spacing internally as the solution progresses and as
high curvature regions within the flow disappear or form. The
major drawback of the unequally spaced mesh is that it does
not permit a second-order accurate solution scheme.
For two-dimensional flow, our finite difference scheme
must be written on the cross-like mesh shown in Figure 1
where a typical first derivative (say, 3vv/9y) and a typical
second derivative f ^— qA ^p-) at the point (yj,Zi) at time
tn would be formed as follows. We define
Am = y. - y. , A=y.,1-y. A, = A + A
m ''j j-1 p Jj+l Jj t m p
Ar, Am
F = ,-2- p = _JJ_
m Vt P Vt
P = J J •*• r< = J ' -1-
Vt "" Vt
and write these two example derivatives as
9v~v i p / nrr rr^
9y i
Equations (5) and (6) may be extended in a straightfor-
ward manner to yield the entire finite difference formulation
of the equations for the twenty means and correlations. The
first major step is now complete. The next step is to decide
10
-------
J-"
H h
j+l
Figure 1. Definition sketch for finite difference scheme
11
-------
how these equations will be solved. We elect to use the
alternating-direction implicit (ADI) tridiagonal algorithm
developed by Peaceman, Rachford and Douglas.23 To do this,
we divide the marching step At into two equal parts (At/2)
and integrate the equations first in the z direction, then
in the y . During the z integration half-step, z
derivatives in the unknowns are treated implicitly, while y
derivatives are evaluated explicitly with the presently known
solutions. The reverse is true during the y integration
half-step. The implicit algorithm requires the construction
of two vector arrays and the demand for sweeping twice through
the grid mesh.
Once the solutions are known at the next step, tn ,
we test these solutions by determining the maximum variable
change normalized by its maximum variable value. Since we
wish to hold the solution to a given tolerance (2% change,
for example), we then have a rather simple scheme for deter-
mining our marching step size At . Since curvature tests
control our mesh in the y and z directions and the ADI
algorithm integrates the equations, our necessary input to
execute the program becomes consistent variable initial and
boundary conditions. More details on the methods used for
mesh control and step size are given in Reference 9.
An examination of the equations that we are solving
should convince the reader that they are strongly coupled. A
strict use of the ADI scheme would require the inversion of
large matrices twice at every grid point per step. The
computing time required would be prohibitive since we must
maintain enough grid points in each direction to resolve the
solution (approximately 30 points). Thus, it becomes obvious
that we must not only linearize the production terms in an
advantageous manner (we evaluate them at their known values
from the previous step), but we must resolve the coupling
inherent in the mean flow equations. For example, for v ,
we have, from Eq. (i) in Part II,
If
Equation (7) shows that_ v is implicitly coupled to not only
u but also vv and vw (p is either prescribed or obtained
from a solution of the Poisson equation found by taking the
divergence of the momentum equation). However, through our
experience with this effect in_the CORE_ program, it has proven
advantageous to replace the vv and vw derivatives by
rewriting (7) as
12
-------
Dv _ 9p . 9
Dt ~ 9y 3y
v
3v
3z
vv + e I*
3_
3z
vw + e r^
3v
9z
(8)
where e = vvA /q and e = wwA /q
J J £i £i
(9)
and the last two terms in (8) are evaluated explicitly. Then
the second and third terms in (8) may be recast into implicit
terms in v , lending a greater stability to the v equation
since these terms are diffusion-like. The v equation is
then coupled only to u by the coriolis body force. Similar
dynamic couplings due to the_buoyancy_force must exist
between w9 and 62 , and cw and c9 . Linearization of
some terms by evaluation of the coefficients from the pre-
vious step permits the rest of the equations to become
uncoupled from each other. The computer time to solve this
equation set is tolerable — on the order of 4 minutes of
META-4 time per step for a 30 * 30 mesh when the specie is
not present, an extra minute when it is. A solution to
Poisson's equation will add 2 or 3 minutes to the solution
time per step for this size mesh.
NUMERICAL TEST COMPUTATIONS
After an extensive time spent debugging the program and
uncovering several new numerical problems associated with the
addition of the cross-dimension y , we have been able to
test the program on two problems: the steady-state solution
of a simple neutral planetary boundary layer, and the steady-
state solution of flow over a simple ramp change in the sur-
face roughness. The primary reason for choosing these two
test cases is that we were able to run the CORE code, which
marches in the y direction, side by side with PEL, giving
us the steady-state solution for all z , and thus giving us
the solutions that PEL should yield in the limit of large
times. The greatest disadvantage of choosing to run PEL to
steady state was that we had to interpret "steady state" some-
what loosely because of the computer time a typical run
requires. Much of the difficulty may be traced to the convec-
tion terms present in the substantial derivative
D_ _ 3_
Dt " 3t
(10)
13
-------
Since our flow for v is from negative to positive y , one
can see rather crudely that our time step is controlled by
so that the solution does not move so fast past the y grid
points that it cannot be resolved in the time step we are
taking. This is the so-called Courant condition. It
restricts our step size in the marching direction even with
a careful selection of v and Ay .
After imposing this cross-velocity step-size restriction,
we were able to approach the "steady" solution predicted by
the CORE program for the neutral planetary boundary layer run.
The second test was the more stringent — that of a ramp
change in surface roughness. For this case, we held the
geostrophic wind at 10 m/sec and zo = 0.0001 m for y <_ -10 m,
z0 = 0.1 m for y >_ +10 m , and a linear ZQ in between. The
solution was begun everywhere with the steady-state CORE
solution for z0 = 0.0001 m. Readjustment occurred rapidly.
The Courant condition, for Ay = 10 m, required us to restrict
At ~ 2 sec (liberally). After several hundred seconds, we
compared the PEL prediction for turbulent energy q2 with
the "steady" solution from CORE over this same ramp change.
The two curves are given in Figure 2. The agreement is rather
encouraging. The slight differences are probably due to the
relatively larger mesh spacing of the PEL program.
Further work with PEL will include more solution compari-
son but also simulation of flows we cannot in fact compute
other ways. Once the pressure loop is implemented, along with
the differential equation solutions of both the cross-plane
velocities v and w , we will be able to simulate a wide
variety of fundamentally interesting flows.
For this test case, we have also included the computation
of diffusion from a line source at x = 0 . Comparison between
the predicted concentration isopleths for the two programs are
given in Figure 3. Again, the results of the CORE program
should be more accurate since it uses a much finer mesh. The
accuracy of the new PEL appears acceptable, but more tests
will be made on it before it is declared operational.
-------
lOOr
VJl
-200
Figure 2.
1000
1200
1400
1-4 x, meters
Comparison of the turbulent velocity fluctuations in the vicinity of a
change of roughness as computed by the steady-state CORE program ( )
and by the new PEL program ( )
-------
cr\
IUU
80
o>
« 60
E
M
40
20
-
—
/
/ .
•^
C/Q =
.0005 sec/m'
X
/
600 800
x, meters
1200
1400
Comparison of the concentration isopleths for a line source located at
x = 0, z = 0 for the turbulent fields shown in Figure 2
-------
V. NEW VERIFICATION OF THE BASIC TURBULENCE MODEL
The wide variety of flow geometries to which the model
has been successfully applied is detailed in Part II. Here
we want only to describe two new validation tests made during
the last year under this EPA contract. One of these is a
comparison with laboratory measurements of temperature fluctua-
tions in the wake of a heated cylinder. The second is a
comparison with laboratory measurements of both the basic
turbulence field and the resulting dispersion of a passive
species which occurs in free convection bounded by a stable
temperature gradient. The principal results of these two
comparisons were included in a recent AIAA paper24 which is
included herein as Appendix B.
Since the modeling involved in the temperature correla-
tion equations had previously been validated only by compari-
son with atmospheric surface layer data5, it was desirable to
check the results for a flow problem which did not involve
either boundaries or stability influences. Such is the case
far downstream of a slightly heated cylinder. Freymuth and
Uberoi25and, more recently, LaRue and Libby26 have made meas-
urements of the temperature variance and mean temperature
increments obtained in the self-similar region of the wake.
The corresponding self-similar profiles predicted by our model
are compared with the data of LaRue and Libby in Figure 1 of
Appendix B. The normalizing length ZCJQ is fixed by compari-
son with the mean velocity profile (see Section VII for a
comparison of the velocity profiles). The agreement is quite
good, particularly in view of the fact that no new constants
or adjustment coefficients are required (or permitted) in
this comparison. The centerline value of the temperature
variance normalized by the mean temperature increment predicted
by the model is 0.276. This agrees with the value of 0.28
observed by LaRue and Libby rather than the value of 0.21
reported by Freymuth and Uberoi.
Free convection is a flow that is quite different from
the flows used to determine the model coefficients. In this
case, the turbulence is produced by the unstable buoyancy force
without being influenced by any mean shear stress. Willis and
Deardorff27have measured the velocity and temperature
17
-------
fluctuations in an unstable layer bounded above by a stable
temperature gradient and below by a surface with a positive
heat flux. Their experiment called for establishing an
initial stable temperature gradient between two surfaces and
then applying constant heat flux to the region through the
lower surface. With mean velocities restricted to a minimum,
free convection occurs in a mixed layer above the lower
surface. The thickness of this mixed layer increased with
time, but the turbulent correlations were observed to exhibit
self-similarity when normalized by the characteristic velocity
1/3
«
" •
where w0 = surface heat flux
o
z. = depth of the mixed layer
Mo_del predictions_for the normalized plots of ww/w*2 ,
uu/w*2 , 97w*2/(w9)o, and wq2/2w*3 are compared with the
laboratory observations in Figures 3, ^» 5, and 6, respect-
ively, of Appendix B. The comparison of the vertical velo-
city fluctuation is very good. The horizontal velocity
observations show modest peaks near the top and bottom of the
mixed layer that are not predicted by the model, but the
prediction for the center half of the layer is in good
agreement. The temperature variance predicted by the model
and observed in the experiment is the same except near the
top of the mixed layer where the observations show a much
stronger local maximum than is predicted. The inertial
oscillations existing in this region do not appear to be
correctly simulated by the model.
The favorable comparisons of Figures 3, ^, and 5 of
Appendix B are in spite of the fact that the modeled vertical
flux of kinetic energy is significantly different from the
experimental observations, as shown in Figure 6. This is
consistent with our basic assumption that the second-order
correlations calculated from their conservation equations are
given more accurately than the accuracy to which the third-
order correlations are approximated.
Deardorff and Willis28 have also simulated diffusion in
a free convection mixed layer by tracking the release of a
large number of small, neutrally buoyant particles in their
convection tank. Figure 4 shows our model predictions for the
ensemble averaged pollutant concentration distribution a
distance x = 0.44zj_U/w* downwind of the continuous point
source. For this calculation, we have taken U to be a
constant between z = 0 and z^ to correspond to the fact
that the experiment was actually performed as an instantaneous
18
-------
y/h
Figure 4. Concentration isopleths at x = O.^z^U/w* downwind
of a continuous point source release into a free
convective layer as simulated in a laboratory test
by Deardorff and Willis28 ( ) and as predicted by
our model ( ). Numbers on the contours are in
nondimensional units of z^UC/Q .
19
-------
line source experiment. Figure 4 also shows the average
concentration distribution from their five experiments. The
asymmetry about the y = 0 axis is apparently a measure of
the sampling error due to basing their average on a relatively
small number of samples. For both the predictions and the
observations, the concentration is normalized by multiplying
by z£u/Q where Q Is the source strength. The corres-
pondence between the model predictions and their observations
is quite encouraging. Even the lateral spread is in reason-
able agreement.
Figure 5 shows our predictions for the centerline concen-
tration isopleths in the plane y = 0 . Their experimental
contours are shown in Figure 6. In both the model predictions
and the experimental contours, the altitude of maximum concen-
tration rises downwind of the ground release. As pointed out
by Deardorff and Willis,28 this result could not be simulated
by either a Gaussian plume model or an eddy coefficient model.
In our model, this effect could be increased slightly by
increasing d to a value significantly greater than 1. Note
also that the experimental contours are taken from only one
experiment and thus should be expected to be somewhat differ-
ent from the curve which would be obtained for an ensemble
average.
Although some discrepancies between model predictions
and experimental observations have been noted in these new
validation tests, taken as a group they provide strong evidence
to the fundamental correctness of the model.
20
-------
xw*7hu
Figure 5- Centerline concentration isopleths downwind of a continuous point source
release into a free convective layer as predicted by our model; notation
as in Figure 4
-------
r\j
ro
°o
Figure 6.
.4
.8 1.0
xw*/hu
1.2
1.4
1.6
Centerline concentration isopleths downwind of a continuous point source
release into a free convective layer as simulated in a laboratory test by
Deardorff and Willis;28 notation as in Figure 4
-------
VI. DISPERSION ESTIMATES COMPARED WITH GAUSSIAN PLUME MODELS
INTRODUCTION
A popular method for estimating pollutant dispersal is
to assume a steady-state Gaussian plume with the spatial
dispersion of the plume given empirically.29 The purpose of
this chapter is to compare the predictions of our model with
the latest recommendations for Gaussian plumes as given by
Pasquill. 30
AMBIENT TURBULENCE FIELD
The details of the turbulence model used herein are
identical to what is used in Reference 6. All of the turbu-
lent velocity fields are computed solving the full set of
equations (Eqs. 2.1-2.3 and 2.24-2.26 of Part II) with the
one-dimensional unsteady CORE program. These equations and
boundary conditions are repeated here in Appendix A. The only
difference from the quasi-equilibrium solution used in Refer-
ence 6 is in the inclusion of all the diffusion terms, since
we will be dealing here with near-equilibrium turbulent
conditions.
The complete velocity field predicted for a neutral,
steady-state, planetary boundary layer is given in Figure 7.
The governing equations show that the only dimensionless para-
meter required to completely define this flow is some form of
the Rossby number. Figure 7 is shown for Ug/zof = 106 . A
more easily measured velocity than the geostrophic wind in
practical cases may be the velocity at a set height, say 10 m.
A popular theoretical choice for a normalizing velocity is the
surface shear stress velocity u*. Here we will use u = 10$
of the wind velocity at the 10 m height. This completely
arbitrary choice has the advantage of being easier to measure
than either Ug or u* and numerically is the same order of
magnitude as u*.
We will not show all of the turbulence variables for all
of the diffusion calculations to be discussed in this chapter.
But in order to show the flavor of the variation with
23
-------
IV)
.03
.03r
-.002
.002 .004 .006
Figure 7. Velocity and turbulence distribution predicted for a neutral, steady-
state, planetary boundary layer for a Rossby number U /z f = 106
B O
-------
Ro = u/z0f , Figure 8 shows (ww)1' /u for three different
values.
The influence of stability introduces at least one more
parameter — some measure of the stability such as a Richardson
number or Monin-Obukhov length. Usually one must also know
something about the time history of the PEL. Certainly in the
unstable case, some specification of the height of the inver-
sion, i.e., depth of the mixed layer, must be given since, in
this case, no steady-state turbulence velocity field is
possible. The depth of the mixed layer and consequently the
whole turbulence field will continue to evolve in time. But
as shown in Section V and in Appendix B, a quasi-steady turbu-
lence can occur for a given inversion height z.j_ . In our
program, Zj_ is not a rigid lid but is defined as the alti-
tude at which the maximum values of F2" occur as the temper-
ature inversion is approached.
Here we have chosen to measure the stability by the
Richardson number based on the potential temperature and wind
velocity taken 10 m above the surface:
e0u2
Again, we have sacrificed theoretical neatness in order to
choose a parameter which may be easily measured. The
Richardson number defined this way is roughly equal to the
gradient Richardson number at the 10 m height. When the
logarithmic profiles for u and 9 are appropriate, the
relationship may be written exactly as
Rl =
0
o
(3u/9z)2 z An z/z
100 (13)
0 z=10m
The correspondence between Ri and Pasquill's stability
parameter P is shown in Figure 9« For convenience, we have
also shown the correspondence between P and Pasquill's
former stability classes.
The implication in the discussion of Figure 6.13 in
Reference 30 is that the influence of surface heat flux on
turbulent fluctuations is uncoupled from the influence of
zo . This does not appear to be the case. The curve in
Figure 9 is drawn for zo = 10 cm. If the definition of P
remained fixed, as defined in Figure 6.13 of Reference 30,
-------
rf.
/v
u
Ro= 3060
.4
8 1.2 1.6
(ww) /u
2.0 2.4
Figure 8. Ratio of the vertical velocity variance to 10% of
the .wind at 10 ra altitude for three different
values of Rossby number
26
-------
O
Q
z0= 10cm
Zj =1
z0 = I cm
Z =3Km
-1.0 -.8 -.6 -.4 -.2
x\
Ri
T
t
D
C
B
\
0 .2
Figure 9- Correspondence between Ri and Pasquill's
stability parameter P or stability classes
A-F
27
-------
for other values of z0 , then the points in Figure 9
indicate a nonuniqueness between P and Ri . However, it
seems more_desirable to have P defined so that at a fixed
altitude ww/U2 is determined uniquely. This would require
the definition sketch of P (Fig. 6.13, Ref. 30) to shift
with z and should make the relationship between P and
R"i contain less scatter.
^Figure 10 shows (ww)1//2/u for four different values
of Ri and two different values of the dimensionless inver-
sion height Zj_f/u . The inversion height represents addi-
tional information supplementary to the surface meteorological
conditions and must be known if estimates of dispersion very
far from the source are to be made. It may typically vary
diurnally from something less than 100 m near sunrise to a
few kilometers near sunset.32 Note tfcat for z/Zj_ £ 0.1, the
two curves for essentially the same Ri, but different Zj_ 's ,
show a significantly different value of (ww)1'2 .
SUMMARY OF DIFFUSION MODEL
As developed in References 3 and 6, the diffusion of a
passive species calls for the solution of the species contin-
uity equation, the species turbulent flux equation, and the
species-temperature correlation equation (Eqs. (l)-(lJ) from
Section IV) with the velocity turbulence field known from a
prior decoupled calculation. In these equations,
qp - qt(Ap/At
where the subscript t represents the value of q or A
obtained from the prior turbulence calculation.
The greatest difficulty in using Eqs. (1) through (4)
is in properly setting the pollutant velocity scale q and
length Ap . When the size of the plume is large in compari-
son to the most energetic turbulent eddies, qp and Ap may
be taken as equal to the same scales used in computing the
turbulence. However, when the plume scale is much less than
the most energetic eddies, then both the velocity and length
scales of the plume may be significantly smaller than the
companion scales for the turbulent velocity field. Our
approach takes the plume scale proportional to the spread of
the species plume
Ap = do (15)
subject to the constraint
Ap < At (16)
28
-------
.8r
.7
Ri=-.222 (z;f/u=.45l)
(ww)l/2/u
Figure 10.
Ratio of the vertical velgcity variance to u
for different values of Ri and mixed-layer
depths for Ro = 106
29
-------
That is, the mixing scale cannot be larger than the ambient
turbulent scale. The proportionality factor d is a model
input coefficient. It is the only new coefficient introduced
in the species equation that has not previously been deter-
mined and used in computing the velocity and temperature
distributions.
A discussion of the relation of our model equations to
K theory has been given in the chapter on Super-Equilibrium
Theory in Reference 3. The validity of K theory requires that
the convective terms (those on the left-hand side of Eq. (2))
and the diffusion term (the term with the second-order
derivative) be negligible in comparison with the remaining
terms. For a one-dimensional problem with neutral stability,
this reduces to
—- wwA 3 C
wo = - -~
i.e. ,
This should be reasonably valid near the surface under neutral
conditions where the expression reduces to Kz = (0.52)u*z .
We wish to emphasize that this reduction is valid only under
these very special conditions. At higher altitudes and/or for
non-neutral conditions, other terms in Eq. (2) are important.
LINE SOURCE RELEASE AT GROUND LEVEL
We first consider cases where sensitivity to the coeffi-
cient d is very weak. In close proximity to the ground,
the scale of the turbulence is proportional to the distance
from the ground so that Eq. (16) rather than (15) will
limit A .
Figure 11 is a plot of the vertical variance a, for
the two-dimensional problem of a line source perpendicular to
the surface wind in a neutral atmosphere with Rossby number
based on the geostrophic wind, Ro = Ug/fzo = 106 . The
difference between az for d = 1 ana an extreme value of
d = 1000 is slight. Also shown for comparison is Pasquill's
latest recommendation for use with Gaussian plume models
obtained from a combination of eddy diffusivity calculations
and empirical data (Ref. 30, Fig. 6.13b).
Figure 12 compares the vertical distribution of species
concentration at several downwind distances with the Gaussian
30
-------
1000
100
m
10
.1
Pasquill prediction
x, km
10
100
Figure 11.
The vertical spread az as a function of distance
line source ground release in a neutral atmosphere
z.
o
= 10 cm, f = ID"* sec'1 )
downwind of a
= 10 m/sec,
-------
c/c
max
Figure 12,
Comparison of species distribution with a Gaussian
distribution for the plume of Figure 11. Normal-
ized profiles are shown for several distances down-
wind of the release point. The numbers give the
exponent in tens of meters;
distribution assumed for x
x = 100, etc.
32
i.e., 0 is the Gaussian
= 0, 2 is the profile at
-------
distribution. In spite of the shear existing in the wind
profile, the normalized results do not depart strongly from
Gaussian. However, Figure 13 shows that this should not be
interpreted as verification of the commonly used Gaussian
plume approximation for determination of the surface concen-
tration, because this approximation ordinarily assumes a
uniform wind. If the 10-meter value of the wind is used in
the Gaussian plume model, Figure 13 shows that it would predict
a lower value of surface concentration close to the source and
a higher value at a distance, even when the empirical az
variation is used.
Pasquill's Figure 6.13d30 shows the growth of az for
various values of zo but with no indication of variation with
geostrophic wind or coriolis parameter. Our modeled turbulence
equations indicate that the proper dimensionless parameter for
this neutral case is the Rossby number. A normalized curve for
a variation of two orders of magnitude in the Rossby number is
shown in Figure 14. These curves should not be used for
estimating az quite close to the source since we arbitrarily
start with az = 1 m at x = 0 . The dependence on Rossby
number is slightly stronger than that implied by Pasquill's
Figure 6.13.
The influence of stability on vertical spread is much
stronger than the influence of Rossby number. In Figure 15
the stability is measured by the^Richardson number defined in
Eq. (1). On the unstable side (Ri < 0), az is capped by the
presence of the temperature inversion layer at the top of the
unstable layer of thickness z^ . This yields an upper bound
on oz of az * 0.6z^ .
An unexpected feature of Figure 15'is that normalization
of the lengths by u/f scales the curves such that they tend
to converge at low x . If the lengths were normalized by
Ug/f instead, there would be more spread between the unstable
and stable cases. Exact comparison between the influence of
stability, as predicted by Figure 15 and that predicted by
Pasquill, is difficult because no normalization is given for
his lengths. The implication in his curves is that for
specified zo and stability class, there is no dependence of
spread rate on wind velocity, as mentioned in the discussion
of Figure 9. If we interpreted his plots as made for a geo-
strophic wind of 10 m/sec and f = 10"''sec'1 so that
Ug/f = 100 km, as was done in Figure 11, then a direct compar-
ison can be made. This yields the comparison shown in Figure
16. On this basis, Figure 16 shows a stronger effect of
stability on az on the unstable side and a weaker effect on
the stable side than that implied by Pasquill.
33
-------
<=>
M
O
E
10
.1
.08
CO
-tr
.04
02
.01
Figure 13.
,05
.1
.5
10
r /^
xf/u
Normalized surface species concentration as a function of downwind
distance ( ) compared with the Gaussian solution ( ) for a wind
assumed constant at its value at a height of 10 meters
-------
.1
OJ
VJ1
.01
.001
Ro=3060
Ro = l05
Ro = 49,000
Ro=l06
Ro = 6l2,000
Ro=l07
.1
10
xf/u
Figure
Comparison of vertical spread in a neutral atmosphere from a line source
ground release for changes in the Rossby number based on surface
roughness z0
-------
uo
cr\
fk!
s\
U
.01
.001
.01
Ri = -.222
Ro = 67,000
xs.
Ri = 0
Ro = 50,000
Ri = .l34
Ro = 33,000
.1
xf/u
10
Figure 15•
Comparison of vertical spread from a line source ground release for
different values of Ri for Ro = 106 (z^f/u = 0.451 for the
unstable case)
-------
OJ
-J
orz,m
Figure 16.
x,m
'Repeat of Figure 15 in unnormalized form for uE/f = 100 km for
comparison with Pasquill's Figure 6.1330 ( ) 6
-------
Neglecting any influence of ug or f on the variation
of az with ZQ in Pasquill's Figure 6.13 is probably not
too serious for practical calculations, because this will lead
to no more than an order of magnitude uncertainty in Ro ,
which translates into less than a factor of two in the estimate
of az . However, neglecting any influence of the mixed layer
depth on az estimates in unstable atmospheres can potentially
lead to more error. The influence of differing mixed layer
heights at similar unstable values of P can lead to signif-
icantly different values of turbulence even relatively near
the ground, and thus to significantly different dispersion
values. This will be the case when
w
When this condition is met, diffusion above z/z^ 2 0.1 can-
not be determined by stability class alone.
The concentration profiles ^re compared with the normal-
ized Gaussian distribution for Ri ? 0 in Figures 17 and 18.
The distributions for stable conditions remain. similar to the
Gaussian but, in unstable conditions there are strong
departures. Note that at x = 10 km, the maximum concentra-
tion occurs significantly above ground level. After the
vertical spread is effectively capped by the inversion, the
concentration asymptotically approaches a uniform distribution
between the ground and the inversion, with a sharp drop-off
at the bottom of the inversion layer. When the Gaussian plume
model is modified by the method of images to account for the
presence of a lid in air quality modeling, it is possible to
approach the uniform distribution between the ground and the
inversion. However, it will not permit the maximum concentra-
tion to rise above ground level. The ratio of the ground
concentration to that which would be predicted by a simple
Gaussian plume transported by a constant wind departs remark-
ably little for Ri t 0 from that shown in Figure 13 for the
neutral case. The plume for the stable case shows a concentra-
tion ratio approximately 10% higher at xf/u = 10 , while the
unstable is approximately 10% lower. Apparently the wind shear
is the dominant variable in determining this concentration
ratio, and u fairly effectively measures this influence.
POINT SOURCE GROUND RELEASE
To this point we have shown results that assumed no varia-
tion in the y direction. The more interesting problem in
practice is usually the three-dimensional plume involving
lateral as well as vertical dispersion. Any theory attempting
38
-------
.6
.8
C/C
max
Figure 17.
Normalized profiles of vertical spread from a line
source ground release into a stable atmosphere
(fii = 0.13^, Ro = 106) as a function of distance
downwind, denoted in powers of 10 meters (0
Gaussian distribution assumed for x = 0)
is
39
-------
c/c
max
Figure 18.
Normalized profiles of vertical spread from a line
source ground release into an unstable atmosphere
(Ri = -0.222, Ro = 106) as a function of distance
downwind, denoted in powers of 10 meters (0 is
Gaussian distribution assumed for x = 0)
-------
to compute three-dimensional dispersion must first be able
to predict any valid differences between vertical and lateral
velocity variances. In the neutral boundary layer, our model
predicts the two variances to be equal. Most investigators
believe the lateral to be larger, but there is considerable
uncertainty. Panofsky32 quotes ratios anywhere from 1.2 to
2.4. Our basic turbulence model probably needs an additional
term representing the influence of mean strain on the model-
ing of the terms involving pressure fluctuations, as has been
suggested by several investigators (for example, Lumley and
Khajeh-Nouri,1!* Naot33) . However, with the large uncertainty
in the empirical data, it is difficult to make any such
adjustments with confidence. Even with a model which correct-
ly predicts the lateral fluctuations in an ideal homogeneous
boundary layer, inhomogeneities in terrain appear to have such
a large influence on the low frequency lateral velocity
fluctuations that a confident prediction in a real situation
will be difficult. Notwithstanding these difficulties, we do
believe it is enlightening to proceed with calculations that
show the influence of stability on the three-dimensional
plume.
In this three-dimensional case, there are two different
values of a : ay and az . Therefore, Eq. (15) must be
generalized to determine the plume. For present purposes,
the A 's appearing in the diffusion terms of Eqs. (2) through
(4) are specified by either day or daz , depending on the
direction of the diffusion, while the A 's appearing in
the tendency-towards-isotropy terms are given by an average
value: ~
2a a
A . a t .v V (19)
(oy + °z'
In any case, the constraint in Eq. 0-6) is still applied and,
therefore, there is very little sensitivity of the results to
the particular form of Eq.(19) in the problem of ground
release.
Concentration isopleths downwind of a point source at
ground level in a neutral atmosphere are shown in Figure 19.
The notable distortion from a Gaussian distribution apparent
in Figure 19 appears to be principally caused by the Ekraan
spiral of the principal wind deviation with altitude.
The influence of the Ekman spiral is much more pronounced
for the case of a point release at ground level in a stable
atmosphere, as shown in Figure 20. This is primarily due to
the fact that the angle between the surface wind and the upper
level geostrophic wind is 30° more for this stable case
41
-------
3
Z/
-------
-2
Figure 20.
Z/o-,
\ I 7 "•7'
,-3 Q
(a)
x = I km
cmax= 1.05X10 --§-
o~y s 28.5m
= 8.54 m
-2
-2 -1
0 1
y/o-y
2
r (b)
x = 50 km
o-y = 5480 m
-------
(Ri = 1.1) than it is for the neutral case. The Ekman layer
is also much shallower. The maximum ground level concentra-
tion occurs immediately downwind of release, even though the
plume is distorted asymmetrically by the Ekman spiral. The
uncertainty present in the lateral velocity fluctuations is
overshadowed in these plots by the enhanced lateral spreading
caused by the wind veering with altitude. Therefore, we
believe Figure 20 is a valid prediction although we have not
yet uncovered reliable data to compare with its predictions.
Figure 21 is a complementary plot of a ground release
into an unstable atmosphere (Ri = -1.1) with an inversion
layer occurring at 3.8 km (z^ = 0.355u*/f). For this level
of instability, there is very little effect of the wind
veering with altitude. Note that the upper level inversion
behaves very much like a reflecting boundary (similar to the
ground boundary) in the distribution at 10 km. This is a
result of the model, since we have not specified any such
boundary condition. As seen in Appendix A, the same boundary
condition on the species is specified for all stability condi-
tions as z ->• °° . Further downwind of the release, the
concentration will continue to be quite homogeneous between
the ground and the inversion layer, but with an ever-increas-
ing lateral spread.
Direct comparisons between vertical spread rates for
two- and three-dimensional plume calculations show relatively
little difference. Estimating az on the basis of two-
dimensional calculations appears valid even when the plume is
three-dimensional.
ELEVATED LINE RELEASE
Results for dispersal downwind of a continuous line
source at an elevation which is a significant fraction of the
neutral planetary boundary layer are shown in Figures 22 and
23. Figure 22 shows the influence of the release height on
the vertical spread rate. The vertical velocity fluctuation
is less at altitude than near the ground, and the mean wind
is larger, but these effects are partly, counterbalanced by
the fact that the turbulent scale is significantly larger.
The net effect is a^reduction in az of approximately a
factor of 2 at xf/u = 10 . This difference would disappear1
asymptotically as az becomes greater than the release
height zr .
In Figure 23, the concentration isopleths have been
normalized to make them directly comparable to the results
of a numerical experiment performed by Lamb et al.,31*i.e.,
-------
r (o)
= 5.33X10 --
-2 - 1 0
1 2
y/o-y
3
r(b
)
x = IO km
max
= 5.01X10 •%-
u
-2
o-y= 1280m,
-------
.001
xf/u
Figure 22.
Comparison of the vertical spread as a function of
distance downwind of a line source for ground and
elevated release into a neutral atmosphere with
Ro = 106
-------
.6
u
.4
f_
*
.2
0
0
8
I0
14
Figure 23. Concentration isopleths downwind of a line source elevated release into
a neutral atmosphere. Concentration levels are consistent with the
numerical experiment of Lamb et al.
3"*
our predictions
16
-------
CUru*/fQ where Ur is the mean velocity at the release
height. In their numerical experiment, they released ideal
numerical particles into the three-dimensional, unsteady
velocity field computed by Deardorff35 for a neutral atmos-
phere. Their results are presented in Figure 24. The
greatest discrepancies between Figure 23 (our prediction) and
Figure 24 occur close to the source where the plume calculated
by Lamb et al.31* shows a wider spread than is predicted by our
model. Because of the time step involved in Deardorff's
original model,35 the smallest resolvable normalized Ax is
Ax = 0.2uVff so that details of the concentration distribu-
tion in the immediate neighborhood of the source are not
expected to be valid. This leaves the unanswered question of
how accurate their concentration distribution is at their
first downstream point. Since the problem is an initial value
problem, it would seem that any uncertainties in the early
distribution would remain for some finite time in the
calculation.
An interesting phenomenon occurs near the source in our
model calculations. The plume appears to split and then
reform, with the centerline concentration value at a position
x < 2 km downwind of the source being much less than the
value above or below it. This is a result of the wave nature
of Eqs. (1) and (2) when the plume scale is much less than the
turbulent scale. Physically, this is a manifestation of the
meandering which is often evident when viewing a small plume
in strong turbulence. Since our distributions represent an
ensemble average of such plumes, the meandering instantaneous
plume spends more time at a distance from the downwind center-
line of the averaged plume than it does on the centerline.
This meandering phase has been previously discussed.6
i
It is interesting to note that the computer output for
the experiment of Lamb et al. (supplied to us by Lamb) shows
a value of concentration on the centerline of the plume at
xf/u* = 0.2 which is only approximately 105? of the value at
the z station above or below the centerline. However, for
the reasons stated earlier, this cannot be used to quantita-
tively verify the meander problem.
It is in determining this early meander phase of the
plume that the coefficient d in Eq. (14) has its largest
influence. The sensitivity disappears when the plume size
reaches the scale of the ambient turbulence. Better data or
numerico-empirical analyses able to determine the early-time
behavior of a plume will be required for further adjustment
of this parameter.
^Private communication from R.G. Lamb
48
-------
.8.-
.6
zf/u*
.4
.2
6 8
xf/u*
10
12
14
16
Figure 24. Concentration isopleths downwind of a line source elevated release into a
neutral atmosphere. Concentration levels are consistent with the
numerical experiment of Lamb et al. 3I* : their results
-------
Model results for elevated point source plumes, for both
stable and unstable conditions, have been presented in Refer-
ence 6. That report also contains model predictions of point
source plume dispersion from different stack heights for
unstable conditions which show the dependence (close to a -2
power law) of maximum ground concentration on effective stack
height.
50
-------
VII. COMPARISON WITH TURBULENCE MODELS OP PSU AND PRI
The Environmental Protection Agency is currently funding
two other research groups on complementary turbulence model-
ing studies. The general nature of second-order closure and
an attempt to develop a completely rational model are being
pursued by a group at The Pennsylvania State University under
the direction of John Lumley. Calculations which will hope-
fully provide validation information for closure models are
being attempted by Flow Research Inc. under the direction of
Steven Orszag. The purpose of this chapter is to show the
interrelationship between these efforts and that at A.R.A.P.
and to provide some comparison calculations.
As outlined in Reference 14, Lumley and his colleagues
are attempting to develop models for second-order closure as
rigorously as possible. The inherent generality of this
approach leads to a large number of coefficients which must
be determined. In fact, as Lumley states in Reference 36,
a "dizzying array" of possible terms appears. In contrast,
our approach has been to use a relatively simple closure
model and only add complexity as comparisons between model
predictions and experimental data prove it desirable to do so.
As a result, the PSU approach and the A.R.A.P. approach are
quite complementary. PSU has explored a greater variety of
closure assumptions, while A.R.A.P. has explored a wider
variety of flows with our simple closure assumptions.
One flow which Lumley has calculated is that of a two-
dimensional incompressible wake.11* Comparison between his
model prediction for the self-similar region of this flow and
that predicted by the A.R.A.P. model is presented in Figure
4.3 of Part II, along with Townsend's37wake measurements.
Both model predictions agree with the mean velocity measure-
ments and predict a maximum shear stress approximately 20%
less than Townsend's values. Both models also predict longi-
tudinal turbulence intensity that is 20% to 30% higher than
Townsend's values. However, Lumley's model, which contains
many more terms and coefficients (some evaluated for this
particular flow), does show marginally better agreement of
the lateral and transverse velocity fluctuations with the
observations of Townsend. It should also be noted that there
51
-------
is some uncertainty in the data since Thomas38 reports measure-
ments of the longitudinal rms intensity which are about one-
third higher than Townsend's reported values.
In Figure 1 of Appendix B we compare our predictions for
the temperature fluctuations in this two-dimensional wake
flow with the observations of LaRue and Libby.26 The agreement
is very good. This is a calculation Lumley and his colleagues
indicate they expect to make in the future.
Although we do expect that Lumley's approach will even-
tually lead to model terms that will be an improvement over
our current simple model, particularly for the pressure-strain
correlation, the preceding comparisons indicate that at present
his complex model has no significant advantage over ours.
Second-order closure modeling requires reliable empirical
information for the numerical evaluation of the modeled terms.
As may be noted by the many discrepancies between data
observed by different investigators, good reliable data are
extremely difficult to obtain. An alternative to physical
experiments is numerical experiments. Exact numerical calcu-
lations of the Navier-Stokes equations at realistic Reynolds
numbers for atmospheric flows are beyond the state of the art
of computer technology for the foreseeable future.39 However,
since turbulence in high Reynolds number flow tends to become
independent of Reynolds number, it does appear possible to
make detailed computations at relatively moderate Reynolds
numbers to simulate atmospheric flow and diffusion. Such
computations are being attempted at "Flow Research Inc.
In consultation with Dr. Arthur Bass of FRI, we defined
a problem for mutual computations and comparison between
A.R.A.P. and FRI. The problem agreed upon is that of the
turbulent field and pollutant dispersal in a simulated, neutral
atmospheric shear layer above the influence of the ground. We
specified the shear layer to be 200 m thick with AU = 4 m/sec
and centered at 0.5 km when the wind is 8 m/sec. A point
source of pollutant was to be released into the turbulence
field generated by this shear layer. Such a neutral shear
layer will, of course, spread in time, but for purposes of
this comparison calculation, we have chosen to hold the velo-
city field as frozen at the equilibrium turbulent conditions
occurring for this AU and Az while the diffusion calcula-
tion is made.
Figures 25 and 26 show the A.R.A.P. predictions for the
mean velocity distribution and the mean values of the second-
order velocity correlations. To date, no results have been
received from FRI for comparison. We have included laboratory
shear layer data as reported by Wygnanski and Fiedler1*0 although
52
-------
A.R.A.P. prediction
Ul
UJ
1.2
1.0
.8
.6
.4
.2
o A Wygnanski & Fiedler
-2.4 -2 -1.6 -1.2 -.8 -.4
Z25 Z75
Figure 25. Predictions for the mean velocity and shear stress distribution
in an equilibrium, self-similar, neutral shear layer; data as
observed by Wygnanski and Fiedler1*0
-------
.2
.16
.12
.08
.04
-°2.4
A.R.A.P. prediction
o A D Wygnonski 81 Fiedler
-2 -1.6 -1.2 -.8 -.4
.4 .8
1.2 1.6
Figure 26.
Z25~Z75
Predictions for the distributions of the second-order correlations
of the normal velocity fluctuations in an equilibrium, self-similar,
neutral shear layer; data as observed by Wygnanski and Fiedler w
-------
there is considerable uncertainty surrounding the data as
reported in Reference ^11.
Two computations of the average distributions of the
diffusing species were made for values of d = 1 and
d = 103. Since d represents a model coefficient which
should be determined empirically, good numerical results by
PRI would provide us with the opportunity for fixing this
coefficient. Figure 27 shows a plot of maximum species con-
centration downwind of the source for the two different
values of d . The vertical and horizontal dispersion of the
species are shown in Figure 28. For ease of computation, we
started with an initial Gaussian distribution with a = 10 m
at x = 0 . At this time, it is unclear as to what the
quality of the FRI diffusion computations will be. High
quality results would permit us to select a more definitive
value of d than the value of 1 used throughout most of this
report. Otherwise, we must await high quality experimental
results before completing this refinement.
55
-------
Figure 27.
x, meters
Maximum species concentration downwind of a
point source released at the center of a 200-m
thick shear layer with the wind going from
6 m/sec at .4 km to 10 m/sec at .6 km
5.6
-------
V)
i_
0 2
o IOZ
Nl
b
10
d=!03
d=l
10'
10'
x, meters
Figure 28.
Vertical and horizontal dispersion as a function of distance downwind
for the same plume as Figure 27
-------
VIII. COMMENTS ON "SOME VIEWS OP MODELING DISPERSION
AND VERTICAL FLUX" BY PASQUILL AND SMITH
In a note published in Reference 22, Pasquill and Smith
raised objections to the use of second-order conservation
equations in turbulent modeling. They expressed belief that
closure assumptions would have to be valid to an unobtainably
high accuracy. After an exchange of letters, this objection
has been partially retracted in Pasquill 's latest report, 1*2
but he still expresses some reservations about the approach.
Their objection is based on the observation that for
some flows, particularly ground-level evaporation during a
typical diurnal cycle, the time rate of change of species
turbulent flux will be much smaller than the flux production
term; i.e., in Eq. (2), the terms on the left-hand side may
be orders of magnitude smaller than the largest terms on the
right. Pasquill1*2 gives an estimate of
(20)
at_ 10 m altitude for the time between sunrise and noon with
ww £ 2 m/sec. Based on this, their original argument was
that the modeled terms which provide the balance on the right-
hand side of Eq. (2) must therefore be modeled to an accuracy
of 1 part in 101* . In the latest report, Pasquill w does
admit that when the flux equation is modeled in the form given
iri Eq. (2), the explicit dependence of the modeled_ terms on
we provides stability to the equation so that we can be
evaluated even in the stationary case.
However, Pasquill still raises the question of whether
adequate "settling down" time will be possible to obtain an
adequately accurate value of we in typical cases. He pro-
ceeds to consider a spatially varying case where A/q is
order 1 minute, while the ratio of the characteristic dimen-
sion of an area source to the streaming velocity may be less
than an hour, and suggests that in such a case the accuracy
required in the modeling may be difficult to obtain due to
inadequate "settling time." A check on the ratio of the
58
-------
left-hand side of Eq. (2) to the last term in Eq. (2) shows
that when this ratio R ceases to bejnuch less than 1 is
exactly when a dynamic equation for we is required.
- u An H JL H. (21}
~ u 3x /Aq A ~ Aq x U1;
That is, whenever the characteristic time scale of the turbu-
lent mixing becomes significant in comparison with a typical
transit time over an area source of pollution, then an
equilibrium balance of the terms on the right-hand side of
Eq. (2) is not adequate, and the dynamic balance including
the_JLeft-hand side must be used to obtain a valid estimate
of we .
By his estimates of the two different conditions,
Pasquill has shown_that a reliable diffusion model must be
able to compute we accurately whether the_t_ime rate of
change of we (or the spatial convection of we) is significant
or not. We believe Eq. (2) has this capability. Whenever
R « 1 , then a computation using Eq. (2) will rapidly approach
an equilibrium between the terms on__the right-hand side and
yield_an appropriate estimate of we even if the initial value
of we used for the computation was grossly in error. On the
other hand, when R XX 1 , the dynamic balance shows the
influence of the initial conditions of we must be retained,
and a simple equilibrium balance between the right-hand side
terms would lead to errors. Stated another way, when R « 1,
superequilibrium diffusion theory which reduces to a type of
K theory3 should be approximately valid, but when R is not
much less than 1, such theories are not valid. Although
R « 1 typically in the homogeneous atmospheric surface layer,
this is not the case in spatially inhomogeneous cases and will
not always be the case at altitude in the typical homogeneous
diurnal boundary layer.
59
-------
IX. REFERENCES
1. Donaldson, C.duP., R.D. Sullivan, and H. Rosenbaum.
"A Theoretical Study of the Generation of Atmospheric
Clear Air Turbulence." AIAA Journal. 10_: 162-170,
February 1972.
2. Donaldson, C.duP. "Calculation of Turbulent Shear Flows
for Atmospheric and Vortex Motions." AIAA Journal.
10_:4-12, January 1972.
3. Donaldson, C.duP. "Atmospheric Turbulence and the
Dispersal of Atmospheric Pollutants." In: AMS Workshop
on Miorometeorology, Haugen, D.A. (ed.). Boston,
Science Press, 1973. pp. 313-390.
4. Lewellen, W.S., M.E. Teske, and C.duP. Donaldson.
"Application of Turbulent Model Equations to Axisymmetric
Wakes." AIAA Journal. 12_: 620-625, May 1974.
5. Lewellen, W.S. and M.E. Teske. "Prediction of the Monin-
Obukhov Similarity Functions from an Invariant Model of
Turbulence." J. Atmos. Sciences. 3_0.:1340-1345, 1973.
6. 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." Aeronautical
Research Associates of Princeton, Inc. Report No. EPA-
650/4-74-035- Environmental Protection Agency, Raleigh,
NC. September 1974.
7. Donaldson, C.duP. "The Relationship Between Eddy Trans-
port and Second-Order Closure Models for Stratified Media
and for Vortices." In: Free Turbulent Shear Flows.
NASA SP-321, Volume 1, 1973. pp. 233-255.
8. Varma, A.K., R.A. Beddini, R.D. Sullivan, and C.duP.
Donaldson. "Application of an Invariant Second-Order
Closure Model to Compressible Turbulent Shear Layers."
(Presented at AIAA 7th Fluid and Plasma Dynamics Confer-
ence, Palo Alto, June 17-19, 1974.) AIAA Paper No. 74-
592.
9. Sullivan, R.D. "A Program to Compute the Behavior of a
Three-Dimensional Turbulent Vortex." Aeronautical
Research Associates of Princeton, Inc. Report No. ARL-
60
-------
TR-74-0009. Aerospace Research Labs., Wright-Patterson
AFB, Ohio. 1973.
10. Harlow, F.H. (ed.). "Turbulence Transport Modeling."
AIAA Selected Reprint Series, Volume XIV, 1973.
11. Launder, B.E. and D.B. Spalding. Mathematical Models
of Turbulence. New York, Academic Press, 1972.
12. Mellor, G.L. "Analytic Prediction of the Properties of
Stratified Planetary Surface Layers." J. Atmos. Sciences.
3_0.:106l-1069. 1973.
13. Launder, B.E. "On the Effects of a Gravitational Field
on the Turbulent Transport of Heat and Momentum."
J. Fluid Mechanics. 6j_ (Part 3):569-58l, 1975.
14. Lumley, J.L. and B. Khajeh-Nouri. "Computational Model-
ing of Turbulent Transport." Advances in Geophysics,
Volume ISA. New York, Academic Press, 1974.
15. Zeman, D. and H. Tennekes. "A Self-Contained Model for
the Pressure Terms in the Turbulent Stress Equations of
the Atmospheric Boundary Layer." J. Atmos. Sciences.
32:1808-1813, 1975.
16. Shir, C.C. "A Preliminary Numerical Study of Atmospheric
Turbulent Flows in the Idealized Planetary Boundary Layer."
J. Atmos. Sciences. 3_p_: 1327-1339. 1973.
17. Mellor G.L. and T. Yamada. "A Hierarchy of Turbulence
Closure Models for Planetary Boundary Layers." J. Atmos.
Sciences. 31:1791-l8o6. 1974.
18. Wyngaard, J.C., O.R. Cote, and K.S. Rao. "Modeling the
Atmospheric Boundary Layer." Advances in Geophysics,
Volume 18A. New York, Academic Press, 1974.
19. Wyngaard, J.C. and O.R. Cote. "The Evolution of a
Convective Planetary Boundary Layer - A Higher-Order
Closure Model Study." Boundary Layer Meteorology.
7_:289. 1974.
20. Lewellen, W.S., M.E. Teske, and C.duP. Donaldson.
"Second-Order Turbulent Modeling Applied to Momentumless
Wakes in Stratified Fluids." Aeronautical Research
Associates of Princeton, Inc. Report No. 206. 1973.
21. Lewellen, W.S., M.E. Teske, and C.duP. Donaldson.
"Turbulent Wakes in a Stratified Fluid." Aeronautical
Research Associates of Princeton, Inc. Report No. 226.
1974.
22. Pasquill, F. and F.B. Smith. "Some Views on Modeling
Dispersion and Vertical Flux." In: Proceedings of the
5th Meeting of the Expert Panel on Air Pollution Modeling.
Roskilde, Denmark, Danish Atomic Energy Commission, 1974.
61
-------
23. Roache, P.J. Computational Fluid Dynamics. Albuquerque,
Hermosa Publishers, 1972. pp. 91-95-
24. Lewellen, W.S., M.E. Teske, and C.duP. Donaldson.
"Variable Density Flows Computed by a Second-Order
Closure Description of Turbulence." AIAA Journal (to be
published).
25. Freymuth, P. and M.S. Uberoi. "Structure of Temperature
Fluctuations in the Turbulent Wake Behind a Heated
Cylinder." Physics of Fluids. l4_:2574-2580, December
1971.
26. LaRue, J.C. and P.A. Libby. "Temperature Fluctuations
in the Plane Turbulent Wake." Physics of Fluids.
17.: 1956-1975, November 1974.
27. Willis, G.E. and J.W. Deardorff. "A Laboratory Model of
the Unstable Planetary Boundary Layer." J. Atmos.
Sciences. 31:1297-1307, 1974.
28. Deardorff, J.W. and G.E. Willis. "Physical Modeling of
Diffusion in the Mixed Layer." In: Proceedings of a
Symposium on Atmospheric Diffusion and Air Pollution.
Santa Barbara, American Meteorological Society, September
1974. pp. 387-391.
29. Turner, D.B. "Workshop on Atmospheric Dispersion Esti-
mates." Office of Air Programs, Environmental Protection
Agency, Publication No. AP-26 (5th printing), 1972.
30. Pasquill, F. Atmospheric Diffusion. 2nd edition.
New York, John Wiley & Sons, Inc., 1974.
31. Lewellen, W.S., M.E. Teske, and C.duP. Donaldson.
"Turbulence Model of Diurnal Variations in the Planetary
Boundary Layer." In: Proceedings of the 1974 Heat
Transfer and Fluid Mechanics Institute, Davis, L.R. and
R.E. Wilson (eds.). Stanford, Stanford University Press,
1974. pp. 301-319.
32. Panofsky, H.A. "The Atmospheric Boundary Layer Below
150 Meters." In: Annual Review of Fluid Mechanics.
Palo Alto, Annual Reviews, Inc., 1974. pp. 147-177.
33. Naot, D., A. Shavit, and M. Wolfshtein. "Two-Point
Correlation Model and the Redistribution of Reynolds
Stress." Physics of Fluids. l6_:738. June 1973-
34. Lamb, R.G., W.H. Chen, and J.H. Seinfeld. "Numerico-
Empirical Analyses of Atmospheric Diffusion Theories."
In: Proceedings of a Symposium on Atmospheric Diffusion
and Air Pollution. Santa Barbara, American Meteorologi-
cal Society, September 1974. pp. 317-324.
62
-------
35. Deardorff, J.W. "Numerical Investigation of Neutral and
Unstable Planetary Boundary Layers." J. Atmos. Sciences,
29.: 91-115, 1972.
36. Lumley. J.L. "Prediction Methods for Turbulent Flows."
The Pennsylvania State University. In a lecture series
presented at Von Karman Institute, Rhode-St,-Genese,
Belgium, March 3-7, 1975.
37- Townsend, A.A. The Structure of Turbulent Shear Flow.
Cambridge, Cambridge University Press, 1956,
38. Thomas, R.M. "Conditional Sampling and Other Measure-
ments in a Plane Turbulent Wake." J. Fluid Mechanics.
51 (Part 3):549-582, 1973.
39. Orszag, S.A. "Numerical Simulation of Turbulent Flows."
Massachusetts Institute of Technology. Presented at
University of Tennessee Space Institute Short Course on
Fundamentals and Applications of Turbulence, April 1975-
40. Wygnanski, I. and H. Fiedler. "The Two-Dimensional
Mixing Region." J. Fluid Mechanics. 4l_ (Part 1):327-
361, 1970.
4l. Birch, S.F. and J.M. Eggers. "A Critical Review of the
Experimental Data for Developed Free Turbulent Shear
Layers." In: Free Turbulent Shear Flows, Volume I.
NASA SP-321, 1973 (also see Volume II).
42. Pasquill, F. "Some Topics Relating to Modeling of
Dispersion in Boundary Layers." Report No. EPA-650/4-
75-015. Environmental Protection Agency, Raleigh, N.C.
April 1975.
63
-------
NOMENCLATURE
A, a, b model constants
C mean species concentration
c species concentration fluctuations about the
mean
D laminar diffusion coefficient
d model constant relating plume spread to plume
turbulent scale, see Eq. (15)
f Coriolis parameter = 2fi sin
g gravitational acceleration
k von Karman's constant = -u*/(z 3u/8z) in
the neutral surface layer °
L Monin-Obukhov length = -0ou*3/kg(w6")o
n decay rate for grid turbulence
P mean pressure (also Pasquill's stability
parameter)
p pressure fluctuation
q square root of twice the turbulent kinetic
energy
Q source strength
Ri Richardson number
Ri characteristic Richardson number defined in
Eq. (12)
Ro Rossby number based on geostrophic conditions
Ro characteristic Rossby number = u/fz
s model constant
S0 ' sl» S2 constants in scale equation
t time
t Brunt -Vaisala period
O
U. mean fluid velocity in the ith direction
U geostrophic wind velocity
o
u. fluctuating fluid velocity in the ith
direction
u* surface shear stress velocity
64
-------
u
vc
W
w
W*
z
z,
z,
z
characteristic surface velocity equal to 0.1
times the wind velocity at 10 meters
model turbulent diffusion coefficient
mean vertical velocity
fluctuating vertical velocity
free convective velocity scale =
[gz. (we")
1 ' " 'O O
Cartesian coordinates
vertical coordinate
temperature inversion height above an unstable
layer
surface roughness height
distance to 1/2 the centerline value
plume release height
e
0
6
A
X
v
P
Superscript
Subscript:
o
g
P
t
,1/2
dissipation function
mean potential temperature
fluctuating potential temperature
model macroscale length
model microscale length = A/[a + bqA/v]"
kinematic viscosity
density
latitude
angular velocity of the coordinate system
plume spread
denotes ensemble average
denotes surface value
denotes geostrophic condition
denotes plume value
denotes ambient turbulence value
65
-------
APPENDIX A
SUMMARY OF ONE-DIMENSIONAL UNSTEADY EQUATIONS AND
BOUNDARY CONDITIONS FOR DIFFUSION IN THE
PLANETARY BOUNDARY LAYER
As a convenience to the reader, the equations solved
for the mean flow variables and turbulence field are summar
ized here for the one-dimensional unsteady case. The mean
flow equations are:
5F=f(v- V
- V+vo ' Iz™ ' (A'2)
- „ ,T7T f*
™r ~ V 5- - -r— W0 (A
Dt O * d dZ
To these are added the equations for the six independent
Reynolds stresses:
_UU. _ Jini^y gj_n _ uw Cos <|>) - 2UW ^ + V •» —
Aq
8uu~|
3z
,2 _
n- UU
2v
66
-------
Dvv
Dt
sin
2vw 9V
3z
V
c 3?
._
Aq
H--V)
32 —
+ v —5- vv
0 3z2
2v
o 2
3X'
(A.5)
Dww _ )iri—
Dt
cos <(>
G
wT + v a
O
Aq
3ww
3z
V
2v
o
0 3z2
ww -
o
3X'
(A.6)
Duv
Dt
= 2fl(w sin - vw cos <|) - uu sin + uu cos ) - ww ^— + •£- u8
Dt 9z 0
oZ
uw (A.8)
67
-------
Dvw
Dt
gy ~
sin + uv cos ) - ww ^— + 7?- v0
d Z U
o
V
A~
Aq
- A VW
+ v —T- vw
0 3z2
(A.9)
the equations for the three components of the heat flux vector
DU6 _ ? ,—
Dt *
sin 4 - w0 cos ) - uw -r— - w6
3z
+ V
c 3z
Aq
9u6
3z
v
9z
(A. 10)
Dve = ^rr^
Dt
— 30 -TT
sin - vw •»— - w6
v
Aq
9^0 "]
9z J
(A.11)
DW6 _ n^^7T
Dt
cos <(> - ww TT-
+ v
c 3z
Aq —=•
3w9~|
3z I
-f
(A.12)
68
-------
the equation for the temperature variance:
—
T-\Q ^ 'N (->
BE- ' -25F H
3z
+ v
2 2
38
^-~k
0 3z2
- 2 v s ^5-
o •> <-
A
(A.13)
and the equation for the turbulent macroscale:
2 X2
gw9
(A.14)
These fourteen equations are solved for the turbulence
field and the results used as inputs into Eqs. (1) through
to solve for the species distribution.
The general boundary conditions applied to these equa-
tions at the top of the boundary layer are that as z •*• °°
U = u ^
^ > geostrophic wind
V = v
gJ
3 0
g— = constant (difference between the
adiabatic lapse rate and the
ambient lapse rate)
u u. = u.6 = 69 = 0
i j i
Rather than compute all the way to the surface with
these equations, which would, in general, require computing
the flow around various roughness elements, the surface
boundary conditions are applied at zo , the effective aero-
dynamic roughness height. At
z = z
o
U = V = 0
0 = Surface (°r w6 = w6surface)
69
-------
3U1UJ _ 9U19 _ 869- _
3z 3z 3z
The flow is now completely specified when initial condi-
tions on all the variables are specified. The dimensionless
parameters governing the flow may be obtained by normalizing
the equations and boundary conditions. For a characteristic
velocity, let us choose Ug , with the coordinate system
aligned so that Vg = 0 . Then, if we normalize all lengths
by z0 and time by f"1 , Eqs. (A.I) and (A.2) introduce two
parameters: a Reynolds number UgZQ/v and a Rossby number
Ug/zof . For the high Reynolds numbers applicable for
atmospheric flows, the Reynolds number dependence can usually
be ignored. When the flow is neutral (i.e., 0 = constant
everywhere), then the only other parameter introduced by the
remaining equations or the boundary conditions is cot $ .
If temperature is normalized by A0, not yet specified,
then the equations introduce one more parameter: a bulk
Richardson number, Ri^ = gA6zo/0oUg . For stable conditions,
A6 may be taken as the temperature difference across the
boundary layer, while for unstable conditions, it would
appear reasonable to take it as the height of the inversion
layer times Oe/Sz)^ . However, there is a rather strong
dependence on initial conditions involved when the influence
of stability is introduced. In fact, our previous work6
shows that for typical daily variations in the planetary
boundary layer, a hysteresis-type pattern is followed for the
turbulence field. In spite.of this complication, in Section
VI we have considered the quasi-steady conditions which will
exist after conditions of roughly the same stability have
been maintained for some time. Under stable conditions, such
a quasi-steady turbulence field may be characterized by a
single stability parameter. In Section VI, we choose a bulk
Richardson number based on the temperature and velocity
differences between the surface and 10 meters. However,
unstable conditions appear to require two stability parameters
even under quasi-steady conditions. To the bulk Richardson
number, we add the ratio of the inversion height, z^ , to u/f ,
For this purpose, z^ is defined as the altitude near the top
of the boundary where the temperature variance reaches a local
maximum (see Fig. B.5).
70
-------
APPENDIX B
AIAA PAPER
No. 75-871
EXAMPLES OF VARIABLE DENSITY FLOWS COMPUTED BY A
SECOND-ORDER CLOSURE DESCRIPTION OF TURBULENCE
by
W. S. LEWELLEN, M. E. TESKE,
and
COLEMAN duP. DONALDSON
Aeronautical Research Associates of Princeton, Inc
Princeton, New Jersey
MM Bin Fluid and Plasma
Dynamics conference
HARTFORD, CONNECTICUT/JUNE 16-18, 1975
For permission to copy or republlsh, contact the American Institute of Aeronautics and Astronautics,
1290Avenueof the Americas, New York, N.Y. 10019.
71
-------
EXAMPLES OF VARIABLE DENSITY FLOWS COMPUTED BY A SECOND-ORDER CLOSURE
DESCRIPTION OF TURBULENCE*
W. S. Lewellen,** M. E. Teske,"1" and Coleman duP. Donaldson"1"1"
Aeronautical Research Associates of Princeton, Inc.
Princeton, New Jersey
Abstract
An incompressible fluid, with the influ-
ence of the variable density coming through
a Boussinesq buoyancy term, is considered
using a second-order closure approach to
turbulent flow. This technique successfully
incorporates the strong influence of strat-
ification for a number of flows. Particular
flows which are computed and compared with
experimental data are: a) stratified shear
flow with homogeneous turbulence; b) free
convection bounded by a stable density
gradient; c) wakes in a medium with a
stable density gradient; and d) the atmos-
pheric boundary layer where the flow
oscillates from stable to unstable condi-
tions in a typical day.
Nomenclature
F,. wake Froude number based on diameter
g, general acceleration vector
k coefficient of laminar diffusion for
temperature perturbation
pressure fluctuation
mean pressure
turbulent Prandtl number
square root of twice the turbulent
kinetic energy
Richardson number
time
Brunt-Vaisala period; defined in (12)
surface shear velocity
fluctuating velocity
mean velocity
P
P
Pr\
Ri
t
'c
u*
ui
Ui
w*
'50
free convection velocity scale;
defined in (11)
direction
temperature inversion height above an
unstable layer
distance to one-half the centerline
velocity value
energy dissipation
fluctuating temperature
0 mean temperature
K von KSrma'n constant
A macroscale length
v coefficient of viscosity
p density
n. angular velocity of the coordinate
3 system
superscript denotes time average
subscript a denotes ambient value
1. Introduction
Small changes in density are sufficient
to have a strong effect on turbulent flows
in the presence of a gravitational acceler-
ation. Such flows are extremely difficult
to predict with eddy viscosity approaches.
This paper describes a second-order closure
approach to turbulent flow which has been
found to successfully incorporate the influ-
ence of stratification for a number of
different flows. Our approach is to use
the Reynolds stress equations to compute
the turbulent fluctuations with invariant
modeling techniques applied to close the
system of equations. Details of this
approach have been given by Donaldson.1'2
The basic assumption is that the dependence
of third-order correlations of the velocity
fluctuations on the second-order correla-
tions, the mean flow properties, and their
derivatives is invariant with respect to
changes in flow geometry. This permits
data from relatively simple flow experiments
to be used in evaluating the necessary
coefficients in the modeled terms.
Similar approaches to stratified flows
have been taken by Mellor,3 Daly and
Barlow," Wyngaard et al., and Launder.6
However, the particular model representa-
tions of the terms required for closure of
the Reynolds stress transport equations are
somewhat different in each case. Limited
comparisons between models have been made
by Meroney.7 This paper will not give
comparisons between various models but will
*This research has been partially funded with Federal funds from the Environmental
Protection Agency under contract number EPA 68-02-1310 and the Defense Advanced Research
Projects Agency of the Department of Defense as monitored by ONR under contract number
N0001i)-72-C-Ol)13. The content of this paper does not necessarily reflect the views or
policies of any agency of the U.S. Government.
**Senior Consultant; Member, AIAA
t Associate Consultant
^President and Senior Consultant; Associate Fellow, AIAA
Copyright® American Institute of Aeronautics and I '-
Astronautics, Inc., 1975. All rights reserved.
-------
be restricted to a brief description of a
model and its application to a number of
different flows involving buoyancy.
2. Model Equations
- v6
(5)
When the flow variables are decomposed
into mean and fluctuating quantities, the
mean flow equations may be written as
3U.
3U
3u7u
T
.
J
3U.
15
dX.
3t
3X
2 _ 36 36
J ~ dK 3x 3x,
J J
"
(6)
P 3x,
3U,
= 0
3x
J
3Ui8
3x
30
(1)
(2)
(3)
The flow has been assumed incompressible
but with small variations in density due to
changes in temperature. With the Boussin-
esq approximation, the effect of the den-
sity variation enters only in the gravita-
tional body force term in the momentum
equation, Eq. (1). The last term in (1) is
due to coriolis forces present in a coordi-
nate system rotating with angular velocity
f2 . It is included here in anticipation of
atmospheric applications. Equation (3) is
a diffusion equation for the temperature
perturbation.
This system of equations is not complete
due to the presence of the Reynolds stress
terms. Exact equations for these correla-
tions may be derived as outlined by
Donaldson.
3U.
- UiUk
3U
g u e g u e
0 + 0
3ue
fl
- f H- +
Closure of the system of equations at
this second-order level requires modeling
the last four terms in Eqs. (4) and (5) and
the last two terms in Eq. (6). The terms
required for (4) have been previously dis-
cussed in Refs. 2 and 8. These remain un-
changed for a variable density flow; the
terms required for Eqs. (5) and (6) are
modeled in a similar fashion. The two new
coefficients introduced with these variable
density terms were evaluated using data from
the atmospheric surface layer.9 With the
introduction of the model terms and the
numerical coefficients, Eqs. CO-(6) are
represented for high Reynolds number flow as
3U
- uiuk
3U
g u 9
2EJJlkVkUi
3u16
3t
3iTJT
3Xj
(7)
30
3U,
g.e' a / 3u e\
-V- 2eijkYke + °-3l7:'qA-i-i
(8)
37
3t
362
0.3
3x,
(9)
To complete closure, it is still necess-
ary to determine the model turbulent scale
A . It could be specified empirically based
on gross features of the particular flow
geometry.2>8 Here it will be predicted
from a semi-empirically modeled, dynamic
73
-------
differential equation. As Bradshaw10 and
Mellor and Herring" have pointed out, the
form of the A equation is similar,
whether one begins with the two-point velo-
city correlation to define an integral
scale, the dissipation function e = q3/8A ,
or the vorticity fluctuations q/A . The
principal difference is in the diffusion
terms. The form of the A equation used
here is
DA _
3U,
0.075q
0.3
0.8A giui6
8A
3x.
0.375
q
(10)
A detailed accounting of these coefficients
may be found in Ref. 12. The coefficient
of the second term on the right-hand side
of (10) is set by comparison with turbu-
lence decay behind a grid; that of the
third term by analogy with the diffusion
terms in (7)-(9); that of the last term by
comparison with atmospheric surface layer
data; and the other two by computer fit to
a number of different simple flows.
Equation (10) cannot command as much
confidence as (7)-(9) because of the need
to model all of the terms in the equation.
This is true whether dealing with equations
for e , A , or q/A . It is also necess-
ary to impose boundary conditions appropri-
ate to the particular flow situation.
3. Example Calculations
With the model set, verification may be
achieved by comparing predictions and
observations for example flows which were
not involved in the determination of the
model coefficients. Several such example
flows are considered here.
3.1 Temperature fluctuations in a two-
dimensional wake
The modeling of Eqs. (8) and (9) may be
independently checked without involving
stability influences by calculating the
temperature fluctuations in a plane turbu-
lent wake when the fluctuations are suffi-
ciently weak for the gravitational terms to
be neglected. Measurements behind a
slightly heated cylinder have been made by
Freymuth and Uberoi13 and LaRue and Libby."1
Figure 1 shows a comparison between the
measured distributions far downstream of the
cylinder and the self-similar profiles
predicted by the model. The normalizing
length ZC-Q is 'fixed by comparison with the
mean velocity profile. The agreement is
quite good, particularly in the concurrence
of the amplitude and position of the maxi-
mum temperature fluctuation. The center-
line value of (62)1/2/A0 is given as 0.28
by LaRue and Libby in contrast to the value
of 0.21 reported by Freymuth and Uberoi.
The present model predicts a self-similarity
1.6
1.4
1.2
1.0
.4
o A LaRue-Libby
data
A.R.A.P. prediction
.4
1.2
Z/Z 50
1.6
2.0
2.4
Figure 1. Temperature similarity profiles
of mean temperature increment and tempera-
ture variance in a two-dimensional wake;
data from Ref. 14
value of 0.276.
3.2 Stratified homogeneous turbulent shear
flow
If an equilibrium flow with constant
gradients in mean velocity and temperature
can be established so that the turbulent
correlations are homogeneous, Eqs. (7)-(9)
reduce to algebraic relations between the
turbulence and the mean flow gradients.
Webster15 attempted to produce such a homo-
geneous, stratified flow in a heated wind
tunnel. Apparently, he did not achieve
this goal completely since his turbulent
correlations show a dependence on the stream-
wise direction. In spite of this acknowl-
edged deficiency, we compare in Figure 2 the
turbulent Prandtl number dependence on
stratification predicted by Eqs. (7)-(9)
with that reported by Webster at his most
downstream station. The prediction of
Launder's second-order closure model6 is
also indicated. Launder uses Webster's data
to estimate some of his model coefficients.
The agreement between theory and observa-
tions in Figure 2 should only be interpreted
as qualitative verification of the depen-
dence on stability. The Prandtl number is
normalized by its neutral value because of
-------
3.Or
2.5
2.0
Pr,/Pr,
i.o
0.5
o Data of Webster
Superequilibrium predictions
A.R.A.R model
Launder model
.2
.3
.4
Ri
Figure 2. Turbulent Prandtl number depen-
dence on stratification in a homogeneous
stratified fluid. Data from Ref. 15;
comparison prediction from Ref. 6
the large uncertainty in its reported values
by Webster, varying by a factor of 3 from
measurements at one station to the next.
Reliable turbulence data under stable
conditions are difficult to obtain. Proba-
bly the best measured conditions are those
existing in the atmospheric surface layer.
The present model quite reasonably predicts.
the observed variations of the normalized
variables with respect to the similarity
variable, the ratio of the vertical height
to the Monin-Obukhov length.9 However, this
cannot be used as independent verification
of the model because of the use made of
these data in determining the model
coefficients.
Another stable flow which has been pre-
dicted by the model is that of shear layer
entrainment in a stratified fluid. Calcu-
lations were made to simulate the experiment
of Kato and Phillips.16 They used dye to
visually measure the rate of spreading of
shear-generated turbulence into a region of
fluid with a stable density gradient.
Comparisons between the present model pre-
dictions and the observations have been
made in Ref. 12. Since this represents.
little more than a qualitative comparison,
it is not repeated here.
3.3 Free convection
Free convection is a flow that is quite
different from the flows used to determine
the model coefficients. The turbulence is
produced by the unstable buoyancy force
without influence from any mean shear stress.
Willis and Deardorff17 have measured the
velocity and temperature fluctuations in an
unstable layer bounded above by a stable
density gradient and below by a surface with
a positive heat flux. Their experiment
called for establishing an initial stable
temperature gradient between two surfaces
and then applying a heat flux to the region
through the lower surface. With mean
velocities restricted to a minimum, free
convection occurs in a mixed layer above the
lower surface. The thickness of this mixed
layer increases with time, but the fluctua-
ting velocity distributions exhibit self-
similarity when normalized by the character-
istic velocity
„* =
(11)
where w6 is the surface heat flux and z^
is the depth of the mixed layer.
Figure 3 compares our model predictions
for the normalized vertical velocity fluc-
tuations with the experimental observations.
For the model predictions, we assumed no
mean velocity and a constant surface heat
flux to simulate the laboratory experiment.
The agreement between predictions and obser-
vations is heartening, particularly since no
model constant adjustments were involved.
.75
Z/Z;
.25
.2
.4
WW/W*2
.6
Figure 3- Similarity profiles of vertical
velocity fluctuation in an unstable mixed
layer; data from Ref. 17
-------
The horizontal velocity fluctuations are
given in Figure 4. Here the comparison with
experiment is not as close since the obser-
vations show a peak near the lower surface
that is not evidenced in the model results.
The distributions of temperature fluctua-
tions are shown in Figure 5- The agreement
between model and experiment is very good
except near the top of the mixed layer
where the observations show a much stronger
local maximum than is predicted. This may
be due to the fact that inertial oscilla-
tions exist in this region. Overall, the
comparison between model predictions and
experimental observations is quite favorable.
A •
I/I
.75
Z/Zj
.5
.25
Figure 5- Similarity profiles of tempera-
ture fluctuations for the conditions of
Figure 3
• A
.2 .4
HU/W*2
.6
Figure 4 . Similarity profiles of horizontal
velocity fluctuations for the conditions of
Figure 3
The favorable comparisons of Figures 3-5
are in spite of the fact that the modeled
vertical flux of kinetic energy is signi-
ficantly different from the experimental
observations, as shown in Figure 6. This
is consistent with our basic assumption that
the second-order correlations calculated
from their conservation equations are given
more accurately than the accuracy to which
the third-order correlations are approxi-
mated.
3.lj Diurnal variations in the planetary
boundary layer
Over homogeneous terrain, the atmos-
pheric boundary layer may be considered a
function of time and altitude only. The
results of calculations for a fixed geo-
strophic wind and upper level temperature
lapse rate with a cyclic surface heat flux
.75
Z/Z
1 0 .1
.2
wq*/2w*3
Figure 6. Similarity profiles of vertical
flux of kinetic energy for the conditions
of Figure 3
76
-------
approximating conditions over the Midwest-
ern plains during summer18 are presented in
Ref. 19. These results were obtained, using
a quasi-equilibrium approach to the equa-
tions and a gross feature approach to the
scale. The results presented here were
obtained using the full equations as presen-
ted in Section 2.
The predicted contours of constant tur-
bulent fluctuations are presented in Figure
7 as a function of time and altitude. The
boundary layer thickness grows during the
day and shrinks during the evening and
morning hours. The turbulent kinetic energy
reaches a maximum of 4.8? of the geostrophic
mean kinetic energy in the midafternoon at
an altitude of approximately 500 meters.
As the sun sets, the turbulence begins to
decay until in the early morning hours its
maximum kinetic energy is sO.25? of the
mean kinetic energy. The biggest differ-
ence between the results presented here and
those in Ref. 19 occurs during the early
morning hours. The full equations with a
is
Figure 7. Isopleths of turbulent energy
in a diurnal planetary boundary layer
dynamic scale predict a slower decay of q2
and consequently considerably larger q2 in
the altitude range from 100 m to 1 km in
the morning hours. But even this larger
value is still quite small compared with
afternoon values. Both model representa-
tions predict such features as the tempera-
ture inversion and local wind jet observed
to occur nocturnally at low levels.
The variation of the angle between the
direction of the geostrophic wind above the
boundary layer and the direction of the
surface wind is given in Figure 8 as a
function of time of day. Some data points
averaged over a summer for another location
by Hoxit20 are also included as a qualita-
tive comparison. Since the average includes
days with atypical surface heat flux and
other influences such as baroclinicity,
geostrophic wind variation, and moisture
transport between the surface and the
atmosphere, it is consistent to expect the
model to predict a wider spread between the
large angle at night and the small angle
around noon. It is not clear whether there
3 6
IZ IS 1821 2«
Figure 8. Diurnal variation of surface wind
angle; data from Ref. 20
is sufficient shift in the surface heat flux
variation to account for the apparent shift
of approximately 3 hours in the rapid trans-
ition from high to low angle at sunrise.
Surface shear velocity is plotted as a
function of the stability variable KuVfL
in Figure 9. Data points for Ro s 10?, as
taken from Tennekes' summary,21 are included
in the plot. The model predictions show a
hysteresis loop with the surface shear stress
significantly larger when the surface heat
flux is decreasing rather than increasing.
The data show considerable scatter but tend
to be biased towards the upper bound. The
factor of two difference in u* at neutral
conditions demonstrates the influence of
unsteady dynamics on the atmospheric bound-
ary layer.
•5 ,
-30
-60
/cuVfL
Figure 9. Diurnal variation of surface
shear velocity as a function of stability;
data from Ref. 21
77
-------
3.5 Stratified wake
The passage of a body through a strati-
fied fluid generates a turbulent wake
containing kinetic and potential energy.
The buildup of potential energy is caused
by the co-mixing of heavier density fluid
above the lighter density fluid. In the
initial stages of wake development in a
weakly stratified medium, the turbulent
kinetic energy dominates the potential
energy, and the wake spreads. This spread-
ing in turn increases the potential energy
until, at a time comparable to the natural
oscillating period of the fluid, inertial
waves are generated. Computations for this
flow using a quasi-equilibrium version of
the present model have been reported in
Ref. 12.
Figure 10 shows typical contours of the
turbulent kinetic energy and secondary flow
stream function at two instances of time
following wake initialization. Time is
normalized by the characteristic Brunt-
Vaisala period of the fluid
2Tr[(g/0)30/8z]
-1/2
(12)
l/tc=0.5
q2mo, =0.0000534 uj
l/lc=0.5
t/lc=I.O
q2mox = O.OOOOI47llJ
l/tc=I.O
Figure 10. Contours of turbulent energy
and secondary flow for two instances of
time after wake initialization for Ri0 =
0.00925. 1 = +10% of maximum value for that
quadrant; 2 = -1058; 3 = +30%; 1 = -30%; etc.
The flow is symmetric about both axes so
that only one quadrant is shown for each
variable at each time. At t/tc = 0.5, the
secondary flow exhibits a simple collapsed
motion, while at t/tc = 1, a second wave
mode has been added. At both times the
secondary flow extends well outside the
region of turbulence.
Figure 11 is a summary plot showing the
sensitivity of the wake shape to initial
Froude number, FD = Utc/ D , or initial
Richardson number of the turbulence, RiQ =
(2nr1/qmaxtc) , along with a qualitative
comparison with laboratory flow visualiza-
tion results. 22
Figure 11. Vertical spread of stratified
wake; data from Ref. 22
1. Concluding Remarks
The explicit appearance of buoyancy terns
in the equations for Reynolds stress and
turbulent heat flux are sufficient to
account for many of the effects observed in
variable density flows. The only direct
influence of buoyancy included in the model-
ed terms required for closure is its influ-
ence on the macroscale parameter A . The
predictions of this relatively simple
second-order closure model appear to be
reasonably consistent with available obser-
vations, although there still exists a need
for unambiguous experiments under relatively
stable flow conditions. One can expect
continued model development, but the model
appears sufficiently accurate to permit
practical calculations to proceed concurr-
ently.
5. References
iDonaldson, C.duP., "Calculation of Turbu-
lent Shear Flows for Atmospheric and Vortex
Motions," AIAA Journal, Vol. 10, No. 1,
January 1972,~pp. 4-12.
2Donaldson, C.duP., "Atmospheric Turbulence
and the Dispersal of Atmospheric Pollutants,"
AMS Workshop on Micrometeorology (D.A.
Haugen, ed.), Science Press, Boston, 1973,
PP. 313-390
3Mellor, G.L., "Analytic Prediction of the
Properties of Stratified Planetary Surface
Layers," J. Atmos. Sciences. Vol. 30, 1973,
pp. 1061-1069.
^Daly, B.J. and Harlow, F.H., "Transport
Equations in Turbulence," Phys. Fluids,
Vol. 13, No. 11, 1970, pp. 2631-2619.
5Wyngaard, J.C., Cote, O.R. and Rao, K.S.,
'Modeling the Atmospheric Boundary Layer,"
Advances in Geophysics, Vol. 18A, Academic
Press, New York, 197.1.
6Launder, B.E., "On the Effects of a Gravi-
tational Field on the Turbulent Transport of
Heat and Momentum," J. Fluid Mech.. Vol. 67,
Part 3, 1975, pp. 569-581.
'Meroney, R.N., "Buoyancy Effects on a Turb-
ulent Shear Flow," Project THEMIS Report No.
28, Colorado State University, April 1971.
8Lewellen, W.S., Teske, M.E., and Donaldson,
C.duP., "Application of Turbulent Model
Equations to Axisymmetric Wakes," AIAA
Journal, Vol. 12, No. 5, 1971, pp.~6~20"-625.
78
-------
9Lewellen, W.S. and Teske, M.E., "Prediction
of the Monin-Obukhov Similarity Functions
from an Invariant Model of Turbulence,"
J. Atmos. Sciences. Vol. 30, 1971, pp. 1340-
iWT.
10Bradshaw, P., "The Understanding and
Prediction of Turbulent Flows," Aeronautical
Journal T 1972, pp. 403-418.
TTMellor, G.L. and Herring, H.J., "A Survey
of Mean Turbulent Field Closure Models,"
AIAA Journal, Vol. 11, No. 5, 1973, pp. 590-
599.
•"•^Lewellen, W.S., Teske, M.E., and Donaldson,
C.duP., "Turbulent Wakes in a Stratified
Fluid," Aeronautical Research Associates of
Princeton, Inc. Report No. 226, September
1974.
13Freymuth, P. and Uberoi, M.S., "Structure
of Temperature Fluctuations in the Turbu-
lent Wake Behind a Heated Cylinder," Phys.
5TC
Fluids. Vol. 14, No. 12, 1971, pp. 25
1ALaRue, J.C. and Libby, P. A., "Temperature
Fluctuations in the Plane Turbulent Wake,"
Phys. Fluids, Vol. 17, No. 11, 1974, pp.
1956-1975.
15Webster, C.A.G., "An Experimental Study
of Turbulence in a Density-Stratified Shear
Flow," J. Fluid Mech., Vol. 19, Part 2,
1964, pp. 221-245.
16Kato, H. and Phillips, O.M., "On the
Penetration of a Turbulent Layer into Strat-
ified Fluid," J. Fluid Mech.. Vol. 37,
Part 4, 1969, pp. 643-656.
17Willis, G.E. and Deardorff, J.W., "A
Laboratory Model of the Unstable Planetary
Boundary Layer," J. Atmos. Sciences, Vol.31,
1974, pp. 1297-1307.
18Wyngaard, J.C., "Notes on Surface Layer
Turbulence," AMS Workshop on Mlcrometeor-
ology (D.A.Haugen, ed.), Science Press,
Boston, 1973, pp. 101-149.
19Lewellen, W.S., Teske, M.E. , and Donaldson,
C.duP., "Turbulence Model of Diurnal Varia-
tions in the Planetary Boundary Layer,"
Proc. 1974 Heat Transfer and Fluid Mechan-
ics Institute (L.R.Davis and R.E.Wilson,
eds. ) , Stanford University Press, Stanford,
1974, pp. 301-319.
20Hoxit, L.R., "Diurnal Variations in
Planetary Boundary-Layer Winds Over Land,"
Boundary Layer Meteorology, Vol. 8, 1975,
BE.- 21-38.
/JTennekes, H. , "Similarity Laws and Scale
Relations in Planetary Boundary Layers,"
AMS Workshop on Mlcrometeorology (D.A.
Haugen, ed. ; , Science Press, Boston, 1973,
pp. 177-216.
22Lin, J.T. and Pao, Y.H., "Turbulent Wake
of a Self-Propelled Slender Body in
Stratified and Nonstratified Fluids: Analy-
sis and Flow Visualizations," Flow Research
Corp. Report No. 11, 1973.
79
------- |