EPA-R4-73-016a
March 1973
Environmental Monitoring Series
Atmospheric Turbulence
and the Dispersal
of Atmospheric Pollutants
Office of Research and Monitoring
U.S. Environmental Protection Agency
Washington, D.C. 20460
-------
EPA-R4-73-016a
Atmospheric Turbulence
and the
Dispersal of Atmospheric
Pollutants
by
Coleman duP. Donaldson
Aeronautical Research Associates of Princeton, Inc.
50 Washington Road
Princeton, New Jersey 08540
Contract No. 68-02-0014
Program Element No. A-11009
EPA Project Officer: Kenneth L. Calder
Meteorology Laboratory
National Environmental Research Center
Research Triangle Park, North Carolina 27711
Prepared for
OFFICE OF RESEARCH AND MONITORING
U. S. ENVIRONMENTAL PROTECTION AGENCY
WASHINGTON, D.C. 20460
March 1973
-------
This report has been reviewed by the Environmental Protection Agency
and approved for publication. Approval does not signify that the
contents necessarily reflect the views and policies of the Agency,
nor does mention of trade names or commercial products constitute
endorsement or recommendation for use.
11
-------
PREFACE
In February 1971, Aeronautical Research Associates of
Princeton, Inc. (A.R.A.P.) was awarded a small contract by the
Environmental Protection Agency to assess the practicality of
developing a theory for the dispersal of pollutants by the
atmosphere - a theory which would be more fundamental than
existing eddy diffusion or K methods. The method to be used
was a second-order closure scheme under development at A.R.A.P.
by Dr. Coleman duP. Donaldson and his colleagues. Since it
appeared, on first application of this method, that a powerful
new technique might be successfully developed for computing the
dispersive power of the atmosphere under arbitrary meteorological
conditions, a further contract was awarded to A.R.A.P. to continue
the development of the method. A portion of this extension of
the effort was to be applied towards a partial funding of the
preparation of a monograph which would bring together a connected
account of the derivation of the appropriate equations for the
model of atmospheric turbulence and transport. The task is a
difficult one for our work is not complete. Although the basic
equations will not change, it is expected that, as our investiga-
tion proceeds, changes will be made in the nature and/or complex-
ity of the models that are used to obtain a closure of the basic
equations. A great deal remains to be done before a completely
satisfactory model with associated programs is in hand. It is
felt, however, that a working paper which gives a detailed
progress report on the present state of A.R.A.P.'s atmospheric
model will be useful. It should be pointed out that the work
reported here has been the result of funding from many sources,
including internal corporate support. Particular acknowledge-
ment should be made of the support received from the National
Aeronautics and Space Administration and the Air Force Office of
Scientific Research, in addition to the Environmental Protection
Agency.
ii
-------
TABLE OP CONTENTS
Preface
Nomenclature
1. Introduction
2. Equations for an Atmospheric Shear Layer
3. Equations for Turbulent Shear Layers
l). Selection of Models
5. Search for Model Parameters
6. An Equation for the Scale A,
7. Relation of Second-Order Models to K Theory
8. Application of Second-Order Modeling to the
Atmospheric Boundary Layer
9. Dispersal of Pollutants from a Line Source
10. Dispersal of Pollutants from a Point Source
11. Critique of Atmospheric Turbulence and Pollution
Dispersal Calculations
12. Invariant Modeling of the Ekman Layer
13. Conclusions and Recommendations
14. References
Appendices
A. Numerical Integration Technique
iii
-------
NOMENCLATURE
a, b, C,, Cp, GO> c', cl constants
A, B, Q general dependent variables
A,, constant tensor
C, C..,, C,, C,. constants
-L > <- J "
C general species concentration
C specific heat at constant pressure
ft- diffusion coefficient
f Coriolis parameter, 2P, sin <(>
g acceleration due to gravity
g. general acceleration vector
h enthalpy
k conductivity
L local integral scale; Monin-
Obukhov length scale
L typical atmospheric length scale
3.
L typical motional length scale
s
M Mach number
p pressure
P superequilibrium solution para-
meter defined in Eq. (7.45)
Pr = yC /k Prandtl number
q =
r, <|> , z cylindrical coordinate system
R gas constant
Re = p u.L /u Reynolds number
O J S O
- //a-.2
Ri = |— -g— / (•9^) Richardson number
o
Ri ., = iiKM + Qh'l critical Richardson number for
superequilibrium flow
Sc = pid/y Schmidt number
t time
T temperature
T.., T. general.tensors of rank two and
J ^-J^ three, respectively
iv
-------
U, V, W
u.
UU, W, WW, etc
'k
a
3
Y
a
6
6
A
e
n
K
X, A, A
char
a
u
u = u/p
T = - pu'w'
Subscripts
A, B
c
char
g
J, k
n
t
o
velocity components
general velocity vector
nondimensional representations of
u'u1, v'v
w'w'
etc . (c .f. ,
Eq. (7.30))
Cartesian coordinate system
general coordinate system
average length vector
constant slope of A, near z = 0
constant
ratio of specific heats
spread
characteristic mean motion scale
Kronecker delta
difference
small distance
constant
von Karman's constant
scalar length scales
first and second coefficients of
viscosity, respectively
kinematic viscosity
length difference vector
density
stress tensor
turbulent shear stress
latitude
Earth's angular rotation
general,vector position notation
concentration correlation value
characteristic value
geostrophic wind condition
distance discretization designation
time discretization designation
background turbulence value
atmospheric equilibrium condition
v
-------
1. INTRODUCTION
In 1967, the author, in order to study some aspects of
the problem of boundary layer transition, constructed a model
of transition [Ref. 1] by a rather simple second-order closure
of the time-averaged equations for an incompressible fluid
whose motion consists of a mean and a fluctuating part. While
the study of the growth of disturbances during transition was
of some interest in itself, the aspect of the transition model
that most intrigued the author was the fact that it produced,
if let run long enough so that transition was complete, a
turbulent boundary layer that had some of the characteristics
of real boundary layers. This result started the writer and
his colleagues on an effort to construct a viable computational
model of turbulent shear layers along the general lines of the
model of transition. It was recognized that the method of
approach was not new. Indeed, the idea was surely understood
by Reynolds, who first formulated the equations upon which it
is based [Ref. 2]. At that time, however, it was quite beyond
practical application since no computational technology was
available to provide solutions to the equations involved, even
had a closure scheme been formulated.
In 19^2, Kolmogorov [Ref. 3] and in 19^5, Prandtl and
Wieghardt [Ref. 4] specifically formulated the method along
the lines that are pursued by a number of investigators at the
present time. Two of the earliest to follow these leads were
Chou [Ref. 5] and Rotta [Ref. 6]. All of these early formula-
tions and discussions fell short of any practical results for,
again, computational facilities and techniques were not avail-
able to solve the sets of equations which resulted.
With the advent of the 1960's, the development of digital
computers had advanced to the point where one could reasonably
expect to develop a practical means of computing turbulent
shear layers by means of some sort of second-order closure
techniques. A rather large group of related methods appeared
within a short time from 1965 through 1970 [Refs. 7 through 15].
-------
1-2
Since 1970 the development of computer programs which exploit
some form of second-order closure scheme to compute turbulent
shear flows or boundary layer transition has proceeded apace,
and no attempt will be made here to summarize 'these efforts.
In January of 1970, the author and Dr. Harold Rosenbaum
presented a paper on the generation of atmospheric clear air
turbulence [Ref. 16]. In this paper, the first invariant model
of an incompressible turbulent shear layer developed at A.R.A.P.
was extended to the case of a thin atmospheric shear layer with
arbitrary stability. It was mentioned in that paper that there
would be no real difficulty in extending the model and program
that had been developed to permit the calculation of the
dispersal of passive atmospheric pollutants in the atmospheric
boundary layer under arbitrary stability conditions.
Subsequently, A.R.A.P. received moderate support from both
EPA and NASA to develop programs for computing the dispersal of
pollutants in the atmospheric boundary layer. The support
from EPA was directed towards the solution of the dispersal of
pollutants from a line source, while that from NASA was directed
towards the dispersal from a point source of pollutant material.
While these efforts were underway, additional studies at
A.R.A.P. were carried out; these were aimed at further under-
standing and improvement of the models used in the closure of
the equation of motion at second order.
In this paper we. will present the derivation of the set of
equations on which the A.R.A.P. model of atmospheric turbulence
and transport is based and, in addition, discuss the process
which was used to obtain the present closure model and the values
of the basic parameters which appear therein. Following these
discussions, the relationship of the present second-order closure
method for computing turbulent transport to eddy diffusivity or
K theory models will be discussed by showing that what is essen-
tially K theory can be obtained from a second-order closure
scheme as a well-defined limiting case. Finally, some typical
examples of computations using the present model will be presented
-------
2-1
2. EQUATIONS FOR AN ATMOSPHERIC SHEAR LAYER
The equations which we shall take to govern the motion of
a compressible, nonchemically-reacting perfect gas are:
the perfect gas law
p = pRT (2.1)
the continuity equation*
J
the momentum equation
3pu.
where the stress due to molecular diffusion is given by
3u
the energy equation
x K 9P j. 3 v
+ u + k
. /s,- N
3t + 3T7 (pujh) =
j
and the species conservation equation
t TT . = x: r a (2-6>
J J < J '
In what follows, we will be interested in the development of
turbulence and the transport of matter in a thin layer of the
atmosphere. We will assume that the matter to be transported in
the atmosphere is not too highly concentrated, so that it makes
no first-order effect upon the heat capacity or gas constant of
the air in which it is carried. In this case, if we designate
the heat capacity as Cn and the gas constant as R^ , we may
r*o O
*In this equation and those that follow, the summation convention
of tensor analysis is used; i.e., repeated indices in any expres-
sion indicate that one is to sum over these indices.
-------
2-2
write (2.1) as
and (2.5) as
CP
"n
9x
P = PR0T
|2- +
9x.
J
3x.
J
k 3T
9x
9u.
(2.7)
(2.8)
Since the Mach number M of the flow in which we will be
interested is small, we may neglect the last term on the right-
hand side of (2.8) since this term represents the heat generated
2
by the dissipation of the motion and is of order M compared
with the other terms in the equation. The final form of the
energy equation is, then,
o
3t
9t
9x
.
J
9
9x
k
9x
(2.9)
Following the usual practice for obtaining the equations
for the motion of an atmospheric shear layer [Refs. 17 and 18],
we will consider the atmosphere to be in a state slightly removed
from an adiabatic atmosphere at rest. We consider then an expan-
sion of the equations presented above according to the following
scheme:
P = P0 + P
P
T
u .
J
P0 + P
T + T
o
0 + u.
J
CL + C.
lo
y =
+ y
-------
2-3
u* = [i* + u*
o
k = ko + k
(2.10)
If we expand the gas law (2.7) according to this scheme,
we have
= PoVo
and
P = RO(POT + PTO + PT)
The continuity equation (2.2) yields
= 0
which agrees with our assumption, and
(2.11)
(2.12)
(2.13)
+ s — P u. + pu.
dxj VKo j p j
The momentum equation (2.3) yields
dp.
= 0
and
p>
o
= - p,
du
(2.11)
(2.15)
PU
+ s—"- I +
dx.
+ u*)
m
m -J
(2.16)
This equation may be simplified further by noting that in the
first two terms on the left-hand side of the equation we may
-------
2-4
neglect p compared to p . Likewise on the right-hand side,
we may neglect y compared to y and y* compared to y* .
In addition, the term on the left-hand side in curly brackets is
zero by virtue of (2.14) so that we may write
9u
i
9u.
'o at + pouj 3x.
9x.
ax.
9u
m
9xm
(2.17)
Expansion of the energy equation (2.9) yields
and
9x.
J
9T
k
o 9x
J
P)
9T
/ v
= (°P0 ~ Ro j po 9t
p)u.
9T
= 0
(2.18)
+
-------
2-5
The two terms on the right-hand side containing p , namely,
+ u
3t
.
represent the heating of the fluid due to the motion and are of
p
order M compared to the other terms in the equation and so
may be neglected. Finally, by neglecting k compared with k
in the last term on the right-hand side, the energy equation
becomes
.
J
^ / -^ '"I r
a /, aT , , ° o
r k r + k r
o Y I O 'i v o Y
OX . I u OX. OX-
J \ J TF
(P0 + P)u-
a
"If * 5
J
[
ac
^0^0
a
ac
a
a
(2.22)
Since the basic atmosphere is all air. dC~ /dx. = 0 . Using
uo J
this result, (2.22) may be written to the same order of accuracy
as the other equations we have derived as
dC
dC
a _
ac -i
(2.23)
We now return to the equations that result from the perfect
gas law, namely, (2.11) and (2.12), and note that we may write
P_ = T_ p_
pT
P.T,
(2.24)
c O- Q ' O r O O
or to the order of accuracy of our other equations
PO To PO
(2.25)
-------
2-6
Now we note that the changes in p will be of the order of
2
J '
so that
T_
T
= o
u
R T
o o,
= OOyM)
(2.26)
We will, therefore, neglect the effects of motion-induced
pressure changes on the variation of density and take
= - T
- - T i
o
(2.27)
This latter assumption is related to, and entirely consistent
with, our neglect of the pressure work terms in the energy
equation.
We may now combine (2.14), (2.21), and (2.27) to find an
expression for the divergence of the velocity field. If (2.27)
is substituted into (2.14), we obtain, after some manipulation,
o j
o
(PnUi}
r
U
(2.28)
Neglecting T compared with T in the third term on the right-
hand side of this expression allows us to write
ou
1 • f 1
+
+ n T U.
PO o j
:j
. .
J
one finds that
(2.29)
,,_..T /Sx.) compared to
If the resulting equation is compared with (2.21),
In (2.29) we may neglect the term (T/T
O
1_
P,
k
.
J
'j J
(2.30)
This equation states that there must be a divergence of velocity
at a point where a volume of air is moving in the vertical
direction (the direction of the gradient of p ) so that air can
expand or contract to the pressure level of the equilibrium
-------
2-7
atmosphere. It also states that any atmospheric element gaining
or losing heat by conduction to the air immediately surrounding
it must expand or contract so as to maintain its pressure at
that of the equilibrium atmosphere.
It may be useful to gather together the equations we have
just discussed. For the atmosphere at rest, we have
p = R p T
o oro o
and
= Po°l
For the disturbed flow, we have
du.
p =r-— +
po St
C o ^
CPopo *t
(2.11)
(2.15)
(2.18)
(2.20)
Li / J
+ a
dx±
M-*
o
um
^.
oxm
(2.17)
k
(2.21)
and
J3
:J
P = - m-
a
R
o 1 B
ST
k
^XJ
+ k
(2.23)
(2.27)
(2.30)
-------
2-8
One recognizes immediately an inconsistency in the
equations for the basic atmosphere. Combination of (2.15) and
(2.20) yields
ST K
^ = - f±- (2.3D
Sxi CP0
which gives, if g. = g6.« ,
0 - ° = 0 (2.32)
and
<•! = _ ^L_ - constant (2.33)
3 °P0
which is the normal adiabatic lapse rate. Since by (2.18) we
must also have
k TT—— = constant
o dx^
the equations for the basic atmosphere are slightly inconsistent
since k is, in general, a function of T . This inconsist-
ency is such that the assumption of no time variation of T is
negated, since (2.18) states that, in general,
dT x / dT \
P ( CD ~ R } ^ = ?— k ^-2-\ (2.18)
^o V P0 o) St ^x, \ o dx, J
The inconsistency, however, is of small practical consequence
since (2.18) yields a time rate of change of T which is
enormously slower than the rates of change involved in practical
computations of atmospheric motion.
The equations given above are greatly simplified if one
considers that the region of the atmosphere in which one is
interested is of small extent compared to the scale of the atmos-
phere. In this case, the variations in the molecular transport
coefficients can, indeed, be neglected in the equations and, in
-------
2-9
addition, as will be shown below, the divergence of the velocity
field does not play a significant role in controlling the
motions that occur.
We begin with an order of magnitude analysis of (2.30).
Let us define a typical scale of the atmosphere as
-1
o
O
and a much smaller scale typical of the motion in question as
L
Equation (2.30) may then be written
k_AT kT
(L
_°_
L
0
0
1 - 0 ^- -2:
o
(2.35)
which may be written
1-0
p u .L
• r\ T c
AT ^a
T L
o s
In general we will have L AT/L T of order ten or less* and,
3. SO
since the Prandtl number u. Cp /k is of order unity and k/k.
is small, (2.36) may be written
chi. u.
1 - 0 I —•
(2.37)
Here the Reynolds number Re is defined as p u.L /u. . We
o j s o
see then that unless the scale in question is so small that Re
is of order ten the primary cause of divergence of the velocity
field is the basic density gradient in the undisturbed atmosphere,
and we may write (2.30) as
u
. .
_ J_ = __ 1
" -
0
(2.38)
* TQ/L^ is of order 1°/100 meters and, in general, we will be
interested in atmospheric flows with lapse rates up to the order
of 10V100 meters
-------
2-10
In the equations of motion, the velocity divergence enters
through the viscous term
du, \ Buk
(2.39)
The order of magnitude of the term on the left in the square
brackets is u. u./L while that on the right is of order
o j s
u. u./L . It is clear then that we may neglect the effect of
o j a
velocity divergence in the viscous terms of the momentum
equation.
It will be convenient to write terms in the equations of
motion such as p u.(du./dx.) and p u,(^T/^x.), etc., in the
ro j i j ro j
.
J
form
u
(QuJ - p Q
= -P,
c>x
J
Qu,
+ Qui^p_0
Pn ^X
(2.40)
Here again it is clear that the first term on the right-hand
side of (2.40) is of order u.Q/L while the second term on
the right-hand side is of order
write
u.Q/L
J
We may, therefore,
PO"J 5*7 ' PC 577 <2-
J J
in manipulating the equations . In view of the foregoing dis-
cussion, we may further simplify (2.38) to
c)u
=0
(2.42)
-------
2-11
In what follows, we will adopt the shear layer approximation
with the result that the coefficients of molecular diffusion can
be taken constant in the equations of motion and the velocity-
field is effectively divergence free. We will make one more
simplification, a simplification that is often made in analyses
of aerodynamic heat and mass transfer phenomena. In view of the
fact that the Prandtl number yC /k and the Schmidt number y/p.0"
for air are of order one (see Table 2.1 below), we will greatly
simplify the equations of motion if we assume that these numbers
are indeed equal to one, and write
ko =
and
Table 2.1 Prandtl and Schmidt Numbers for Several Gases at 0°C
and 1 Atmosphere
Gas. Pr = yC /k Sc = y/p.0*
Ne 0.66 0.73
A 0.67 0.75
N2 0.71 0.7^
CH* 0.74 0.70
02 0.72 0.74
C02 0.75 0.71
H2 0.71 0.73
Au 0.71 0.74
-------
3-1
3. EQUATIONS FOR TURBULENT SHEAR LAYERS
If the assumptions listed in the previous section are
adopted, the equations for an atmospheric shear layer may be
written:
PO = poRoTo (3'1}
(3.2)
(3.3)
(3.4)
(3.5)
(3.6)
(3.7)
9u.
r. i T P U .
o at o j
9T
P f\ j_ ~^~ P U .
O ot O J
90
n a 1 n n
P0 3t P0Uj
1
9p
o
ax.
9u.
i
9x.
j
9T
3xj
90
a
9XJ
- posi
9T
P n .9 — no-
CP0Po ax± ' posi
a P
P i O ~. rp i ,.
I g T + y
1 0 ,
92T
U° 9x2
J
92C
a
° 9x2
J
2
a"u.
9x2
J
o
and
9u.
alE*1 = ° (3>8)
j
A set of equations, generally referred to as the Boussinesq
equations for a stratified.fluid, is usually found in treatises
on atmospheric motion. That set reduces to (3.1) through (3.8)
when it is assumed that the Prandtl number and the Schmidt number
are both equal to one, as discussed in Section 2.
Following the technique used by Reynolds [Ref. 2], we now
develop the equations for the mean properties of a turbulent
atmospheric shear layer. To do this, we express the physical
variables as the sum of a mean value of the variable and an
-------
3-2
instantaneous fluctuation about that mean according to the
following scheme:
p = p + p'
T = T + T'
p = p + p' (3.9)
u. = 5. + uj
c = c" + c1
a a a
In these expressions, a bar over a quantity indicates the
average value of that quantity, while the prime denotes the
instantaneous fluctuation of the quantity about its mean.
The continuity equation (3.8) immediately yields the two
results:
du.
^-J- = 0 (3-10)
and
ou!
^ = ° (3.11)
(3.11) implies the further result
ui JT- = IT-^ (3.12)
tJ \J J\. -? **JJ\. -3
V (J
This relation is used frequently in the manipulations that follow,
If (3.9) is used in the momentum equation, the result is
/ aiL du|
Po ^ + XT^ ) + Pn ( Ui %T- + U." ^ + U1 ^ +
N N - . -
Bt ^t / Ko \ j 3xn- J Sx, J ox, j ox
J
x . ox .
J J
-------
3-3
If this equation is averaged, we obtain
2
9u. 9u. a- p _ 9u. 9u!u!
»o 31T + PO=J 55T- ' - £7 + T2 %T + "o 7TT - "o IxP (3
J 1 O o X J
This is the usual equation for the mean motion of an atmos-
pheric layer. It is seen that the turbulence has introduced
an effective stress, the Reynolds stress, which is of magnitude
pouiuj '
We may obtain an equation for the velocity fluctuations by
subtracting (3.14) from (3.13); thus,
9u!
o x i
Po 9t '
/- 3ui
n 1 11 . H
o J 9xi
\ J
9u.
hn- x
- 9p' 1
' ~ 9x± +
9u|
Uj 9Xj
p
o
9u! ^
u '
J 9XJ i
^L
in
•V° 3x2
J
From this equation we can obtain an equation for the Reynolds
stress correlation u-'u,*. . If (3-15) is multiplied through by
u
' and time averaged, one obtains
9u! / 9u! 9u. 9u!
p u,f -—• + p u.u.' v-i- + u'.u/ •*—=- + u'ul T—-
o k 9t Mo \ j k 9x. j k 9x. k j 9x ,
\ J J J /
J\n t K_
= - 11' ^ + — —
U +
92u!
11 n '
yu
k 9x. T ik ok , 2
1 O oX .
J
If the indices i and k in (3.16) are interchanged and the
resulting equation added to (3-16), there results, after some
rearrangement
-------
3-4
3u]u,'_ 9u_!u' 9u 9u.
K - - --tut £. _ n ii'iit —
Po W^- + poUJ 93^ = - Pouiuj 9T7 - poukuj
o 9x.
J
(P'u,
9u;
+
9u'\ p
—— I 4- — I e u'T' + e u'T'l
9v / T ( °i lr °k i J
OA^y J. \ J. K K.1 /
9x,
92u.'u; 9u! 9u'
+ u i k _ i k ( ,7}
uo » 2 ^yo 9x. 9x. U.in
3x, J J
This is the standard equation for the Reynolds stress correla-
tion and has been discussed at great length in various texts on
atmospheric motion [see, for example, Ref. 17].
We will briefly review the nature of this equation. The
two terms on the left-hand side represent p times the
substantive derivative (rate of change following the motion) of
the correlation u'uil • On the right-hand side, the first two
terms.
9u, 9u.
i ,1 •**• —. » t i ~^~
K,, ",•",• aY Kr> b- i ^Y
O1J dX. OKJ OX
J J
are called the production terms. They represent the production
of new correlations by the interaction of the turbulence with
the variation of the mean velocity field as well as the modifica-
tion of existing correlations by the variation of the mean field
(stream tube stretching effects).
The term
represents the diffusion of the correlation u'uil by the
turbulent velocity fluctuations and is called, appropriately,
the velocity diffusion term.
-------
3-5
The next two terms
j\ g
0 X . k O X, 1
i k
are called the pressure diffusion terms. They are the least
understood of all the terms in (3-17) because of the great
difficulty associated with experimentally measuring a pressure-
velocity correlation. In general, they have been found not to
play a dominant role in the development of a two-dimensional,
truly incompressible, turbulent shear layer.
The sixth expression on the right-hand side of (3-17)
is called the tendency-towards-isotropy term. This nomencla-
ture stems from the fact that, while the term appears in each
of the equations for the separate energy components in a truly
2
incompressible turbulent flow, i.e., in the equations for u] ,
22
u' , and u' , it drops out of the equation for the total
22 2
energy u' + u' + u' by virtue of the vanishing of the diver-
gence of the turbulent velocity field. The term must then, for
incompressible flows, represent a rearrangement of the turbulent
energy among the various components of velocity.
The seventh term on the right-hand side of (3.17), namely,
pn
=2. (g u'T' + g,u!T')
1Q i k k i
represents the production or reduction of the- correlation u'uw
due to the interaction of fluctuations in velocity with fluctua-
tions in temperature.
The term ? _
3 u!u'
_
in (3.17) is obviously the diffusion of the correlation u.'u'
1 ±C
by the action of molecular viscosity.
Finally, the term 9u! 3u,
'-).. _ 1_ _ _ l£
Mo 9x. 3x.
J J
-------
3-6
is called the dissipation term and represents, in the equations
222
of U-! , ui , and ul , the conversion of these energies to
heat by the action of molecular viscosity.
We will return to a discussion of the nature of these
terms later in this paper when we discuss how some of these
terms might be modeled.
If the Reynolds scheme (3-9) is applied to the energy
equation, the result is
^m Ci rn f \ / Cirn ^ i~P T Ci T
d ]. o L \ | — o 1 — o 1 , f a 1
po\9T 3t~ j po \ .j 3x. uj ¥xT uj 8x
j j j
When this equation is averaged, one obtains
.
J
Here again, the turbulence produces an increased heat transfer
whose magnitude is p CD u'.T' . The equation for the temperature
o "
fluctuation is
_ u,
U + U - U - y
o 3t o j Zx. j 9X j 9X - j 8x - o 3X2
(3.20)
We can obtain an equation for the heat transfer correlation
u'T" in the following way. First, multiply (3.20) by u^ and
time average to obtain
P u/ T- + P u,u' -- + u'.u/ - + u'.u'
ro k 8t o \ j k 9x. j k 3x. j k
-------
3-7
Next multiply the equation for 9u'/9t [i.e., (3.15)1 with i
replaced by k, by 1" and obtain
D T, 9uk + (
Po1 at po \
= _ T'
9u' 9u,
, , m i **• i , , t m i K- .
3 9x. j 3x.
J J J J
||L + ^gkT'2 + y '
+ u'.T'
J
a2 ,
uil
rn f ^_ „.._ _L _ Jy"
Xj
8uk
ax.
j
If (3.21) and (3.22) are added, one obtains, after some
rearrangement,
9T
PO St (uiT>> + pouj feT (ukT'> = - poujuk 577 - 'oujT' 33T
" po 3x ujuk 3xk p
3T
i 9. o- T11 4. n 9n ££ f3 o':5^
T~ SkX + yo , 2 ~ 2yo 3x. 3x. U.^3;
o 9x j j
J-
By analogy with the Reynolds stress equation, we may identify
the various terms in this equation for the heat transport corre-
lation. The two terms on the left-hand side represent P
times the derivative following a particle trajectory of the
correlation uilT' • The first two terms on the right-hand side
are production terms.' The first of these represents production
of u'1" by the interaction of the turbulent velocities and the
zC
variation of the mean temperature field, while the other term
represents.production due to the interaction of heat transport
correlations and the variation of the mean velocity field.
The third term on the right-hand side of (3.23) is the
diffusion of the correlation u'T' by the turbulent velocity.
The fourth term is the equivalent of the pressure diffu-.
sion term in the Reynolds stress correlation equation and will
be referred to as the pressure diffusion term. Similarly, the
-------
3-8
term p'dT'/dx represents a "tendency towards isotropy" or a
tendency for the pressure fluctuations to eliminate any corre-
lation between u' and T' .
The term
T^gkT'
represents the production of temperature transport correlation
in the direction of gravitational acceleration due to the
fluctuations in temperature of the flow.
The term
^ -^-
is obviously the diffusion of u'T' by molecular action, while
the final term
" o dx, ox,
j j
represents the role played by molecular viscosity in correlating
or uncorrelating the velocity and temperature fluctuations. In
most cases, the effect is.one of uncorrelation. In spite of the
fact that this last term is not always dissipative, it is none-
theless referred to as the dissipative term.
2
In (3.23), the correlation T' occurs. It is easy to
obtain an equation for this second-order correlation. To do so,
we multiply (3.20) by 2T' and average. The result is
2
o
This equation is similar to the two previous equations for
second-order correlations with the exception that there are no
pressure diffusion or tendency-to-isotropy terms. It is clear
-------
3-9
then that the only way that temperature fluctuations can be
destroyed in a flow of constant temperature is by the action
of molecular diffusivity.
When the Reynolds scheme is applied to the species
conservation equation, one obtains
a a a
a a a
-
p + 5 + - , ,
po I dt dt / Ho I j dx. J dx j dx. J dx,
* / \ J J J J
d2C!
"
o2C
J
Averaging the equation results in
a . - ~~a a
(3.26)
72
As before, we can obtain equations for C'u,' and C' . They
Gu Jv C6
are
' n '
po at— + pouj axT (ukca) = - POUJU
p1
m2 g, C'T' + u 0K a - 2u. ^ ^ (3.27)
TQ k a o Sx2 o OXJ dxj
and
PO 3^ + PoSJ %(Ca2)=-2Po«jC^
dc
~(uic'2) + ^n —^ - ^o ^T-9 (3.28)
X. V j a y ^o 2 , ^0 ^X £X
-------
3-10
It will be seen that because of the gravitational term
(p /T )g.C'T' in (3.2?) it will be necessary to have , an
O O it Ot _______
equation for C'T' if a closed set of equations is to be
realized. This equation is derived in the usual manner and is
p0 It '^ + POGJ - rr =
J
oir
n
= - n n'P' _ n n "T" _ — - n - _ ( n ' f ' T '
poj a 9x. pouj 9x. po 3x . Uj a
J J J
9"C'T' 30' 3T
"o -T- - 2lJo 93T 35
We need not discuss these last three equations in detail
for they are almost identical in form to the equations for
u~VTT and T'2 .
K
Let us now collect the equations we have derived for the
production of turbulence and the dispersal of pollutants in an
atmospheric shear layer.
Making use of the notation
arid writing y /p = u , we have
9-
Du. -13- a u. ~ .,
K^ = - ±- |f- + u —^ - |— uTIJ + i- g T (3.31)
Dt PQ ^x± o 9x2 3Xj i J To i
DT _ 9_T_ _ 9— ^TTT ,_ _ s
Dt o 9 2 3x, j v^o^J
•-)
rir1 9 P
Dt^^or^-fcr11?! <3.33)
o A . J
-------
3-11
and
dx,
= 0
(3.34)
Duiuk
Dt
du.
u'u' 's-
j k dx.
j
u'.u'
..
J i dx
^— (unXu^
dxn. i J k
(P'uj
(P'u') +
+ v
duj
'o dx, dx
4
J
(3.35)
K.
' 11 '
Dt
- u'.T'
J
. ,
dx, J k
1_
PO
TfT) + El ^1
+ V
- 2v
dT-
(3.36)
DT'2
Dt
f rp |
... -------- .
dX j
u!T'2
v°5?
- 2V
dT' dT'
o dx. dxj
(3.37)
-------
3-12
Du'C' 3C _ 3u
_ ,
a _ rrr '
_ nr
" ~ Uu
Dt " ~ k 3x " ja 3x. 3x
jk 3x " ja 3x. 3x
J J J
+ P'
P0 3xk P0
9u; sc1
(3'38)
DP ' T ' —
_jx__ _ —j—r 3T
Dt = ~ UJ g 3x7 "j" 3x. ~ 3x.
J J J
im t
3 C 'T1 9C ' 31"
vo-2*-- 2vo inr 3F- (3>39)
9 X U O A .
-------
3-13
First, let us assume that the motion in question is laminar.
In that case, all the fluctuations are zero and the equations
needed to solve for the motion alone, i.e., for the five unknowns
u1, u2 , u- , T , and p are (3.3D, (3-32), and (3-3*0. Since
(3.3D is a vector equation, it yields 3 scalar equations for u,
Up, and u_ so that, together with (3.32) and (3-3*0, there are
five equations for the five unknowns. Once the motion is known,
the mass fraction (concentration) field for each species a may
be found from the a equations given by (3-33).
If the motion is turbulent and a viable second-order closure
model can be found, the unknowns involved in the solution for the
motion are u u, u, T, p, u-u-,
_
-
'2
uj = u'u| , u£u' = u'u£ , u|T' , u£T" , u'T1 , and T' . There
are then 15 unknowns. One obtains the 15 equations required to
solve for the motion as follows:
from (3.31) 3 equations
from (3-32) 1 equation
from (3.3*1) 1 equation
from (3-35) 6 equations
from (3.36) 3 equations
from (3.37) 1 equation
for a total of 15 equations.
Once the motion is solved for the concentration field for
any of the a , each species involved may be solved by solving
five simultaneous equations for the unknowns C , C'u-! ,
01 Ot J-
CjuJ , CjiU , and CjT7" obtained
from (3-33) 1 equation
from (3.38) 3 equations
and from (3.39) 1 equation
Once this solution has been obtained, the variance of the
fluctuations of mass fraction of the a species may be obtained
from a solution of (3.*IO).
-------
4-1
4. SELECTION OF MODELS
In order to construct models of the four types of terms
referred to in Section 3 so as to form a closed set of equations
for solving for the motion of the atmosphere and its ability to
disperse pollutants, one can be guided by some very general
principles which reduce the number of models that might be
investigated. Once models that adhere to these general principles
have been selected, resort must be had to a comparison of computed
results with experimental data in order to refine the models or
to determine the numerical values of certain parameters which
occur in the construction of the models. At the present time,
there are more models for closure of the equations of motion at
the second-order correlations than there are principal investiga-
tors working on the problem. This occurs because each investiga-
tor generally has an option or two that he is currently trying
out. In our own work, we have tried.to start from the simplest
model that would satisfy the general principles and increase the
complexity of the model only if experimental data indicate that
more complexity is needed.
The basic principles that must be observed in selecting a
model are the following:
(1) The model must be of tensor form so that it is invar-
iant under an arbitrary transformation of coordinate systems.
The model must have all the tensor properties and, in addition,
all the symmetries of the term which it replaces.
(2) The model must be invariant under a Galilean transform-
ation.
(3) The model must have the dimensional properties of the
term it replaces.
(4) The model must be such as to satisfy all the conserva-
tion relationships known to govern the variables in question.
With these basic principles, let us try to determine the
simplest models possible for the various terms that must be
modeled in (3.35) through (3-40).
-------
4-2
We will start with the velocity diffusion terms. We wish
to model the tensor u.'ulu/ in terms of the second-order
1 J K
correlations u.'u,1 . The simplest tensor of rank three that
1 K
can be obtained from the second-order correlation uJuil that
is of the form T. ., is d (u! u') /d x . . The tensor u'u'u' is
symmetric in all three indices so that our model must be
symmetric in these indices. We therefore choose
S
u'ulu' ~ i — (u'.u/) + 5 — (uju.1) + £ — (u'u«)
i j k dx- J k dX 1 k' dx i J
This expression has all the tensor character necessary for the
model but is riot dimensionally correct. To make the dimensions
correct, we must multiply the right-hand side of (4.1) by a
scalar with dimensions of length times velocity. The simplest
scalar velocity we can form from the second-order velocity
correlations is
q = / u'u'
mm
If we multiply this expression by a scalar length, say, A~ >
which is to be related to the scale of the mean motion or to
the scale of the turbulence, we can form a simple model of the
triple correlation uju'.u' . Thus,
:uiuj>]
(4.2)
In (4.2) we have included a minus sign to ensure that turbulent
energy will be diffused down the gradient or from regions of
high turbulent intensity to regions of low intensity.
Let us consider a model for the triple correlation u'u'T'.
J K
From the second-order correlation uilT' we can form a tensor of
form T., with the necessary symmetry. Thus,
J •"£
m _ O / t m | \ ^. O ( U ' T ' ) (4 ^ )
Jk dx, k Sxv j
-------
4-3
This expression does not have the dimensions of u'.u'T', so
J K
we write
u'u'T' = - A q
j
•**-
(u'TM + ^— (u'T')
Here we might have chosen a scale length other than A2 . But,
as we shall see later, in order to keep the model simple it
will be necessary to choose as few scale lengths as possible.
We therefore take A? as the scale for all velocity diffusion
terms unless experimental results later show that the model must
be more complicated.
In modeling the term u'u'.T1 we might have chosen a term
of the following form:
const u.'u! /T1'
-L J
We have chosen not to do this since we want this diffusion term
to be a gradient type of diffusion.
By analogy with the velocity diffusion expressions (4.2)
and (4.4), we will choose to model the four remaining triple
correlations that appear in the velocity diffusion terms as
u'.T'2 = - A0q |— T'2 (4.6)
J 2 dX
u'.C'T' = - A0q ^— C'T' (4.7)
J a <* dx« a
Consider now the tendency-towards-isotropy term in the
equation for u-uj!- • We wish to model
-------
4-4
po \9xk
This expression is a second-order tensor of form T. . The
1K.
simplest tensor of this form that can be formed from the second-
order correlations is obviously the second-order correlation
u.!u' itself. We therefore start by taking
9uk'
- u,
The correlation u'uil has all the tensor and symmetry properties
we wish, but it does not have the proper dimensions. We can.
remedy this by multiplying u.'u' by the scalar q/A, where A,
is a new scale length. We then have
This cannot be correct, for if we set i equal to k , the left-
hand side of the supposed equation vanishes in an incompressible
fluid while the right-hand side does not. This can be remedied
by -taking
(4.10)
The minus sign in this expression is chosen so that in the
absence of other influences, the turbulent energy will be equil-
ibrated between the various components of velocity. The model
was first given by Rotta [Ref. 6]. We could have included a
term in (4.10) proportional to the mean strain [Crow, Ref. 43],
but for our initial simple model, we choose not to do so.
We must now model'the "tendency-towards-isotropy" terms
that appear in the equations for u'T' and u'C' . By analogy
with the tendency towards isotropy of the stress correlation for
i ? k , we adopt the following models:
-------
4-5
po
and
, 3C1
Po 3xk lka '
Here again, more generality can be put into the model by selecting
the scales in these last two equations to be different than the
A, of the isotropy term for the stress correlation equation.
Whether this should be done or not will depend upon the agreement
between computed results and experimental data. For the time
being, we will use the same scale in each of the isotropy terms.
As mentioned previously, the pressure diffusion terms are
difficult to model since so little is known concerning them.* We
will, however, choose a model of the terms which will be essen-
tially nonproductive; i.e., the model will depend upon the grad-
ient of the quantity in the equation for which we are trying to
model the pressure diffusion term. Thus, we choose
1
3
(4.14)
and p-c; = - pQqA3 J^ujc; (4.15)
In these expressions, since we. expect the pressure diffusion
effect to be small, we expect that A_ will be small compared,
say, to A . We have chosen a negative sign so that the model,
as written, will represent a transfer of the quantity in question
down the gradient.t
*0ur knowledge of these terms is limited for, to date, it has
not been possible to simultaneously measure the fluctuation of
pressure and velocity at a point in a flowing medium.
tit should be noted that in some recent studies using a model
which is essentially that described here, Prof. Paul Libby has
reported somewhat better agreement of computations with experi-
mental results if a small Ao is used and the signs of the
right-hand sides of (4.13), ^4.14), and (4.15) are reversed.
-------
4-6
To model the correlations appearing in the dissipation
terms, namely,
9ui
9Xj
30'
a
9x.
J
9uk
9x. '
J
9T'
9Xj ''
9u' 9T'
1 d -v ^ v '
O A . 0 A .
J J
9C'
nnrt a
, ana a
8XJ
^\ rn | f\ m f '
O 1 O 1 (
' 9x. 9x. J !
J J
9Ca
3X.
}u,' 9C'
k a
)x. 9x.
J J
it would appear that the simplest model that can account for
all of these terms is
9A' 9B' A'B' ,, M
9x. 9x. ,2 . ^•lt)^
J J X
where X is another scalar length.
Here again we have the same question that was raised
before; i.e., should not. the X's in each term be different.?
The answer to this must again be "yes" if comparison between
computed results and experimental data requires it. As .in the
other cases, we will, for the present, try to develop the
simplest possible model and will therefore make use of only one
X in all the dissipation terms. In later sections of this
report, we will return to the question of•modeling these dissi-
pation terms in light of the results found with the model
adopted here.
If the models we have described above are placed into the
equations for a turbulent atmospheric shear layer, (3.31)
through (3.40), we obtain the following set of modeled equations
Du. 9 u. „ ., -,3-
_i. = v ±-$ iTuT^ + i- K T - — &- (417)
Dt o,2 9x. i j T si p 9x. ^'x';
9x. j J o Ko i
§f - -o rf - k: ^ <4-18)
-------
3x
u!u,' +
A.
ujui
Du^.T'
Dt
= - ulu
«111
pirn ""V ~\
°J: u t rn , K + J^_ rp ,
j"k 3x, V 3x, T V
4-7
DCa v !!fa
Dt - Vo 9x2
|— u'.C'
3x. j a
(4.19)
3u.
= 0
(4.20)
Duiuk
Dt
3u.
3 a,
• - ujui ixy - ujui
+ V
3X2
J
uiuk
(4.21)
3 Yl
£ ,, irn | I
3x. UkX
J /J
A.
3x.
u'.T
**
AkJ
9x2
J
UkT'
(4.22)
,2
Dt
.
J J
T
,2
0
-------
4-8
Du'C' 3C 3u,
k a = -n—r __a TTPT _j£
Dt ~ ujuk 3x. j a 3x.
J J
~ A, ukCd
DC'T' .FfT 3C
g = u ' C ' _ u ' T ' —-
Dt i a 9x. 1 9x.
d J J J
9x
J
DC'2
a
Dt
u!(
J
92(
8x'
own a / a
-. i a J (An
"a 9x. 9x. I 2q 3x.
J J \ J
P P
:,^ c,^
a _ a
? ° X2
C'2
a
+ v^ -^- - 2v -4- (4.26)
In the pressure diffusion terms, p has been taken out of the
derivative by the same reasoning that led to (2.4l) from (2.40).
This set of modeled equations is solved according to the proced-
ure described at the end of Section 3-
For thin shear layers, these equations can be further
simplified by making"boundary layer" assumptions. First we
introduce the notation
x.^ = x x2 = y x_ = z
u. = u Up = v u_ =-w (4.27)
-------
4-9
Then we assume that the vertical component of the mean velocity
is small compared to the other components
w « u or v (4.28)
and that derivatives with respect to z are large
H»H • If <«•*»
The equations are of particularly simple form if, for
example, the problem is time dependent and all mean variables
(excepting p) are independent of x and y ; that is, if
u = u(z,t), T = T(z,t), u'w' = u'w'(Zjt), etc. In this case we
must have (1/p )(3p/3x) and (1/p )Op/3y) (which are
functions of z alone) specified as functions of time, and we
find that the vertical momentum.equation becomes
H = - Pg (4-30)
while continuity gives w = 0.
For this problem, since the pollutants involved have been
assumed to be passive, we have a coupled set of thirteen
partial differential equations for the thirteen variables:
u(z,t) , v(z,t) , T(z,t) , u'u.'(z,t) , v'v'(z,t)
w'w'(z,t) , u'w'(z,t) , v'w'(z,t) , u'v'(z,t)
u.'T'(z,t) } v«Tt(z,t) , w'T'(z,t) , T'2(z,t) .
If initial conditions on these quantities as functions of z
are given at some time, say, t = 0, the equations can be solved
for the development in time of the mean quantities and the
correlations. Once the turbulence characteristics of the atmos-
phere are known as functions of z and t , one is then able to
solve for the transport of the pollutants by solving-a coupled
set of five equations for C (z,t) , C'uT(z,t) , C'v'(z,t) ,
LX LX UC
C'w~r(z,t) , and C'T'(z,t) . Once these quantities are known,
the variance of the fluctuation in C can be determined by
-------
4-10
2~
solving the final equation for C' (z,t). These solutions can
be obtained, of course, providing that initial conditions on
these quantities as functions of z are provided at some time.
Of particular interest is the steady dispersal of passive
pollutants in an atmosphere that has reached a given steady
state of turbulence. Since the equations for the atmospheric
motion and for the dispersal of pollutants are decoupled, one
can first solve for the atmospheric turbulence that would be
generated by known profiles of velocity and temperature, say,
u=u(z) , v = 0 , w = 0 , and T = T(z), and then solve for
the steady-dispersal of pollutant material in this atmosphere.
The turbulence 'characteristics of the atmosphere used for these
latter computations are those found by holding u(z) and T(z)
fixed in the equations for the second-order correlations u'uil >
ujT7" , and T'2 (these equations are (4.21), (4.22), and (4.23),
respectively) and observing the resulting distributions at large
times. These distributions, as we shall see, are independent of
the initial conditions on these variables used to state the
calculations at t = 0 .
It must be borne in mind that in order to perform the
calculations just mentioned, it is necessary to have equations
for the scales A, , A? , A_ , and A or have these scales
given as functions of z and t .
Before we can proceed further, it is necessary to make more
precise the details of the model of turbulent atmospheric flow
that has just been presented. To do this, we must construct
models for some well-known flows and compare the results with
experimental data. In this way we will determine relationships
between the various scales A and X , by choosing these
parameters so as to give the best results for the widest number
of experimental observations and determine in this way the
property of these flows which is related to these scales.
To make the calculations necessary to carry out the compar-
ison between experimental results and calculations suggested
-------
above, it is necessary to have computer programs which accurately
solve the sets of equations we have discussed. A general tem-
plate for the solution of such systems of equations has been
developed at A.R.A.P. over the past few years. A description of
these numerical techniques is contained in Appendix A.
-------
5-1
5. SEARCH FOR MODEL PARAMETERS
It is easy to verify that equations (4.17) through (4.26)
reduce to the equations for a low Mach number, constant density
medium if we assume T E 03 T' E 0, C' E 0. They may there-
fore be used to calculate the generation of turbulence in a
classical incompressible shear layer for which a great deal of
detailed experimental data are available. The detailed model
of turbulent flow that we will use in later sections of this
report was developed by comparing computations with experi-
mental results for three shear flows :
(1) the axially symmetric free jet;
(2) the two-dimensional free shear layer; and
(3) the flat plate boundary layer.
To calculate the first of these flows in cylindrical
coordinates (r,<|>,z) with velocities (u,v,w), one finds that
in the appropriate boundary layer form the equations are
- 3u'u' , - 9u'u' 9
U 9r- + W 9^ = 9?
(3A2 + 2A3)q
24 /
r I
2 9"
r 9r
0 9uT
J 9r
A3q(
+ 2A3
u1
u'u'
)q
(\
„ 9v'v
<- <\
9r
- v'v'
1 1 ,1 1
t \
1
>]
^r I ir 1
3u'u.'
(5.3)
-------
5-2
u
-dv'v' a /, dv'v'\ . 3A2q dv'v1
w 3 = K— A Oq 3
az dr \ 2^ dr
A0q(u'u' - v'v1
2A0q
(4A + 2A?)q
+ £ —I— (u'u1 - v'v') - ^-l v'v
2v'v
X2
, i a w' w' , i_ aw
2w
X
(5.5)
u
' w
, -
+ w
' w
2A0q
_
r dr
+ -%- (A
r dr
(2A2 + A,)q
- p—^— u'w' - ^— u'w'
r 1
.d'u'w' !_ du'w' u'w' 2u'w'
(5.6)
-------
5-3
For the second and third flows, if u is taken as the
free stream velocity in the x direction and z the direction
normal to x in which large gradients in u exist, the
appropriate equations are
+ dw
d~x~ d"z~
(5.7)
- du , - du , d u du' w'
U 5— + W s— = V
ox dz
dz2 dz
(5.8)
u du'u' + - Su'u' =
dx dz
- 2u'w
''
du'u
''
2u'u
(5.9)
u
dx
dz
T~
dz
d v'v1 2v'v'
,dz'
X
(5.10)
- dw'w' , - dw'w' d I/,,, , OA
u — + w — = ^ 2
\ dw ' w ' I
3 ^ — J
''
w' 2w'w'
dz
(5.H)
u
dx
dz
= - w
'w' ^ + ^-
dz dz
A3),
u ' w ' + V
,. / d u'w' 2u'w'
z2 " X2
(5.12)
-------
5-4
It is a fact established by careful experimental observa-
tion that both the free jet and the two-dimensional free shear
layer become, after a transient due to initial conditions, self-
similar turbulent flows at high Reynolds numbers. If (5.1)
through (5.6) and (5.7) through (5.12) are to permit self-
similar solutions for the free jet and the free shear layer,
respectively, it is necessary that certain relationships exist
between the scale length A , A , and the characteristic•scale
of the mean motion 6 , .It can be shown that these rela-
char
tions are
A = c o
"1 \ r* n Q T°
_1_ -L w lictX
A2 = c^ = c£
A0 = cnAn = c'6 , (5-15)
3 31 3 char
and
A = A^/a + bReA (5-16)
where Re. = pqA-,/y (5.17)
1
The relation (5-16) between the dissipative scale A and the
isotropy scale A, is required so 'that as the Reynolds number
of the flow increases, dissipation will keep pace with the other
productive and diffusive terms in the equation so that a self-
similar flow can result. This form of the relationship between
A and A, has been used in the past by Glushko [Ref. 7].
For self-similar free turbulent flows, the structure given
above is all that is needed to compute a turbulent shear layer'
or a free jet, provided the five constants, c., , c? , c_ , a ,
and b , are given. To find these constants, we must resort to
the comparison of calculated flow fields with experimental
results.
If we wish to compute a boundary layer flow, we must
consider an additional problem. When a wall is present in a
shear flow, we wish to apply the boundary condition at the wall
that
-------
5-5
• 0, we'must have
X = az (5.19)
2
where a = 2/(l + n)n . Thus, near a solid surface, we always
assume, in applying our model, that (5.19) holds in the region
near the wall.
It is convenient to express this result in terms of A, .
Near a wall, (5.16) becomes
X = A^/a (5.20)
Using (5.19) we may write
A1 = a/az (5.21)
-------
5-6
Thus, for boundary layer flows, a Is another number which must
be found from experimental results.
In our first attempts to construct a model of turbulent
shear flows [Refs. 10 and 16], the following assumptions were
made in order to construct the simplest possible model of
boundary layer flows:
(1) It was assumed that all the- large lambdas.associated
with inviscid modeling were equal; i.e., A, = A? = A_ = A.
(2) It was assumed that a was equal to one.
(3) In the outer portion of a boundary layer, A was
taken to be a constant c-, times 5 nn (& nn is the value
_ -L «yy • yy
of z for which u is 99% of the free stream velocity). This
value was assumed to hold, independent of z , as the wall was
approached, until A became equal to /a times z . For
smaller values of z, A was taken equal to /a z.
With these assumptions, (5.-7) through (5.12) were solved
with various choices for the parameters a, b, and G-, = A/6 QQ
to produce a developing turbulent boundary layer on a flat plate.
It was determined that the following choice of parameters
cl = A/6 99 = °'°6^
a = 2.5 (5.22)
b = 0.125
yielded a fair representation of a turbulent boundary layer. The
mean velocity profile and the.behavior of skin friction with
Reynolds number were adequately represented. The distributions
of the second-order correlations within the boundary.layer were
reasonable.
The results of this original parameter search were used to
compute a number of other turbulent flows in order to demonstrate
the method [Refs. 19 and 20].
Before proceeding with-further applications, it was
considered necessary that ^a more detailed parameter search be
made. In particular, two free turbulent flows - the free jet
and the free shear layer - should be calculated to determine
the values of the parameters Cp, c_, a and b that would best
-------
5-7
fit the experimental results for both flows. The value of c,
being the ratio of A, to some arbitrarily defined character-
istic length in each case is not an invariant of the problem
and was to be chosen with fixed values of the other parameters
to obtain best results in each case. Once these studies were
complete, the model would be used to compute turbulent boundary
layer flows so that, by comparison with experimental results,
values for G]_ and a could be made for this flow. Hopefully,
all flows could be described in a reasonable way by a single
choice of the basic model parameters c?, c^, a, and b, and
(where appropriate) a . The values of local A, determined
from the values of c, in each case were then to be compared
with the local magnitude of the integral scale L in each case.
If it was found that the value of c-, represented a choice that
amounted to
A = const L = 3L (5.23)
then it would be assumed that a reasonably invariant model had
been determined.
Application of Model to Free Shear Flows
Our search for a new model of turbulent shear layers began
with an attempt to describe the axially symmetric free jet with
the original turbulence model obtained for a boundary layer flow.
This model, as mentioned previously, was one for which A = A =
A-, = A . This choice leaves three parameters to be determined.
They are c, = A/6 , and the two constants a and b in the
expression•
X = A// a + b • ReA
The method of searching for values of these parameter was as
follows. The equations for a free .jet were programmed so as to
solve the system of equations for a free jet developing in the
axial direction. At an arbitrary initial station in the axial
direction, a mean velocity profile and profiles of the pertinent
second-order correlations were arbitrarily assumed. For a given
choice of model parameters (in this case, a, b, and c, = A/r
-------
5-8
where r c is the radius for which u is one-half the center-
. b
line value), the free jet equations were solved for the develop-
ment of the jet downstream of the initial distributions. In all
cases, essentially self-similar solutions were obtained far down-
stream of the start of the calculation. If a.set of parameters
could be found so that the resulting self-similar flow agreed
with experimental measurements with respect to the rate of
spread, as well as with respect to mean velocity and correlation
distributions, it would then be assumed that a reasonable•turbu-
lence model had been achieved.
Actually such calculations were carried out for both free
jets and two-dimensional free shear layers. With the single A
model, it was found that ^no combinations of parameters a, b,
and c, could produce an adequate description of either a free
jet or a free shear layer. In general, it was found that if the
parameters were adjusted so as to give an adequate rate of
spread of the mean profile (i.e., if the level of the turbulent
shear correlation was large enough) the spread of the correla-
tions u'uil by diffusion was always too large. This general
result is illustrated in Figure 5-1 where it is seen that, if the
general level of the shear correlation u'w1 were to match the
experimental data of Wygnanski and Fiedler [Ref. 21] in the region
of maximum shear, it is clear that far too long a tail of u'w'
at large r would result. This was a very general result for
free shear flows and forces us to consider a more complicated
model.
The difficulty.that was experienced with the constant A
model was the existence of too much diffusion relative to the
rate of loss of correlations, either by dissipation or the
tendency towards isotropy. To correct this difficulty in the
studies reported here, the diffusion lengths A~ and A_ were
made smaller than A.. . An idea of the effect of reducing the
diffusion lengths relative to the isotropy length can be seen
from Figure 5.2. Here the rms value of the longitudinal velo-
city fluctuation w1 that has been calculated for several
choices of model parameters is plotted versus radius in a self-
-------
5-9
.020
.016
U'W7
W2
wmox
.012
.008
.004
Wygnanskl 8t Fiedler
\ a = 4
b=0.005
A,=0.3 r>5
\
\
r/r
.5
Figure 5-1. Result of a free jet computation with
a single A model of turbulence
-------
5-10
.32
.28
.24
.20
yw'w'
w2,ax .16
.12
.08
.04
Wygnanski
Fiedler
r/r
.5
Figure 5.2. Behavior of solutions for a self-similar free
jet when the parameter c? =
is varied
-------
5-11
similar free jet. Note that as the diffusion lengths A~ and
AO (which are c~ times A ) are reduced, the amount of diffusion
is obviously reduced and the levels of turbulence on the jet
centerline are appreciably increased.
The effect of the choice of the scale of the isotropy
length A, can be seen from Figure 5.3. The distribution of
longitudinal turbulence intensity is shown as a function of
radius for two choices of A_ relative to the local value of
r ,. . It is seen that the levels are much lower for the smaller
• o
A-, than for the larger value. This is what one might expect
because of the increased dissipation as well as the increased
loss of shear correlation by the tendency towards isotropy when
the scale A-, and, hence, \ is made smaller.
The effect of neglecting pressure diffusion can be seen in
Figure 5.4; the longitudinal velocity fluctuations in a free jet
are shown as a function of radial position for a given choice of
model parameters a , b , c, , and c~ for two choices of c~ .
One choice is c~ = Cp and the second is c- = 0 , i.e., neglect
of pressure diffusion. It is seen that for this choice of the
other parameters the effect of neglecting pressure diffusion is
not large.
Having given some idea of how some of the various para-
meters entering the model for turbulent shear layers affect the
solutions, we must now discuss the selection of an actual set of
parameters. If one considers only a single type of shear flow
that one wishes to model, say, the free jet, it is possible to
choose a whole spectrum of models which will give a good
description of the mean spread of the free Jet and the distribu-
tion of, say, the longitudinal turbulent velocity field. To
illustrate this point, we may refer to Figure 5-5. Here we see
that two profiles of longitudinal velocity fluctuation can be
obtained with radically different choices of b and A-, . It
is observed that if one chooses small b one must also choose
a small value of A relative to a characteristic scale of the
jet. What then is the basic difference between these two solu-
tions? It is this. For the solution with small b and small
-------
5-12
wV
w2
wmax
28
.24
.20
.16
.12
.08
.04
Wygnanski 8t Fiedler
a=2.5
b=0.05
=0.2
=0.2
Figure 5.3. Effect of variation of the isotropy scale A,
on characteristics of a self-similar free
jet
-------
5-13
Iw'w'
w*
wmax
.28
.24
.20
.16
.12
.08
.04
Wygnanski 8 Fiedler
r/r
.5
Figure 5.4.
The effect of neglecting pressure diffusion
when calculating a self-similar free jet
-------
.28
.24
.20
.16
.12
.08
.04
Wygnanski 8 Fiedler
a=2.5 b=0.05
A,=0.02r5 c3=0
C2=0.l
a=2.5 b = O.I25
A,=0.5r5 C3=0.l
'".6
Figure 5-5. Two choices of model parameters that yield
almost the same distributions of w'w' for
a self-similar free jet
-------
5-15
A, , the balance of the production of turbulence is more by
dissipation and less by diffusion than for the other case. Also
for the case of small b and small A-, , the solutions are more
isotropic on the jet centerline than for the other case.
The choice between the two models exhibited in Figure 5.5
must be made on the basis of the degree of diffusion and the
degree of isotropy desired in the calculated result. This is a
difficult decision to make, for existing experimental data do not
agree as to how isotropic free jets are on their centerline, as
will be seen later. There is another way that one can decide
between two different models. If one uses the same model to
compute two different turbulent flows having essentially differ-
ent geometries, the model which gives the best results for both
flows is, since we are seeking an invariant model, the one to
choose .
As mentioned previously, we have computed self-similar
solutions for a free shear layer as well as for an axially
symmetric free jet. Actually a search for model parameters for
each type of flow was carried out. As a result of these studies,
it was determined that, insofar as the parameter studies have
proceeded at this point, the following model for free turbulent
shear flows gave the best results:
a = 2.5
b = 0.125
-0.1 (5
Also, the value
2
= 0.1
l = V6char - °'50 (5'25)
was found best for both flows, although it was not part of the
plan to have a common value of c, . As mentioned above, for
the free jet,
6char = r.5 (5'26)
The characteristic length for the free shear layer was taken as
6char = z.25 ' Z.75 (5'2?)
-------
5-16
which is the distance normal to the flow in the shear layer from
the point where the velocity is one-quarter the external driving
velocity to where it is three-quarters this velocity.
In Figures 5-6 through 5.13, we show comparisons with
experimental data of the velocity correlation profiles computed
for both a free jet and a free shear layer using the model para-
meters given above. The experimental results are taken from the
work of Wygnanski and Fiedler [Refs. 21 and 22], Gibson [231,
and Donaldson, Snedeker, and Margolis [24].
Figures 5-6 and 5.7 show the longitudinal fluctuations in
a free jet and free shear layer, respectively. The agreement
between model calculations arid experiment is good in both cases.
For the free jet in Figure 5.6, it would perhaps have been
desirable to have a little more diffusion (larger A-, and larger
b) in the model in an attempt to reduce the overshoot in w'w'
near the centerline of the jet.
Figures 5-8 and 5-9 show distributions of normal fluctua-
tions in both the free jet and the free shear layer. Here we note
the agreement with experimental data is not so good. There
appears to be a little too much diffusion for these cases. Also
note the very large discrepancy between measured normal fluctua-
tions on the centerline, as reported in three separate experiments.
The data of Gibson show the components of turbulent velocity to be
essentially isotropic on the jet centerline, while those of
Wygnanski and Fiedler and Donaldson, Snedeker, and Margolis do
not. From the results shown in Figure 5.8, it would appear that
if one were to desire more isotropy, one would wish to choose a
smaller value of A-, and, hence, a smaller value of b . This is
opposite to the conclusion drawn from Figure 5.6.
Figures 5-10 and 5.11 show the sidewise components of turbu-
lence for the free jet and free shear layer, respectively. The
agreement between experiment and computed results is better for
the free jet than for the free shear layer. The reason for this
behavior is not known.
In Figures 5.12 and 5.13, we show the shear correlations for
the free jet and the free shear layer. The agreement in both cases
-------
5-17
o Wygnanski 8 Fiedler
Gibson
A Donaldson, Snedeker
8 Margolis
Figure 5.6. Comparison of experimental results and model
predictions for the longitudinal velocity
correlations in a free jet
-------
5-18
uV
u
max
a=2.5
c2=0.l
b = O.I25 c3=0.l
A,a0.50(z25-z78)
o Wygnanski 8 Fiedler
*25~z.75
Figure 5.7
Comparison of experimental results and model predictions
for the longitudinal velocity correlations in a free
shear layer
-------
5-19
.24
.20
.16
.12
.08
.04
o
I- o
0 = 2.5
b =0.125
A,=0.50r«
I .9
c2=0.l c3 = 0.l
a Gibson
A Donaldson, Snedeker
8 Margolis
o Wygnanski 8 Fiedler
r/r
.5
Figure 5-8. Comparison of experimental results and model
predictions for the radial velocity fluctuations
in a free jet
-------
5-20
w'w'
.05
.04
.03
u
max
.02 -
.01 -
a = 2.5 C2 = 0.l
b =0.125 c = 0.l
= 0.50(z25-zT5)
o Wygnanski 8 Fiedler
o o
z.25"z.75
Figure 5-9
Comparison of experimental results and model predictions
for the normal velocity fluctuations in a free shear layer
-------
5-21
v'v'
w2
mox
a = 2.5
b=O.I25
A, = 0.50 r 5
c2=0.l C3 = 0.l
o Wygnanski ft Fiedler
a Gibson
A Donaldson, Snedeker
ft Margolis
r/r
.5
Figure 5-10. Comparison of experimental results and model
predictions for the sidewise fluctuations in
a free jet
-------
5-22
vV
"max
.05
.04
.03
.02
.01
-2
o o
a = 2.5
b =0.125
A , = 0.50(225-775)
c2=0.l C3=0.l
o Wygnanski S Fiedler
-1.5
-I
.5 0
z-z
Z.25~Z.75
.5
1.5
Figure 5.11. Comparison of experimental results and model predictions
for the sidewise fluctuations in a free shear layer
-------
5-23
u'w'
w2
wmax
.028 r
.024-
.020-
.016-
.012-
.008-
.004 -,
a = 2.5
b =0.125
A|=0.50r.5
o Wygnanski & Fiedler
Figure 5-12. Comparison of experimental results and model
predictions for the shear correlation in a
free jet
-------
5-24
.025
.020
u'w'
Umax
.015
.010
.005
oWygnanski d Fiedler
as measured
• Corrected Wygnanski
8 Fiedler
a = 2.5 c2 = 0.l
b =0.125 C3=0.l
A,=0.50(z25-z.75)
Tollmien
-1.5 -I -.5 _0
z -z
z.25"z.75
.5
1.5
Figure 5.13.
Comparison of experimental results and model computations
for the shear correlation in a free shear layer
-------
5-25
is fair. It should be noted that the experimental values of
shear correlation from Ref. 21 have been shown as reported (the
open circles) in Figure 5.13 and also as corrected by us (the
solid symbols) so as to agree with the measured rate of spread
of the free shear layer. A comparison of the measured shear and
that inferred from the mean velocity profiles was reported by
Wygnanski and Fiedler but apparently their computations contained
an error. Also shown in Figure 5-13 is the level of shear that
may be inferred from the mean spread of the free shear layer
studied by Tollmien [Ref. 25] and Prandtl [Ref. 26] many years
ago. It is seen from the results presented in Figures 5.12 and
5.13 that the model gives a fairly good representation of the
shear in both the free jet and the free shear layer.
A careful study of Figures 5-6 through 5.13 shows that it
really.is necessary to study further the problem of choice of
model parameters. However, before this is done, it appears
desirable to have at hand experimental data,which one can rely
on to be truly representative of the basic flow which is being
calculated. It is difficult to choose a more sophisticated
model until the question of the degree of isotropy on the center-
line of a free jet'is settled. In addition, one should at this
point determine if the model just found for free shear layers can
be used for a model of the outer regions of a boundary layer and
give reasonable results.
Before turning to the problem of the turbulent boundary
layer, it will be instructive to find a relationship between the
values of A, used in the free shear layer and the free jet
calculations and the general magnitude of the integral scales
measured for such flows. In the computations that have been
made, it has been assumed that A, is constant across a free jet
or a free shear layer at any given longitudinal position and,
in magnitude, proportional to the local scale of the mean flow.
It is well known that the'integral scales of such flows are, in
general, proportional to the local mean scales, but the actual'
value of the integral scale varies across the layer.
-------
5-26
In Table 5.1 we present the values of integral scale within
a free jet, as reported by Wygnanski and Fiedler. The integral
scale tabulated is the longitudinal integral scale
L = w'(z,)w'(z0) d(zD - z, ) (5.28)
. . / \ f X £- C. J-
w'w' ( z-, ) JQ
for the free jet.
Table 5.1. Integral Scales in a Free Jet [Ref. 21]
Radial Position Dimensionless Scale Scale Ratio
r/x L/r
0 • 0.448 1.12
.05 0.595 0.84
.10 0.726 0.69
.15 0.850 0.59
.20 0.855 0.58
Also shown in Table 5.1 is the ratio of the computational scale
A, to the local integral scale L . Thus a typical value for
this ratio for the free jet is
A1/L = 0.69 (5.29)
For the free shear layer, similar results are given in
Table 5.2. These experimental values are also due to Wygnanski
and Fiedler. The longitudinal integral scale is, in this case,
defined by
L = u'(x )u'(x?) d(x? - x ) (5.30)
u'u'(x1)Jo J-
Table 5.2. Integral Scales in a Free Shear Layer [Ref. 22]
Location in Jet Dimensionless Scales Scale Ratio
L/x L/(z 25~z 75) A-^/L
Inner Region
Center
Outer Region
0.098
0.103
0.147
0.846
0.883
1.27
0.59
0.57
0.39
-------
5-27
A typical value of A,/L for a free shear layer appears to be
approximately
A1/L = 0.55 (5.3D
Before proceeding further, it must be demonstrated that if
the present model is applied to a boundary layer, useful results
will be obtained for the same choice of model parameters that has
been made for free turbulent shear flows .
Application of Model to Boundary Layers
If the model of turbulent shear flows is to be applied to
a boundary layer, the parameters c?, c_, a , and b are known.
But, since the characteristic length in a boundary layer is
arbitrary (as it is in the free jet and the free shear layer),
we are at liberty to choose c, , i.e., the ratio between A
and the characteristic length (which, in this case, we take
equal to 6 , the thickness of the layer in which the
. y y
velocity reaches 99% of its free stream value).
As discussed previously, one other parameter enters the
problem, namely, a , the coefficient appearing in (5.19). We
have, then,
AI = av/~a z (5.32)
for 0 < z < c, 6 OQ/(ov/a) and
-L « y y
A! = c16>g9 (5.33)
for z > c-^6 99/(o/a)
With only these two parameters ex and c, to determine,
the search is not difficult. The model that has been found is
the following:
a = 2.5
b = 0.125
C2 = °'1 (5.34)
c,, = 0.1
= Al/6.99
= 0.443
-------
5-28
The ability of this model of a turbulent shear layer to
predict the known mean properties of turbulent boundary layers is
shown in Figures 5.1^ through 5.16. In Figure 5.1^, we show the
skin friction developed by our model as it proceeds from a
disturbed laminar layer to a fully turbulent layer. Also shown
•
are the laminar skin friction law and the turbulent law proposed
by Coles [Ref. 27] which is a good.fit to experimental data. It
is no great surprise that the general levels of skin friction we
computed agree well with experimental findings inasmuch as the
values of a and c, were chosen to get these levels correct.
Of more importance is the nearly exact following of the trend of
skin friction with Reynolds number by the model computations.
Figure 5.15 shows a comparison of the computed mean velocity
profiles developed by the model in the vicinity of the wall and
the well-known law of the wall as proposed by Coles [27]. It may
be seen that the law of the wall is not quite achieved by the
present selection of model parameters. However, the results are
sufficiently accurate to be encouraging.
In Figure 5.16, we compare the experimentally determined
velocity defect law proposed by Coles [271 with the results of
our model computations. It is seen that once the turbulent
boundary layer is well established, the computational model
gives a fairly good representation of the outer regions of the
turbulent boundary layer.
With these results in hand, we must now consider the
relationship of the computational scales used to the longitudi-
nal integral scales that are found in the outer regions of
turbulent boundary layers. For this purpose, we may use the
measurements of Grant [Ref. 28]. When the experimental correla-
ft
tions reported by Grant for y/6 =0.66 in a turbulent bound-
ary layer are integrated to give the longitudinal integral scale,
one obtains L/6 <*• 0.3- Since for our calculations, 6 /6 QQ =*
0.83, we find that L/6 QQ =* 0.25. Since the computational scale
*Grant defined 6Q as that height in the boundary layer where the
velocity defect was equal to the friction velocity.
-------
5-29
,-2
10
5XIO"3
2XIO'3
10-31
2XI04
COMPUTED VALUES FOR
a=2.5 b=O.I25 A,=O.I5S99
:9 = 0.l c, = 0.l a = 0.443
TABLE H OF COLES
I
I
I
sxio4 io9 2xio5
I
5XI05 IO6 2XI06
5XI06
Figure
Rex
Computed variation of coefficient of friction with
Reynolds number
-------
5-30
U/Ur
30
25
20
15
10
COMPUTED VALUES FOR
0 = 2.5 b = O.I25 A,=O.I5899
C9 = 0.l c, = 0.l a = 0.443
Cm O
o Rex=0.277X|Q6
<
o Rex = l.62Xl06 0
A Rex=6.76X|Q6
TABLE
OF COLES
10
io
10-
Figure 5.15-
Computed velocity profiles for three Reynolds numbers
Computation started at Rex = 20,000. The reference
velocity u is defined by T ,, = pu2 .
T waxj. r T
-------
5-31
u-Uedge
COMPUTED VALUES FOR
= 2.5 b = O.I25 A2 = 0.15899
C2=0.l c3=0.l a=0.443
Rex=0.277X|06
= I.62X|06
6
Re* =6.76X10
TABLE
OF COLES
-10 -
-12
0.02
0.05
Z/Zedge
Figure 5.16. Computed velocity defects for three Reynolds numbers
Computation started at Rex = 20,000. The reference
=
velocity UT is defined by
= Pu
-------
5-32
used was A,/6 = 0.15, we find that
A-J/L -0.6 (5.35)
This is a most welcome result since it shows that, for all the
turbulent flows we have investigated, the ratio of the proper
computational scale to the longitudinal integral scale is approx-
imately the same.
Some Comments on Second-Order Modeling Techniques
The method of modeling turbulent shear flows which we have
just described was developed, as we have previously pointed out,
in order to attempt calculations of turbulent flows other than
the classical shear layers that were discussed in the previous
sections. The author and his colleagues have applied the model
to the calculation of the decay of a turbulent line vortex [Ref.
29], to the generation of turbulence in the earth's atmosphere
[Refs. 16 and 20], and to the dispersal of pollutants by the
atmosphere [Ref. 30]. Since these computations were carried out
with the original oversimplified constant A model discussed
previously, one must not take the numerical values obtained too
seriously; nevertheless, these computations did give some most
interesting results and insights. Certainly the utility of the
method was demonstrated and it appears that it, and others like
it, should be carefully studied and refined in the next few years
The first order of business should be a continuance of the types
of parameter studies that we have just described, for as large a
spectrum of shear flows as can be reliably measured. In this
way, we would hope to develop a model with the broadest capa-
bility possible.
It is here that one runs into the difficulties that were
touched upon previously in connection with the free jet and the
free shear layer. It would be very helpful indeed if the
research community could agree on canonical free jet and free
shear layer experiments which could be performed by several
investigators. The purpose of these experiments would be to
-------
5-33
strive for agreement between experimentalists as to what the
characteristics of such flows were and to explain any discrep-
ancies that might exist. Before very much refinement of
turbulent shear flow models can be accomplished, it appears that
we are going to have to have a more precise definition of what
the models must predict.
-------
6-1
6. AN EQUATION FOR THE SCALE h
Although our work has not proceeded so far as to couple
an equation for the scale A to the rest of the model equations
when making shear layer computations, it is our intention to do
so in the near future. For the sake of completeness, then, we
outline here the derivation of an equation for the scale length
Ai •
We have seen in Section 5 that A-, is related to the
integral scale L defined by (5-30). Since this scale depends
on an integral of the two-point correlation
we will start with an equation for this quantity.
Starting with (3.15), equations for (u.J), and (U/.)R
1 1\ f\. -D
may be written. Much as in the derivation of (3-17)» these may
be combined to give
'du,\ /du,c]
~'B
. .
77L(ui)A(uj)BB
\ J ' B
-------
6-2
In this notation, a variable at one point is considered a
constant with respect to differentiation at another point. For
example,
/-\ \ / -x \
C(u').(u')D] (6.3)
We introduce new variables. Let
(xk}B]
(6.4)
and
'k
(6.5)
In terms of these variables, the contraction of (6.2) can be
written
d ,-
5t (U:
:U.)A(U.)B
(UJ)B
£(u\)
A
i'B
A+4
2 dy
(u{)A(uj)B(ui)B
- (ui)A(uj)A(ui)B
'B
pA(u,
A
vo S2
I e T "I
e>j -
CiiM CnM + —=- T'CiiM + T'fuM
••^^j i;A^Ui;B TQ[_iAlUi;B iB^Ui;AJ
(6.6)
-------
6-3
This equation may now be 'integrated with respect to £,
from 0 to +°° . The process is discussed in Ref. 31 and the
reader is referred to this paper for more detail. The result-
ing equation is very complex and it is not necessary to repro-
duce it here, since it is clear from the structure of (6.6)
that the simplest model equation 'for A, , for the case of a
simple parallel shear flow, will be of the form
77TT7TA
U W A
Dt 1 l 9z 2 9z
2
q A
_
+ —? (qA ) - C 2v — jy + CM - w'T'A (6.7)
^ 3z^ -1 ^ X^ q o L
Here the C's are constants that remain to be determined by
matching computed solutions with experimental data.
This equation is almost identical to that used by Rod! and
Spalding . [Ref . 15] in their calculation of turbulence in free
jets. There is, of course, an additional term, namely,
^ w'T'A
o
which represents the tendency of atmospheric stability to
increase or decrease the scale of atmospheric turbulence
depending on whether the atmosphere is unstable, with w'T' > 0,
or stable, with w'T1 < 0.
It is clear that one might derive equations .for a whole
host of scales. For example, equations might be derived for
the integral scale of the temperature fluctuations by writing
the equation for T'T' and proceeding as above, or the integral
scale of the concentration fluctuation by use of the equation
How might these equations be used? Since, as we have seen,
we always need to have available in our calculations a length
A, (or several A,'s, if a more complicated model of turbulence
-------
6-4
is desired), we could, rather than relate these A's directly
to some locally defined mean scale of the motion, carry along
equations, such as (6.7) for the required A 's and add these
equations to the group of equations that are to be solved
simultaneously.
Although such a procedure should be the ultimate aim of
any second-order closure scheme, we will not try, in this report,
to couple an equation for the scale A-. to the equations for the
second-order correlations. We have chosen to proceed in this way,
as mentioned previously, because it permits one to observe the
sensitivity of the modeled equations to various choices of A-, .
An ability to observe this sensitivity has been and continues to
be useful in studying the general characteristics of solutions to
turbulent shear flows by invariant modeling.
-------
7-1
7. RELATION OP SECOND-ORDER MODELS TO K THEORY
Before we begin a discussion of the calculation of
atmospheric flows using the invariant second-order modeling
scheme that has been described here, it is instructive to relate
this model to the classical eddy diffusivity or K models of
turbulent transport. To do this, we first write out the model
equations [(4.17) through (4.26)] for an essentially parallel
shear flow in the atmosphere; that is, we assume u = u(z) ,
v = d/dy = 0 , w « u , and d/dx « d/dz . The result is
(7.:
DT aiT . dw'T
r\4~ ~" V r\ "~ N.
DC c cT
_ ^ - v _ — . _ —
Dt o 2 5z
(7.5)
Dt
_
0 dz2 ° X2
(7.6)
''
Dv'v
Dt
(7.7)
-------
'w1
w
Dt
T
W'T
'o
0 9
2
9w'w
'w'
Du'T'
Dt
= - u'w
''
- w'T' . + - A q
oz oz \ ^ dz
du'T
'1
u'T
+ 2 |_(A q %L
dz \ ^ dz
1 ^ LA
'N— P A0
p. dz \ Ko 3
w'T
7-2
- 2v
''
W'W
(7.8)
_ 2v
X
(7.9)
+ v
0 Sz2
- 2v
u'T1
0 X2
(7.10)
'o ^_2
- 2v
w'T
0 x2
(7.11)
-------
Dt
= -2*-^^-+ 4- A0q
Du'ca
Dt
= - u'w'
dC
^ - w'C' ^ + ^- I A9q x
dz a dz dz \ ^ dz
a
A U'C'
An a
Dw'C!
dC
= - w'w'
a
1 d
dw'C'
OL>
Sz~
o N_2
- 2v
w'Ca
In addition, we have u'v' = v'w' = v'T' = v'C'.= 0
\JL
7-3
+ v
d2T'2
0 dz2
- 2v
t *-
(7.12)
+ v
g
0 dz'
- 2v
u'C'
a
0 x2
DC'T'
cz
Dt
d2C'T'
C'T'
-V
X2
dz
(7.15)
DC'
Dt~
dC
= _ 2w'C'
a
a dz
S
^
dC
,2
a
,2
- 2V
(7.16)
-------
7-1
Equations (7.1) through (7-5) define the mean variables
u, w, T, C , and p . In these equations the second-order
transport correlations u'w1 , w'T1 , and w'C' appear. In
\JL
older methods of turbulent transport calculations, these correla-
tions are related, by analogy to molecular transport, to the
gradients of the mean velocity, temperature, and concentration
fields. This assumption is equivalent to saying that the second-
order transport correlations have no "memory," i.e., do not have
dynamic or rate equations which govern their behavior. This
assumption also implies that the second-order correlations at a
point have no knowledge of what is going on in contiguous elements
in the flow. Since we may look at each of (7.6) through (7.16) as
a rate equation for the particular species in question (in the
sense of chemical kinetic equations), we should be able to recover
traditional turbulent transport theory (eddy transport theory or
K theory) by a suitable limiting process applied to (7.6) through
(7.16).
The suitable limiting process requires the following
assumptions:
(1) Since the dynamics of the correlations (species) are
not important, we will make the assumption that the correlations
(species) are in local equilibrium. That is, we will assume, in
a manner precisely analogous to the assumption of local chemical
equilibrium in a flowing gas, that
§t30 (7'17)
(2) In view of the fact that in eddy diffusivity theories
the turbulent transport at a point has no knowledge of the value
of transport at another point, we must neglect all the diffusive
terms in the dynamic equations for the correlations.
(3) Since eddy diffusivity or K theory is a large
Reynolds number theory, we must assume in arriving at our equa-
tions that the Reynolds number is large. We will assume, there-
fore, that
-------
7-5
Al _ VoAl
a +
o
Thus, to relate our present model to K theories, we
must seek the equilibrium, nondiffusive limit of (7.6) through
(7»l6). Because of the high Reynolds number, nondiffusive, and
equilibrium assumptions, we refer to the resulting equations as
the "superequilibrium" equations. This word aptly describes the
flow we are studying.
The superequilibrium equations for the correlations are
0 = - 2uTw~T" |n. - 9-- (i + 2b)uTuT + a - (7.19)
oz A-j^ 3A
Q _ Q3
0 = - 9- (1 + 2b)v'v' + * - (7.20)
Al 3A1
0 = %& wTTT - 9L (i + 2b)wTwr + 3 - (1.21)
To Al 3A
0 = - v^w1" + - \TnFr - - (1 + 2b)uTvTr (7.22)
0 = - UTWT - w7^ - - (1 + 2b)uTTT (7.23)
0 = - w^w"5" -+~T'--(l + 2b)wTT~r
dz T A
0 = - 217^0 - 2b - T'2 (7.25)
dz A]_
0 = - iT'lv7" ^—^ - vTrCT |^- - SL_ (i + 2b)CTuT (7.26)
5z a 5z A.J_ a
_ o'C" _ _
0 = - w'w1 r— ^ + |- C'T1 - 9- (i + 2b)w'C' (7-?7)
dz TQ a A1 a
-------
7-6
m
0 = - w'C' 2L _ wiTi a _ 2b q_ C,T, (7.28)
a dz dz A-^ a
O O rt~
0 = - 2wTCT Tr-2. _ 2b 3- C' (7.29)
a dz A1 a
(7.19) through (7-29) are algebraic equations for all of
the correlations. It is easily seen that (7.19) through (7.25)
are a coupled set of equations for the local state of turbulence
and turbulence transport in the absence of pollutants. This
occurs because we have assumed, in the derivation of the basic
equations from which the present set was obtained, that the
pollutants are passive. Knowing the state of turbulence from a
solution of (7.19) through (7.25), we can obtain the vertical
contaminant transport correlation C'w1 from a simultaneous
solution of (7.27) and (7.28). The correlation u'C^ is obtained
directly from (7-26) once w'C T is known and, finally, the
variance of the contaminant fluctuation is found directly from
(7.29).
The solution of (7.19) through (7-29) can be obtained in
a particularly elegant fashion if the following definitions are
introduced into the equations. Let (for the case when du/dz > 0)
= UUA2 l~f}2 u^"= UTA2 |^T
1 \ dz / 1 dz
(7.30)
= WWA? |H
1 V dz
u'w1 = UWA" ^
1 \ dz
-------
7-7
arid likewise for the concentration correlations
• o an *^ni P am '
= UCA? 4^ vr^- C'T' = CTA, ^
-5-- - .—
a 1 dZ dZ a 1 dZ dZ
(7.3D
""
uT^T = WCA2 |u o c,2 = CCA2 a
a 1 9z 9z a 1 I 9z
Note that
QQ = Q2 = UU + VV + WW (7-32)
Substitution of the definitions given in (7-30) and (7.31) into
(7.19) through (7-29) results in
o3
Q(l + 2b)UU = | 2UW (7-33)
O3
Q(l + 2b)VV = |- (7.3^)
O3
Q(l + 2b)WW = |- + 2R1WT (7-35)
Q(l + 2b)UW = - WW + RiUT (7.36)
Q(l + 2b)UT = - UW - WT (7.37)
Q(l + 2b)WT = - WW + RiTT (7.38)
Q(2b)TT = - 2WT (7-39)
for the atmospheric equations while for the contaminant equations,
we have
Q(l + 2b)UC = - UW - WC (7.^0)
Q(l + 2b)WC = - WW + RiCT (7.^1)
Q(2b)CT = - WT - WC (7.^2)
Q(2b)CC = - 2WC
-------
7-8
In these equations, Ri is the Richardson number
Ri = ?- £•/ (sr- 1 (7.H)
o
It is immediately obvious from (7-33) through (7.43) that all the
nondimensional second-order correlations are a function only of
the Richardson number and the parameter b from the second-order
closure model. It will be remembered that the value of b
determined in the parameter search reported previously is 0,125-
It is convenient to express the solution of (7-33) through
(7.^3) in terms of the parameter
P = 1 - (4 + 15b)Ri + [1 + 2(2 - 9b)Ri + (4 + 9b)2Ri2]1/2
(7.45)
In terms of this parameter, the various correlations may be
written
(7.^6)
o2 -
^
T1TT = (*
1 r
? r
;i + 2br
3 + bRi)[P + (1 + 4b)Ri] + 2b[P + (
1 + b)Rl]
3(1 + 2b)(P + bRi)[P
2b)
uw = _
u
3(1 + 2b)[P + (1 + T
P + (1 + b)Ri
3 (P + bRi)[P + (1 + 4b)Rl] *
TIT = b 2P + (1 + 2b)Rl
u 3(1 -i- 2b) (P + bRi)[P + (1 + 4b)Ri]
-------
7-9
WT " ~ Q (7-52)
b
3LP + (!•+
1
3CP + (1 + 1
b
Q3
_2
b)Ri] y
2P + (1 + 2b)Ri
TT " Q (7.53)
3(1 + 2b) (P + bRi)[P + (1 + 4b)Ri]
WC = - 3[P + (1 + Q (7'55)
CT = 3[P + (1 + ^b)RlJ Q (7'56)
CC = 3[P + (1 + i|b)Ri] Q (7.57)
It is clear from these equations that when the parameter
p
P = 0, there is no turbulence (Q = 0) and all the second-order
correlations vanish. The critical value of the Richardson number
for which this occurs is a function of b and is given by
R1crit = 4b(l (7>58)
For b = 0.125, we find that the critical Richardson number is
Ricrlt = 1.636 (7.59)
The behavior of all the nondimensional second-order correl-
ations as functions of the Richardson number are plotted in
Figures 7.1 through 7.5. From these figures the profound differ-
ence between turbulence and turbulent transport in stable and
unstable atmospheres is obvious. Note particularly that the non-
dimensional vertical transport of matter and heat fall off far
more rapidly than do the nondimensional turbulent energy compon-
ents when a stable atmospheric situation is approached. In fact,
-------
7-10
lOr
oB
WW
uu=
U'U'
vv=
ww=
Critical Ri
1.636
Richardson No., Ri
Figure 7.1.
Superequlllbrium values of the turbulence
components as a function of the Richardson
number for b = 0.125
-------
7-11
Q(
20,
16
12
8
-2
~2_ U7U7 +V/V/+W/W/
Critical Ri
1.636
-1
Figure 7.2.
Richardson No., Ri
Superequllibrium value of the sum of the
squares of the turbulent velocity fluctuations
as a function of the Richardson number for
b = 0.125
-------
7-12
2.Or
o
z>
05
ID
I
UW=
-uw
UT=
u'w'
u'
A 2 du dJ
A
UT a uc
Critical Ri
1.636
Richardson No., Ri
Figure 7-3
Superequilibrium values for shear and longitu-
dinal heat and mass transfer correlations as
functions of the Richardson number for b = 0.125
Note UT and UC approach the value 1.328 for
Ri = - oo.
-------
7-13
o
£
i
(0
WT =
A <3T
A
wc=
w'c'
dz dz
-WT a -we
Critical Ri
1.636
Figure
Richardson No., Ri
Superequilibrium values for vertical heat
transfer and vertical mass transfer correla-
tions as a function of the Richardson number
for b = 0.125
-------
7-14
O
O
05
O
•*
h-
I-
9.387
TT =
T'1
A2^)2
AlV37/
CT =
dz
cTT
A
A
CC =
c'c'
A
•\dz/
Critical Ri
1.636
Richardson No., Ri
Figure 7.5
Superequilibrium values for TT, CT, and CC
as functions of the Richardson number for
b = 0.125. Note that these functions approach
the limit 9.387 for Ri = - «>.
-------
7-15
above a Richardson number of one, vertical turbulent transport
has almost ceased to exist although there is still some atmos-
pheric turbulence.
It should be noted that the superequilibrium results just
obtained specify the nondimensional values of second-order
correlations. For example, assuming the value of b = 0.125 to
be correct, a Richardson number of 0.10 would give
UW = - 0.2712 (7.60)
WT = WC = - 0.2631 (7.61)
The transports of momentum, heat, and matter would then be given
by (for 3u/9z > 0)
(7.62)
'o°P0 ^ ' 0-2631P0A (7-63)
-po w^-= 0.2631P0A ^L (7.64)
It is clear from these expressions that the actual transport is
not defined until the length scale A, is known. This is a
difficulty with atmospheric flows, for unless A, is determined
at a given altitude and the local Richardson number specified
there, the transports are not known. In general, A will depend
at a given altitude on the Richardson number but can assume a
range of values depending on the past history of the motion.
While this range of values is limited so that the order of magni-
tude of the transport might be determined, there will always be a
variation in transport proportional to the square of the variation
in A, at any fixed Richardson number.
For classical laboratory flows, this problem does not
exist. In this case, it is generally found that A, is propor-
tional to the characteristic breadth of the layer under consider-
ation while 'the gradients are proportional to a characteristic
-------
7-16
velocity, temperature, or concentration difference divided by
this characteristic breadth. Thus, for the classical shear
flows
o / ^ - \ 2 o / Au , \ 2 , . „
A2/^M = 62 -r^1 = const AU 2 (7.65)
char \ dchar/ l char;
and likewise
A2 |u |T = congt A~ AT (7.66)
1 dz oz char char
and
A-, ^~ x—21 = const Au , AC" , (7.67)
1 dz dz char char
For each type of flow, these constants are well defined. This
type of simplicity is, alas, not true of the atmosphere,as we
will try to demonstrate presently.
Before considering the problem of the scale A-, , it is
instructive to compare the results of superequilibrium theory
with certain well-known results from classical turbulent trans-
port theory Tor^the case when no gravitational effects are
involved. To do this, we place the Richardson number equal to
zero in the expressions given in (7-^5) through (7-57). We
obtain, for b = 0.125,
P = 0.3333 (7.68)
Q2 = - - - ? = 1.7066 (7.69)
3b(l + 2b)
UU = — a + 6b - - = 0.7964 (7.70)
9b(l + 2b)3
VV = WW = -- - - 5- = 0.4551 (7.7D
9b(l + 2bK
UW = TW = CW = -- . - _ = -0.2786 (7.72)
9(1 + 2b)J
-------
7-17
UT = UC = — T = 0.3413 (7.73)
3(1 + 21))-3
TT = CT = CC = =- = 1.7066 (7.74)
3b(l + 2b)
First we note that superequilibrium theory indicates that
v'v' = w'w' and, further, that
6b =
v'v' w'w' VV WW
Second, we note that the value of - u'w'/q , which
Bradshaw, Ferriss, and Atwell [Ref-. 8] assume to be a constant
equal to 0.15, is defined by superequilibrium theory to be
- ™ = 1 I gb/l = °'163 (7.76)
Q
This is a rather surprisingly accurate result in view of the
fact that the value of b was determined from very different
considerations during the parameter search reported in Section 5
We may also derive from superequilibrium theory the value
of von Karman' s constant /c in his expression for the turbulent
shear near a surface, namely,
~ 9 / >, - \2
T = -pu'w' =p/c
-------
7-18
a = 0.7. Letting
. A]_ = 0.7z (7.79)
we find that (7.78) gives
T = 0.^9/lTb
t 9(1 + 2b)3
Comparison of (1.11) and (7-80) reveals that
9(1 + 2b).:S
or
K = 0.37 (7.82)
The value of von Karman's constant is actually 0.4. Again, the
agreement between results obtained by taking the equilibrium,
nondiffusive limit of our second-order closure model of turbu-
lent shear flow and classical mixing length theory is rather
remarkable.
We now return to a discussion of the scale A^ . We ask
whether superequilibrium theory might give some useful informa-
tion about this scale. The answer is obtained by considering
the equilibrium, high Reynolds number, nondiffusive form of
(6.7) . The result is
n _ ,-, —;——r . oU „, .-, 3 A i r> S , i m i A I n Q o \
0 = - C.u'w'A-j r~- 2bC_q A-, + C^ S— w'T'A-, (7.oj)
1 J- 3Z J -L Q
Substitution in (7.83) of the definitions given in (7.30) results
in
We note that the scale A-, may be cancelled from this equation
We may imply from this that the scale of turbulence cannot
-------
7-19
depend upon purely local conditions but must depend to some
extent upon the nearby elements of the flow. In the classical
laboratory free shear flows (which are self-similar), the scale
at two equivalent points is tied to the scale of the mean motion
at these points, because the neighboring regions about these two
points are equivalent; thus the effects of past history and
diffusion are equivalent. This is not so in the atmosphere and,
hence, the scale A is not well defined.
In dealing with the atmosphere with classical K theories,
we must then not only know the effect of local Richardson number-
on the nondimensional parameters UW , WT , WC, etc., but must
also have some general information about nearby elements of the
flow so that a scale can be determined.
The fact that the scale is undetermined under superequili-
brium conditions is consistent with the experimental results of
Rose [Refs. 33 and 3^] who examined the behavior of turbulence
introduced into regions of uniform shear. In these experiments
it was found that the scale length tended to increase contin-
uously as the observation station was moved downstream [Ref. 33]
and that the level of turbulence, when the scale was large
enough so that the turbulence was not dissipated, followed the
scale of turbulence introduced into the shear layer.
A characteristic length used in similarity theory of atmos-
pheric flows is the Monin-Obukhov length defined [see Ref. 32]
as
T
L = - -3L
I T.T I
u'w
3/2
(7-85)
K gW ' T '
Substitutions of the definition given as (7-30) into (7.85)
gives
IIW|3/2 An
MJ -- L_ (7.86)
WT KRi
If we substitute the expressions for UW and WT from (7.50)
and (7.52) into (7.86), we obtain an equation for the ratio of
-------
7-20
the Monin-Obukhov length to the scale A-, , namely,
L _ 1 /b \1/2 [P + (1 + b)Ri]3/2
T~ = 7 ( T ) - T75 - T7? (7 .87)
Al K ^/ Ri(p + bRO^CP + (1 + 4b)Ri]1/£f
This relationship is plotted in Figure 7.6.
The singular character of the relationship between L and
A-, at Ri = 0 and the very rapid variation of L/A, in the
range of a Richardson number of general interest is apparent.
The results presented are interesting but do not, alas, permit
any further insight into that elusive parameter - the Monin-
Obkuhov length.
For those familiar with atmospheric turbulence, the critical
Richardson number of 1.636 found in this superequilibrium study
is obviously • too high. We have investigated the cause of this
result in considerable detail. It seems clear from the excell-
ent results of superequilibrium theory for neutral flows that
this method of looking at the model equations is a powerful one.
Indeed, it has been found that it is an excellent basis from
which to decide upon the basic modeling that is to be used. We
are presently changing our model so that the superequilibrium
limit of the modified model will yield critical Richardson
numbers more in agreement with experimental observations .
-------
7-21
L/A,
I4r
12
10
Critical Ri
1.636
-I
Figure 7.6.
Richardson No., Ri
Relationship between Monin-Obukhov length L
and the isotropy scale A^ as a function of
the Richardson number as given by superequili-
brium theory for b = 0.125
-------
8-1
8. APPLICATION OP SECOND-ORDER MODELING TO THE ATMOSPHERIC
BOUNDARY LAYER
We will now take up the application of the second-order
closure model of turbulent flow developed in Sections 4 and 5
to the problem of computing the atmospheric boundary layer.
Before doing this, however, it is important to point out once
again that the model used contains a minimal number of scale
lengths. It is felt that the results obtained with this simple
model are most instructive. One of the goals of our present
studies is to obtain an answer to the question of just how
sophisticated in the specification of the scale lengths one
must get to obtain reliable predictions of the properties of
turbulent flows in which the profiles of the mean velocity,
temperature, and concentration are not related in a simple
fashion as they are in the classical free turbulent flows.
To demonstrate the method, we will compute the state of
atmospheric turbulence for several cases of motion in the
atmospheric boundary layer for which measurements have been
made of the mean Velocity and the mean temperature profiles
and of the average values of the second-order correlations
u'u', v'v1, w'w', u'w', u'T', w'T', and 1" . These data
obtained for the Air Force Cambridge Research Laboratories
were provided by Messrs. Wyngaard and Cote of APCRL. The pro-
files of mean velocity and mean temperature for the three cases
considered here are shown in Figures 8.1, 8.2, and 8.3.
To calculate by means of our model the turbulent energy
and the fluxes of momentum and heat associated with these
profiles, we will assume the mean velocity to be parallel to
the surface in the xz plane. Thus, v = w = 0. We will also
consider the- motion to be independent of x and y and hence
a function only of z (height) and t (time). Consequently,
d/dy = c3/dx = 0, and it is seen that the substantial (total)
derivative D/Dt = d/JH . The mean velocity and temperature
equations are then analogous to Eqs. (7.1) and (7.2) and are
-------
8-2
120
100
80
60
e
«*
N
40
20
Figure 8.1.
\
\
\
o Measured velocity
o Measured temperature
468
u, meters/sec
-i
-2
-3 -4
T, °K
-5
-6
Mean velocity and mean temperature profiles for a stable
atmospheric boundary layer (Wyngaard & Cote, Run 25).
The solid lines are faired through the data points; the
dashed curves are assumed continuations of the
measured profiles for computational purposes.
-------
8-3
120
100
80
tn
-------
8-4
120
100
80
60
-------
8-5
_ _ _
u o* u du'w' 1
_ _ y - _ _ _ - _ -
t o ^Z2 3Z PO
_ „
= v — 5- - ^ - + Q
In (8.1) the pressure gradient dp/Sx is included. Addi-
tionally, in (8.2) a term Q(z,t) has been added which repre-
sents any local source of heat that might be present in the
atmosphere. The purpose of this term will be explained below.
The equations for the turbulent quantities are identical
with (7.6) through (7-12) except of course that D/Dt = d/dt
on the left-hand side of all the equations.
In order to calculate the turbulence produced for the
three cases of atmospheric motion given in Figures 8.1, 8.2,
and 8.3, we may proceed as follows. First, we note that the
profiles given were steady, at least on a time scale large
compared to the time required for turbulent adjustment accord-
ing to our model, as will be demonstrated later. If this is
so, we may, in the equations for the double correlations, assume
that u(z) and T(z) are the given profiles. These equations
can then be solved for the time variations of these correlations
assuming any initial distributions for the turbulent correla-
tions themselves. The solution is continued until the deriva-
tives of the correlations with respect to time approach zero.
The resulting distributions of the second-order correlations are
then taken to be the appropriate distributions for an atmosphere
whose mean profiles of u(z) and T(z) are those given.
Equations (8.1) and (8.2) may then be solved, assuming c^u/dt = 0
and df/dt = 0 for the distributions in z of the driving
force (1/p ) Op/c)x) and the heat source Q necessary to maintain
the atmosphere in its assumed state.
It is found that the steady state or equilibrium distribu-
tions of the second-order correlations found in this manner are
independent of the assumed initial conditions on these variables.
-------
8-6
For this reason, computations are generally run with a small
level of isotropic turbulence introduced into the system, while
2
the•other correlations u'w1, u'T", w'T', and T" are assumed
to be zero initially.
These calculations can be illustrated by considering the
case of the almost neutrally stable atmosphere shown in Figure
8.2. For this case, we choose the scale length A-, in the
following way. Near the surface (the flow is unstable only
very near the surface), we choose A, = 0.7z in accordance
with the findings of our parameter search for the incompressible
boundary layer. Also in accordance with the parameter search
for the flat plate boundary layer, we take A constant in the
outer layers of the atmospheric boundary layer and, as a first
guess, equal to 0.15 times a characteristic height of this layer
This is a somewhat difficult and arbitrary dimension to choose
since the actual measurements of the mean velocity profile do
not permit such a height to be chosen, and one must do the best
one can with a continuation of the measured data. For this
case it would appear that a characteristic height might be
chosen as something like 110 to 120 meters. For this reason,
we take, as a first guess at the scale length in the outer
region, A, = 17 meters. With these assumptions, a computation
of the turbulence and transport properties of the'atmospheric
boundary layer may be made according to the scheme outlined
above. The results obtained are shown by the solid lines in
Figures 8.4 and 8.5. It is seen that the order of magnitude of
all quantities is estimated correctly. The agreement between
the calculated and measured values of u'w1 and w'w' is
reasonably good. The agreement for u'u' and v'v' is not
so good. It is interesting to note in this regard that the
magnitudes of u'u' and v'v' as measured are alike, with v'v
more than twice w'w1 . In view of the structure of the govern-
ing equation, this latter experimental result is certainly
surprising. The agreement between the measured results and
computed values for the temperature correlations u'T', w'T'
-------
8-7
30
-------
8-8
30r
g 20
o>
10
.02
.04
.06
.08
.10
T"
Figure 8.5.
Steady state distributions of the temperature correlations
for profiles of mean velocity and temperature given in
Figure 8.2. Near the surface, A]_/z = a = 0.7. The solid
lines are for Aj max = 17 meters. The dashed lines are
for A, m^v = 30 meters.
-L ITlciX
-------
8-9
2
and T1 is poor. The measured correlations are, however,
very small and the accuracy of these results may be somewhat
questionable. For example, the measured vertical heat flux
correlation is constant with height and extends far past the
region where any significant potential temperature gradient
exists.
An idea of the sensitivity of the computed results to
the choice of the outer or maximum value of A, can be had by
referring to the dashed lines in Figures 8.4 and 8.5. These
results were obtained for an outer A, of 30 meters. The
a = A /z near the wall was maintained at 0.7. It is seen that
there is not a great change in the computed results below a
height of 20 meters, as might be expected, since the A 's in
both cases are equal below 17/0.7 = 24.3 meters. For the
temperature correlations, there is essentially no chance since
all potential temperature gradients are confined to the region
below 24.3 meters.
An idea of the sensitivity of the computed results to the
choice of a = A,/z near the surface can be had by reference
to Figures 8.6 and 8.7. Here the results shown by dashed curves
were obtained with A../Z = a =1.0 near the surface and
A, = 17 meters. For reference the results obtained with
1 max
a = 0.7 and A, = 17 meters are again shown as solid
_1_ ITlcLX
curves. The effect of increasing a is very pronounced, as
might be expected. It is interesting to note that even with
this choice of parameters it was impossible to match the magni-
tude of v'v? and the 'temperature correlations. The magnitudes
of u'w' and . w'w' are now greatly overpredicted. Before making
a comment on the accuracy of these results, it will be well to
compare computed results for the profiles shown in Figures 8.1
and 8.3 with experimental data. First, however, some further
discussion of the results that have been presented is in order.
First we note that the effect of increasing a is to in-
crease the level of the correlations near the surface. For a
fixed A.^ , the effect of increasing a is generally to
_L ITlcLX
-------
8-10
en
-------
8-11
CO
N
20
10
0.
0
.2
0
.02
.06
.08
CO
k_
CD
N
30
20
10
0
0
.02
.04
.06
.08
.10
T
/2
Figure 8.7.
Steady state distributions of the temperature correla-
tions for the profiles of mean velocity and temperature
given in Figure 8.2. Near the surface, A,/z = a = 1.0
for the dashed curves and for the solid curves a = 0.7,
In both cases A
1 max
= 17 meters.
-------
8-12
draw out and sharpen the maxima of the correlation distributions
near the surface. Second, the effect of increasing A-, is
X ulcLX
to increase the correlations in the region above the height
given by A /a . Both of these results are expected in view
-L fflcix
of the very strong role played by the scale A-, in the super-
equilibrium theory discussed in Section 7-
We may obtain information about the dynamics of atmospheric
turbulence from the details of the steady state computations that
are carried out to obtain the curves shown in Figures 8.4 through
8.7. Figure 8.8 shows the behavior of the shear correlation with
time at two heights for the computations shown in Figures 8.4 and
8.5 with a = 0.7 and A-, = 17 meters. It is seen that
j- max
the shear correlation, starting from a zero value, has been able
to reach a near-equilibrium value in about 20 to 40 seconds.
This characteristic time of the atmospheric profile in question
is most interesting, for it indicates that the correlations should
be able to track, in the sense of approximating the equilibrium
solutions, changes in the mean profile which are accomplished in
times large compared to 20 seconds. Since the measured mean
profiles under consideration were essentially constant over a
period of the order of thousands of seconds, one would expect
that equilibrium calculations such as those we have made should
be permissible.
Let us now turn to the effect of coupling the equations for
the mean profiles to the equations for the second-order correla-
tions. We may investigate this for the case we have been
considering by choosing an appropriate A and a, computing
_L max
out to equilibrium with the mean velocity and temperature pro-
files fixed, and then using the resulting turbulence profiles as
initial conditions for considering the coupled equations. As an
example of such a run, we consider the results shown in Figures
8.9 and 8.10. Here an equilibrium solution for the case shown
in Figure 8.2 was obtained for the case when a = 0.7 and
A, = 15 meters. This result is shown as the reference lines
in Figures 8.9 and 8.10. With these initial conditions, we then
-------
8-13
.28,
24-
.20-
-u'w'at z=5 meters
-u'w' at z=20 meters
Figure 8.8,
t, seconds
Approach of the shear correlation to its equilibrium
value. Mean profiles are given in Figure 8.2.
A, = 17 meters and a = 0.7.
_L
-------
120
100
en
v
o>
N
80
6°
40
20
Starting (measured)
profile
u
Figure 8.9.
Effect of letting u profile run free for 10 seconds
with - dp/dx = 2 millibars/100 kilometers
-------
CD
E
100-
80
100-
-u'w'
8-15
w'w'
.4
u'u'
Figure 8.10
Effect of letting u profile run free for 10
with - dp/dx = 2 millibars/100 kilometers
seconds
-------
8-16
ran with all of the equations solved in a coupled mode for a
time of the order of the characteristic time of the turbulence,
namely, 10 seconds. A uniform driving force of 2 millibars/
100 kilometers was assumed to act. The results of this computa-
tion are shown as the solid lines in Figures 8.9 and 8.10. It
is seen that for this case a very small change in the velocity
profile.caused a rather significant change in the distributions
of the turbulent correlations close to the surface, i.e., within
the first 5 meters. There was essentially no effect on the .
correlation profiles above a height of 10 meters. The results
are identical for the mean temperature profile and the profiles
of the temperature correlations although plots of these results
are not presented. Additional calculations similar to tho'se just
discussed have been carried out for somewhat different conditions,
including an attempt to add effects of surface roughness to the
model. While these studies are not yet complete, it appears that,
except for a region very close to the surface, it is possible to
estimate the turbulent structure of an atmospheric boundary layer
on the basis of an equilibrium solution, knowing the mean
velocity and temperature profiles. Near the surface, this is
also possible but, in this case, it is absolutely necessary to
have some adequate model of the effect of .surface roughness. A
method for including the effects of roughness will be developed
in Section 12 of this report.
We turn now to the calculation of the state'of turbulence
for the velocity and temperature profiles shown in Figure 8.1.
For this case of a stable atmosphere, it is found that, if one
assumes a = 0.7 as seemed appropriate for neutral atmospheres,
and if, further, one assumes A, of the order of 0,15
_1_ IQcLX
times an estimate of the shear layer thickness, one greatly
overestimates the level of turbulence and the magnitudes of the
turbulent transport correlations. The answer is, as might be
-------
8-17
expected, that the stable character of the atmosphere has greatly
reduced the scale of the turbulent motion compared to the scale
of the mean profiles. A search for values of a and An
-"-max
which gave reasonable agreement of the computed results with
experimental results yielded the values of a = 0.5 and
AI =5 meters. The equilibrium distributions computed for
ITlclX
these values are shown in Figures 8.11 and 8.12. It can be seen
that this choice of parameters gives the correct magnitude of
all the correlations. In this case, the vertical heat transfer
correlation w'T' as well as the turbulent shear and vertical
velocity fluctuations are in agreement with experimental results.
2
The temperature correlation u'T1 is underpredicted, while Tf
is overpredicted.
The mean profiles for the case of an unstable atmosphere are
given in Figure 8.3- Results of such calculations are shown in
Figures 8.13 and 8.14. Here, reasonable agreement between experi-
mental results and computations was found when a = 0.7 and
A-i was taken to be 30 meters. The results show that all
xmax
correlations,are predicted within an order of magnitude. The
predictions for w'w' , w'T' , and T' are good. That for u'w'
is of the right magnitude, but the predicted distribution is poor.
The predictions of u'u1 and v'v' are low, as they have been in
all cases. Of particular note here is the fact that the experi-
mental results show v'v' to be larger than w'w' for an un-
stable shear layer! This is difficult to understand on the basis
of the equations of motion, and there is certainly no way the
model equations can yield this result. It is conceivable that the
large values of u'u' and v'v1 that were measured might have
been due to a very slight torsional oscillation of the measuring
tower. A comment is in order concerning the distribution of
u'w' . For this unstable case, the fairing of the mean velocity
profile between the data points has a profound effect upon the
outcome of the computations in view of the very large changes in
the value of 9u/3z that can result from the different fairings.
-------
(
E
40
30
N 20
10
95
.03
.04 .05
.06
-uxw'
8-18
.08
w'w'
.16
en
i_
o>
"CD
•»
M
0
.08
V'V'
.16
Figure 8.11.
Steady state distributions of the velocity correlations
for the profiles of mean velocity and temperature
given in Figure 8.1. Near the surface, A,/z = a = 0.5
and A, =5 meters .
-------
8-19
50
c/> 40
k_
o>
E
N
30
20
10
0
0 .01
.05 .06
.04 .08 .12 .16 .20 .24
50
en 40
E
N"
30
20
10
.06 .07 .08
Figure 8.12.
Steady state distributions of the temperature correla-
tions for the profiles of mean velocity and temperature
given in Figure 8.1. Near the surface, A,/z = a = 0.5
and Al max = 5 meters-
-------
8-20
30 r
2 20
-------
8-21
.12
.16
I T /
w'T
.20
CD
N
r 10-
Figure
T'2
Steady state distributions of the temperature correla-
tions for the profiles of mean velocity and temperature
given in Figure 8.3. Near the surface, A,/z = a = 0.7
and Al max - 30 -^ —
-------
8-22
Such fine points are rather academic in view of the profound
effect of stability on the scale A, . We note here that
j_ rricix
A, is, for this unstable case, far larger in relation to
1 max ' jo
the characteristic thickness of the mean profile than it was for
the neutral or stable case.
Prom the results presented above, it is clear that one will
need a great deal of reliable and detailed experimental data on
the atmospheric boundary layer if a completely satisfactory
model is to be developed. Some of the work of developing this
model to its ultimate form can be done by considering laboratory
flows in which different gases are mixed in classical free jet
and free shear layer experiments. Some headway can be made by
considering the case of the compressible turbulent boundary
layer. In the long run, however, it'is going to be necessary to
have many good sets of data from real atmospheric boundary layers
if a high degree of confidence in the detailed prediction of any
such model is to be developed.
-------
9-1
9. DISPERSAL OF POLLUTANTS PROM A LINE SOURCE
We have already seen in Section 3 that the equations for
the dispersal of a pollutant in a turbulent medium, under the
assumptions made in this study, are decoupled from the
equations for the generation of the turbulence field itself.
If, then, one is given or can compute the turbulent structure
of the earth's boundary layer, one can use the equations given
in Section 3, namely, (3-33), (3.38), and (3.39) to compute the
mean concentration field of the species a downstream from
some source of pollutant. The uncoupled equation (3.MO) may be
solved at the same time for the variance of the concentration
2
fluctuations C1 . In this section we will demonstrate the
nature of such calculations by discussing the dispersal of
pollutants downstream of a line source for a few specially
chosen cases of turbulent motion.
Let us assume that the mean and turbulent velocity and
temperature fields are steady and of the form u = u(z),
v = 0, w = 0, T = T(z), u~X" = uju/~(z), uTfT = un^U), and
~" Q Q IK IK 1 1
T" = T' (z). These are turbulent fields similar to the equi-
librium atmospheric motions discussed in the previous section.
Let us further assume that at x = 0 and in the vicinity of
the height z = z , there exists a line source of a single
pollutant species a = p . Thus, at x = 0 we will assume
that pollutant is introduced into the velocity field u(z)
through the existence- of an initial distribution of pollutant
C_(x,z) given by
C (0,z) = C (z) (9.1)
p po
This assumption requires us to specify, since £ (T = 1, that
a
the background or atmospheric species has, at x = 0 , the
distribution
Ca = C (0,z)=l - C (z) (9.2)
o po
-------
9-2
We should note here that it is quite possible, using the equations
we have developed, to assume that at x = 0 the distribution of
u is disturbed by the introduction of pollutants. In this case,
we must solve for the mean and turbulent velocity and tempera-
ture fields downstream of x = 0 which are now of the form
u(x,z), T(x,z), u!u'(x,z), etc., before calculating the dispersal
of pollutants in this turbulent field. In what follows, in order
to keep the exposition of the method as simple as possible, we
will neglect the effect upon the background turbulent field of
the introduction of the pollutants. In this case our problem
b_e_c_omes that of solving for the mean fields C (x,z) and
C^2(x,z) by using (3.33), (3-38), (3.39) and.(3.40) when the
mean and turbulent velocity and temperature fields are known
functions of z . Before displaying the forms that (3-33) and
(3.38) through (3-^0) take for the particular problem we have
set, we will make one more simplifying assumption.
This final assumption is the usual thin layer or boundary
layer assumption; i.e., that the gradients of variables in which
we are interested are far larger in the z direction than in
the x direction. Thus, in the equations, we will neglect
derivatives with respect to x compared to those with respect
to z . Since we are dealing with the dispersal of a line source
located on the line x=0,z=z in a flow for which deriva-
tives with respect to y are zero, we may also set 3/3y =0 in
our equations. With these assumptions, our equations become
3C"
(9.3)
32C
= v°^
3
3z
W'Cp
u
3w'C'
p _
3x
"ir
9 3 / 9 \ 1 3
I , , ^ J?... f T,T f n ! \ . ;. , { i
z 3z\ p/ p3zj
o
0
a'-w'c1
.-Lrtff-pi, U
J-,, P ° ^^
, 3C'
P P0 9z
o
3w' 3C'
">v p
o 3x . 9x .
(9.4)
-------
9-3
u
3x
arn . 3C a
_ w'C' —- - w'T1 «- - —— fw'C'T'
w p 3z w x 3z 3z ^w p1
+ V
32C'T'
P
o 3x 3x
J J
(9.5)
u
3C
3x
3C
= -2w'C'
P
w'C
o
1
(9.6)
To this point, our equations are exact. In order to solve
our problem, we introduce the models discussed in Section 4
into (9.4), (9.5), and (9.6) to obtain a closed set of equations
The modeled equations are
u
u
3w'C' 3C , (
P - w ' w ' - P + -_...< ( 1 A
*\ ~ — W W fN ' ~ 1 \ t. .1 L
c)X 3Z 3z \^
T
0 P
J\ /-» f m f _ J\ /"^
P — T.T f P ' i-i- ur f T f ,P
3x ~W p 3z W T 3z
3C'2 3C" . f
. P -= OT.T 1 P » , P 4. .... J A
3w'C' ^
i * \ P I Q f ^i i
2 3 q 3z j Ax p
n
3"w'C' w'C'
• i. P Oil P
0 3z2 ° X2
3 ( 3CpTt\
• + —'-• < A q ^ • 7
0
9P t Tl 1 *^* f fp t '
OX OX
i \j P °-j P
~v~ i ^-v_ i
0 3z2 . ° X2
2 > 22 2
3C' 1 3 C' C'
P I I M P - OM P
3x """ "p 3z
(9.7)
(9.8)
o A2 ».9)
-------
In these equations, we recognize the scales A,j Ap, A,,,
and X that were used in the equations which generated the
turbulent velocity and temperature fields. In what follows we
will use the same general model form as was used in the turbulence
calculations; namely, A? = A_ = 0.1A, and
A = A // 2.5 + 0.125Re. (9.10)
i AI
but we will investigate the consequence of various choices of the
scale A . In Section 4 it was pointed out that to choose A,
the same for both the calculation of turbulent velocity and
concentration fields would, of course, be the simplest choice.
Thus, following our rule of trying to find the simplest model,
the consequence of such a choice should be studied. It is not
hard to argue that this might cause some difficulty in the
calculations. Consider, for example, the introduction of a thin
stream of pollutant, say of the order of 0.1 meters thick, into
a region of the atmosphere where the integral scale L of the
velocity fluctuations was known to be of the order of 30 meters.
Since A, for this velocity field is of the same order as L ,
we might guess that some difficulty would be experienced in the
use of so large a A, in the equations for the concentration
correlations. On the other hand, if the thickness of the pollu-
tant layer was much larger than the scale of the atmospheric
turbulence, it is difficult to see how the scale of, say, the
velocity-concentration fluctuations could exceed by very much
the scale of the background turbulence. With these preliminary
thoughts in mind, let us turn to some actual calculations to see
just how the choice.of A, affects the solution of our equations.
To do this, let us consider a special case of turbulent
diffusion. Let us consider the diffusion of a line source of
pollutant into a region of homogeneous isotropic turbulence in
a neutral atmosphere that is moving by the source at uniform
velocity, i.e., u(z) = constant.' We will consider the conse-
quences of making three choices of A, :
-------
9-5
(a) A is a constant and equal to the A, of the back-
ground turbulence into which the pollutant is dispersed.
(b) A, is proportional to the half-breadth of the
concentration profile. In the calculations presented here, A,
is taken equal to the half-breadth defined as the distance
between the z for C = 3/^Cn and the z for G = Cn /4 .
P Pmax P pmax
(c) A is proportional (equal in these calculations) to
the half-breadth of the concentration profile while the A, so
defined is less than the background A, and thereafter is
constant and equal to the background A, .
For the purposes of these expository calculations, we have
chosen a background A equal to 15 meters. The turbulent field
is assumed to pass the line source of pollutant material at 10
meters per second. The initial mean concentration of pollutant
at x = 0 is taken to be a Gaussian profile with C. (0,z ) = 1 .
The ratio of the square root of the second moment of this distri-
bution to the zeroth moment was chosen to be 1 meter. This
results in an initial half -breadth of .929 meters. It 'was
assumed that at x = 0 the background turbulent field extended
throughout the region occupied by the initial distribution of
pollutant . In addition, it was assumed that initially u.'C' = 0
2 1 P
and C' = 0 .
P
. Some results of solving for the mean concentration field
downstream of the source described above in a turbulent field
-
having ( u1 ) = ( v' ) = ( w ' ) = 0 . lu are given in
\ / \ / \ /
Figures 9.1 and 9-2. These calculations were made under the
assumption that A, was constant and equal to the A, of the
background turbulence, namely, A, = 15 meters.
Figure 9-1 exhibits the behavior of the maximum concentra-
tion of pollutants C-. downstream of the source at x = 0 .
^max
Figure 9-2 depicts the mean concentration profiles at various
distances downstream of the line source. It will be noted in
Figure 9-.! that very little reduction of concentration is appar-
ent for some ten to twenty meters downstream of the source.
-------
max
10
H
10-2
Aj=constant = 15 meters
U = 10 meters/second
12
Cn ~ x~ 2
Kmax
b = 15 meters
w = I meter/second
10
Figure 9•1•
%
x , ( meters)
10
10'
Decay of maximum pollution concentration with distance when the scale length
A-
is assumed constant
-------
x=2000 (m)
Aj=constant = 15 meters
U = 10 meters/second
u
2
•
w' = I meter/second
Figure 9.2. Concentration profiles at various distances downstream of a line source when the
scale Aj is assumed constant
-------
9-8
This is because we have assumed that no turbulent transport
w'C' existed at x = 0, and the dynamic equations for the
production of this turbulent transport are such that for a mean
concentration profile with a spread of approximately 1 meter
and for a turbulence level of roughly 1 meter/second, it requires
of the order of 1 second for the turbulent transport terms to
adjust to their proper values. In the region from 60 to 300
meters downstream of the source, the decay of the centerline
distribution is rapid and approaches a rate somewhat greater
than x"'1 . At some 420 meters downstream of the source, the
half-breadth b of the concentration profile is equal to 15
meters. For larger distances than this, the rate of decay of
the centerline concentration continues with the x l behavior.
Eventually, when the half-breadth of the profile is considerably
larger than the A, used in the computation, the rate of decay
of the centerline concentration'falls off as x~1/2. This sort
of behavior is to be expected at large distances from the source
and has been observed'in experimental studies [Ref. 35]. We
will defer a detailed discussion of this phenomenon to the next
section when we will discuss dispersion from a point source.
Here we will take up the mid-range behavior of pollutant disper-
sal following the period of initial formation of the'turbulent
exchange correlation w'C' .
To what do we ascribe the unrealistic behavior of Cn
i-'max
in the mid-range 100 < x < 300 meters where the decay of
p
is greater than x ? This behavior has been traced to
the unrealistic profiles of C that are initially generated
due to the choice of A-. = 15 meters in Eqs . (9.7), (9.8), and
(9-9) when the half-breadth of the concentration' prof lie is •
consderably less than this figure. These unrealistic profiles
are evident in Figure 9.2. Note in this figure that for x = 100
meters and for x = 336 meters, the concentration profiles have
a bump or "nose" near the region of maximum concentration. It
is the elimination of this unrealistic nose in the region between
-------
9-9
100 meters, where it is very pronounced, and 336 meters, where
it has almost been eliminated, that causes Cp to decay
faster than x"1. Once the half-breadth of the pollution profile
has exceeded 15 meters, this difficulty is eliminated and the
concentration, profiles behave properly.
A close examination of the details of the computation of
the profiles in the formative and mid-range regions of concentra-
tion decay was made. This examination showed that the cause of
the unrealistic profiles and decay behavior was the existence of
too much diffusion in the equations for the correlations in the
formative and mid-range regions of decay. The amount of this
error was dependent on the ratio of the A, 'used in the calcula-
tions to the local half-breadth of the concentration profiles.
For this reason, a second solution of this dispersal problem was
run under the assumption that the A, in the calculations was
always equal to the local half-breadth b of the concentration
profile. The results of such a computation are shown in Figures
9.3 and 9-4.
From Figure 9-3, we see immediately that after a formative
distance of the order of tens of meters, the decay of the maximum
concentration of pollutant approaches a decay rate exactly
proportional to x"1. By extrapolating this decay rate back to
Cp =1, one obtains an effective source position x of
26.5 meters which is really indicative in these calculations of
a characteristic time for the formation of turbulent transport
correlations, namely, x /u as 26.5/10 = 2.65 seconds.
Examination of Figure 9-4 shows that with the present choice
of A, no difficulty in profile shape is encountered. We do
note, however, that this choice of A, will never permit a far
— -1/2
region decay where Cp is proportional to x . Since we
wish our decay rates to ultimately approach x~1/2 as well as
behave properly in the formative and mid-ranges of decay, it
appears that the proper choice for A, is that A, be propor-
tional to the concentration half-breadth b for all cases when
-------
'Pma
10-2
x0= 26.5 meters
I I
A,= b
TT = 10 meters /second
w =1 meter/second
i i i i i i
10
Figure 9.3
10'
x, ( meters)
Decay of maximum concentration of pollutant with distance when the scale length
AQ_ is assumed equal to the half-breadth b of the concentration profile
-------
x =2000 meters
u = 10 meters/second
,2
w' = I meter/second
Concentration profiles at various distances downstream of a line source when
the scale A]_ is assumed equal to the half-breadth b of the concentration profile
-------
9-12
b is less than the A, used In the background turbulence
computations and equal to this background A, when b is
much larger than A .
At this point it would appear that it is necessary to
abandon the assumption of a universal A for all the correla-
tion equations and distinguish between the A, 's in the
equations for the various correlations. Therefore, we will
identify a different A, for the concentration correlation
equations, say, A]_ , and for the background turbulence, say
O
A]_ . In terms of these new definitions 3 we may write the
u
choice of A]_ that has been suggested previously as follows:
AT = .b for b <_ kAi and A]_ = kAn for b > kA]_ .
C v C u L/
In Figure 9-5 we show the behavior of the maximum concentra-
tion for the problem that has previously been discussed, using
the above-noted choice of A]_ when the constant of proportion-
\s
ality k between A]_ and Aj_ is one. This figure shows the
kind of decay that we desire. After a formative period of decay
which depends entirely on the initial conditions applied to the
problem, a mid-region of decay is found where Cp falls off
as x~ . At very large distances, CD approaches a decay
max
-
rate of x .We will not present here a separate figure
showing the concentration profiles for this case since, as might
be expected, the profiles are identical to those of Figure 9 • ^
for most of the range of x considered.
It may be instructive to compare the results of all three
of the calculations we have just discussed. Figure 9.6 shows
the behavior of Cp (x) for each of these calculations. The
results shown in this figure may be compared with those shown
in Figure 9.7 where we have plotted the results of pollution
dispersal computations using the three A]_ models for the case
C
of a higher relative level of turbulence, namely, u = 10 meters/
second and /u'2 =/v'2 =/w'2 = 3 meters/second. It is easy to
see that while the effect of higher turbulence level affects the
rate of spread, it does not affect the general character of the
solutions .
-------
max
I0~'
10
-2
u =10 meters/second
nax
,2
,2
,2
w = I meter/second
b=A.. = 15 meters
= 15 meters
i i i
i i
i i i i i i
I i i
I I I
10
10'
x , ( meters )
10'
Figure 9.5- Decay of maximum pollution concentration with distance when
b <_ Aj, and A^ = A]_, for b > A]_, = 15 meters
= b for
10
-------
max
u = 10 meters/second
I^BB^BHBHB
w' = I meter/second
b= A(. = 15 meters
A. = 15 meters
x, ( meters)
Comparison of rates of concentration decay for three choices of A-,
i
M
J=-
-------
Cp
max
10'
Figure 9-7-
u = 10 meters/second
u
,2
w1 = 3 meters/second
0 A = 15 meters (Run* 4)
v A = half -breadth (Run* 5)
* A=half-breadth < 15 meters
(Run *6)
= 15 meters
ICT 10'
x , meters
Comparison of rates of concentration decay for three choices of
MD
I
VJ1
-------
9-l6
Before going on to present a calculation of the dispersal
of pollutants for a more realistic atmospheric condition, a
presentation of the nature of the solutions for the second-order
• 2
correlations C'w' and C' is in order. Figure 9.8 exhibits
plots of both C and C'w' at 100 meters from the source for
P P
the calculations presented in Figure 9-3- For this same condi-
;ti_on, the variance of the fluctuation in pollutant mass fraction
2
CT is shown in Figure 9.9- Note that the actual variance is a
maximum near where the transport of pollutant is a maximum. If
one computes the ratio of / C'2 to C , it is interesting to
not_e that this ratio has a value of 0.84? at the point where
2
C' is maximum. As the mean concentration falls off towards
the outer edges of the sheet of pollutants, / C^,2 falls off
less rapidly than C so that the ratio of / C ' /C tends to
approach large values. This is what is to be expected when
there is an occasional "blob" of pollutant in an otherwise
unpolluted background. This behavior of /C'2/C is also
sr C*
shown in Figure 9.9-
Having studied the nature of diffusion calculations by
second-order modeling for a very simple case, let us now turn to
the calculation of pollutant dispersal from a line source under
more realistic conditions. Let us consider an atmospheric
boundary layer having the mean velocity and mean temperature
distributions shown in Figure 9-10. We have assumed that a
stabilized region exists in the region between 88 and 140 m
above the earth's surface. Also shown in Figure 9.10 is an
assumed distribution of the background turbulence scale A]_ .
In order to calculate the dispersal of pollutants in this layer,
we first compute the field of atmospheric turbulence produced by
this combination of conditions. To do this, we use the methods
described in Section 8 to obtain profiles of all the turbulent
correlations necessary for the solution of Eqs.(9.3)> (9.7),
(9.8), and (9.9). The results of such calculations are shown in
Figures 9.11 and 9.12. With these results in hand, Eqs. (9.3),
-------
25r-
9-17
0.05
0.10
0.15
0.20
0.25
Figure 9.8
0.01
0.02
0.03
0.04
0.05
0.06
w Cp , ( meters /second )
Distributions of C and w'C' as functions of
at x = 100 meters ior dispersal problem shown in
Figure 9-3
-------
20r
9-18
CO
N
C/p2/Cp
0
0
.1
.4
.5
2 ~~
Figure 9.9. Distributions of C^ and / C£ /Cp as functions
of z at x = 100 meters for the dispersal problem
shown in Figure 9-3
-------
9-19
240r
24 6 8 10 12
u, (meters/second)
240
200
160
v>
u.
0>
1 120
Nl
80
40
0.
Neutral
Stable
region
Neutral
4 8 12
A|t .(meters)
16
0
Figure 9.10.
Assumed mean velocity and temperature profiles
for an atmospheric boundary layer. Also shown
is an assumed profile of the background turbulence
scale An
-------
240
9-20
240
-u'w'fmVs")
.4 .8 1.2 1.6 2.0 2.4
uV(m2/s2)
.4 .8
v'v'(m2/s2)
Figure 9.11. Velocity correlations computed for assumed atmospheric
layer shown in Figure 9.8
-------
9-21
z, meters
240r
200
I60V-
120
80,
40
-.04 0 .08 .12 .16 .20
u7T7(m-°K/s)
-.003 -.002 -.001 0 .001 .002 .003
w7T7(m-°K/s)
24 Or
.01 .02 .03 .04 .05 .06
f7r(°K2)
Figure 9.12. Temperature correlations computed for assumed atmospheric
layer shown in Figure 9-8
-------
9-22
concentration of pollutant C (x,z) downstream of a line source
(9.7) and (9.8) may be solved simultaneously for the mean
concentration of pollutant C
located at x = 0 and z = z
o
As in our previous calculations, we assume an initial
source of pollutants that has a Gaussian distribution with
C" (0,z ) = 1 . The half-breadth of this distribution is again
Pmax °
assumed to be .929 meters.
In view of our experience with the spread of pollutant into
uniform turbulence where, for pollutant spreads less than Aj_ ,
we found that A-j_ must be chosen proportional to the spread b
\^
rather than equal to A^_ , we may anticipate a difficulty with
the computation we are contemplating. Since the spread above
and below the maximum concentration will not remain equal because
the problem is no longer symmetric, we must have a way of choosing
A} differently on the two sides of the concentration profile.
We have'selected a method for doing this which appears to be
appropriate for the problem under discussion here, as well as for
more complex pollutant concentration distributions. Consider for
the moment the concentration distribution shown in Figure 9-13-
This distribution has five extrema: 2 maxima, 1 minimum, and 2
zeroes. The range of z is broken into ^ regions which are
defined by the ranges of z separating these extreme values. In
each region, AQ_ is assumed to be constant and equal in magni-
O
tu.de to the distance between where the change in C is 0.25
and where the change in C is 0.75 times the change between the
two extreme values that define the region under consideration.
The A-, so defined are identical to the half breadth b
-"-c
when this method is used for the symmetric dispersals we have
discussed previously, and we will use it for the computation under
consideration here.
Figure 9-1^ shows the behavior of solutions for the mean
concentration profiles downstream of the line source we have
discussed when the pollutants are injected into the atmospheric
boundary layer described in detail by Figures 9.10 through 9.12.
Note that the spread is initially almost symmetric, but at
-------
9-23
c,(z)
Figure 9.13-
Hypothetical concentration distribution showing
regions of constant An and values of Ai
xc c
to be used in each region
-------
I20r
Stable
Start of stable region
Figure 9
Computation of the dispersal of pollutant downstream of a line source
below a layer of stable air
-------
9-25
distances of the order of 500 meters from the source which is at
a height of 50 meters, the ground has played an important role
in increasing the levels of concentration. Between 500 and 2000
meters from the source, it is easy to see that the region of
stable air has formed an effective lid on the dispersal of
pollutants, and the distribution of pollutants has become almost
uniform below this stable region at 2000 meters from the source.
In order to see a significant decrease in this level of pollutant
concentration, the computations would have to be run to enormous
distances downstream from the source because of the-very slow
leakage of pollutants through the stabilized region. Indeed, for
strongly stabilized regions, calculations such as those described
here show that the only real hope for reduction of pollution
levels is a change in atmospheric conditions.
It.is hoped that this brief description of the calculation
of pollutant dispersal by invariant modeling of the second-order
correlation'equations has served to demonstrate the power of this
new method. While there are several improvements that must be
added before this second-order closure method can be considered
complete, there is no question but what the model will enable one
to make computations of pollutant'dispersal for a number of cases
for which-no reliable methods have been available in the past.
Before discussing what some of these improvements in the- overall
model might be, we'will complete our discussion of pollutant
dispersal by considering the problem of dispersal from a point
source of pollutant-.
-------
10-1
10. DISPERSAL OF POLLUTANTS FROM A POINT SOURCE
In this section we will consider the dispersal of pollutants
from a point source. We will assume that the background atmos-
phere is given. For atmospheric boundary layers, we may again
use solutions of the form used in the previous section, namely,
u f. u(z), v = 0, w = 0, u.'u,' = u.'u'(z), u!T' = u.'T'(z), and
Q O 1K1K 1 1
T' = T' (z). If we consider a point source of pollution for
such flows, the distribution of pollutant downstream of the
source is steady (if the source is steady) and is of the form
C~ = C (x,y,z). If we again make the boundary layer or, in this
case, the thin plume assumption, we may neglect derivatives of. quan-
tities with respect to x in comparison with the same deriva-
tives with respect to y and z . In this case the equations
which govern our problem are
92C
9z'
_ ^- (C'v') - ~ (C'w1)
9y P 9z P
(10.1)
u.
9x
9C 3
vtv' _£ + f-
9y 9y
+ 1)
90
9y
£1]
(10.2)
-------
10-2
u
3C'w'
P
,.3x
3C
- w'w'
_
A ) q
9z
3C'v'
£r I A2CQ a^
3z
3C'v'
P
3y
C'W1
P
3 C'w
P
C'w1
- 2-u
(10.3)
u
3C'T'
P
3x
5— - W'C' ^~
3z p 3z
3C
b
3C'T'
P
3z
|32C'T' 32C'T'
.P + V-
o
C'T'
-2—-
(10.4)
u
3C
I
3x
'2
3C
3C
3C
3C
P ? ? ?
3 C' 3 C'
- 2-u
3z
0 X2
c
(10.5)
-------
10-3
A digital computer program has been developed which can
simultaneously solve (10.1), (10.2), (10.3), and (10.4) for the
mean concentration C (x,y,z) downstream of a steady source when
the background mean and turbulence fields are known. In the
process, of course, the distributions C'v'(x,y,z), C'w'(x,y,z)
and C'T'(x,y,z) are obtained. With these results in hand,
P 2
(10.5) may be solved for the distribution CT (x,y,z).
As is the case for the other digital programs mentioned in
this report, an outline of the numerical techniques used to
construct the program is given in Appendix A which is dedicated
to this subject. Here we must note that the program required to
accomplish this task is a very large one and, in its development
and check out, an uncommonly large number of "bugs" emerged.
Most of these bugs have been eliminated, but the program is not
yet operational in all possible modes and options. In particular,
the existence of the cross-derivative terms in the equations for
the second-order velocity-concentration correlations still leads
to certain difficulties in obtaining smooth solutions.
We may demonstrate the nature of the solutions that are
generated by second-order modeling at this time, however, by the
choice of a special turbulence model. It is clear from an exam-
ination of (10.2) and (10.3) that if we choose a turbulence model
in which Ap = 0.1A-, and A-, = - Ap , the cross-derivative
c -"-c Jc c
terms in (10.2) and (10.3) are eliminated. This assumption, for
the two-dimensional concentration layers studied in the previous
section, has an effect almost analogous to making ^2 and ^3
somewhat smaller compared to A]_ than the usual value of 0.1.
For this reason, it is felt that the solutions that will be
discussed below are representative solutions and should exhibit
the nature of pollution dispersal solutions generated by invariant
modeling techniques. At present, every effort is being made to
clear up the difficulty that is occasioned by the retention of the
cross-derivative terms in the equations. It is hoped that this
difficulty will be resolved in the very near future.
-------
10-4
To examine the nature of these solutions and to enable one
to compare the results of computations of both line and point
source dispersal, we shall compute the dispersal of an initially
Gaussian concentration of pollutants into the same homogeneous
and isotropic turbulent background atmosphere that was used in
the calculations presented in Figure 9-5•_ For this case,
- / 2 / 2 / 2
u = 10 meters/second, / u1 = / v' = / w' =1 meter/second,
and A]_ = 15 meters . We will again assume that no transport
correlations C'v' or C'w' exist at x = 0 . The initial
Gaussian profile used has the same half-breadth as before, namely,
0.929 meters .
In Figure 10.1 we show the behavior of the maximum concentra-
tion Cn as a function of distance from the source of
Pmax
pollutants. For the case shown here, in order for the results to
be comparable to one of the line source dispersals presented in
Figure 9-5, we have again chosen A]_ = b for A]_ <_ AT, and
C C w
AI = A]_ for b > A], . It is seen from Figure 10.1 that after
C tx G
an initial formative stage of essentially the same duration as that
_2
for the line source, the maximum concentration drops off as x
In this mid-region of decay, the half-breadth of the plume is less
than A]_ . At a point x = 429 meters from the source, the half-
breadth becomes equal to AT . Shortly thereafter, the fall off
- -1
in Cp becomes proportional to x . This is analogous to the
behavior found for the line source when CD as initially propor-
-1 max_l/2
tional to x and finally proportional to x .In both cases,
this behavior is due to the fact that the half-breadth of the sheet
1/2
or plume in question grows initially as x and finally as x
In Figure 10.2, we show this behavior for both the line source and
the point source by plotting the effective spread a of the sheet
z
or plume of pollutants as a function of x for the same background
turbulent conditions used for Figure 10.1. Here a is defined as
Z
/-°°(Z - zQ)2C (x,0,z) dz
'+°°C (x,0,z) dz
-m f
-------
10'
10
-I
X
o
E
Q.
10
10
-2
10
-3
10
-2
u =10 meters/sec
-------
£ 10
O)
N
b
Po nt source
10'
Figure 10.2.
10'
10
x, meters
Comparison of the spread of pollutants from point and line sources for
= O.lu
and A]_ = 15 meters
/u'2 = /v' 2 = / W 2 =
O
I
-------
10-7
The results shown in Figure 10.2 show the very great similarity
between the spreading rates for line and point sources that is
predicted by the present model of turbulent diffusion.
In Figure 10.3 are shown the central mean concentration
profiles for this case of point source dispersal. Figure 10.4
shows the'behavior on a traverse through the center of the plume.
of C (z) and C'w'(z) at 100 meters from the source.
We may compare the results of computations such as those
just presented with experimental data. In Figure 10.5, we have
plotted the behavior of the spread of pollutants downstream from
the source we have described for three cases of atmospheric
turbulence: aQu = 2, aQu = 1, and aau = 0.2. Here
2— / 2 — / 2 —
aQ = / u' /u = / v' /u = / w' /u . The results of these computa-
u
tions are compared in Figure 10.5 with the experimental results
of Fuquay, Simpson, and Hinds [Ref. 35]. It is seen that the
rates of spread predicted by the computations are within a factor
of two of the experimental results in the mid-range of decay. The
departure of the computed results from experiment for breadths
greater than 15 meters might be due to two causes. It might be
that the background scale of turbulence A]_, was greater than
15 meters or it might be, if A-j_ is actually of the order of
15 meters, that one must assume that the proper A]_ to use in
O
the computations should be equal to the half-breadth of the plume
until A]_ is quite a bit larger than A], before setting A]_
c ^ c
equal to a constant times A}_ . In view of our calculations of
U
atmospheric turbulence in Section 83 it would appear that Ai ,
L>
for the experiments in questions was not greater than 15 meters.
We must therefore conclude that the second' reason for the depar-
ture is the correct one. In future studies it will be necessary
to pin down the precise point at which-the scale A]_ departs
O
from the half-breadth b and becomes proportional to the scale
of the background turbulence and to determine the magnitude of
this constant of proportionality.
-------
50r
40
S 30
0>
20
10
0
= IOOm
= 43.3 m
x = 0 m
10"
10
-2
10
-I
Cn
rmax
Figure 10.3. Central mean concentration profiles at various distances downstream of a
10'
o
I
oo
source when
u'
2 _
V'
2 _
w'2 ~ 0. lu and A]_, = 15 meters
-------
10-9
0>
-------
io
10'
Q)
0)
E
10'
IOV
U = A/V' 2 > 1.0 (m/s)
<0.2 (m/s)
Measured data
(From Ref. 35)
Extrapolated data
(From Ref. 35)
Computed
10'
10"
t, seconds
ICT
10'
o
I
Figure 10.5-
Comparison of measured and calculated lateral spreads
a as a function of time. Data are from Ref. 35-
-------
11-1 .
11. CRITIQUE OP ATMOSPHERIC TURBULENCE AND POLLUTION DISPERSAL
CALCULATIONS
In the previous three sections, we have outlined the charac-
teristics of a new method of handling the computation of the
structure of atmospheric turbulence near the surface of the earth
and the dispersal of a passive pollutant in such a turbulent field,
It is clear that in order to carry out such calculations, it is
necessary to know the magnitude of certain scalar lengths (the
A, 's) that enter into the formulation of. the equations for the
second-order correlations necessary to calculate the mean fields
of velocity, temperature, and species•concentration. It was
pointed out in Section 6 that it is possible to obtain equations
for these lengths that can be solved simultaneously with the sets
of equations we have solved in the previous three sections so as
to have a closure of the turbulent problem. The construction of
such equations and their incorporation into our computational
programs is, of-course, an item of high priority in the further
development of our method. For the time being, however, we will
continue to study the nature of turbulent flows by invariant
second-order modeling under the assumption that we can find suit-
able scales associated with the distributions of the mean
variables u , T , and C that will'enable us to make'meaning-
ful calculations of turbulent flows that are of importance in
engineering and which have heretofore been intractable. In all
the problems, except one, that we have discussed in the previous
sections, there have always been lengths associated with the mean
profiles that were available on which to base the A, 's that are'
required for the computations. For a free jet-and a two-dimen-
sional shear layer, there were the half-breadths of the layers in
question. For the classical boundary layer, the boundary layer
thickness was available. Near a solid surface it was shown that
the A, 's involved must behave in a particular manner, namely,
proportional to the'distance from the surface. For the spread of
a pollutant, we always have available a local breadth associated
-------
11-2
with the mean concentration profile. Only in the case of the
atmospheric surface layer, as studied in Section 8, was there no
scale available on which to base a calculation as one moved away
from the surface to a region where A, was no longer proportional
to z , the distance from the earth's surface. Surely there must
be some appropriate length. It is not the Monin-Obukhov length,
for this length does not determine a solution but merely correlates
solutions once the magnitudes of heat transfer and shear are known.
The answer to this question is not to be found in the equations we
have studied to date. It must be sought in the more complete
equations that govern the motion of the atmosphere on our rotating
earth. Indeed, the missing outer scale required for the surface
layer calculations discussed in Section 8 must be determined by
the balance of centrifugal (Coriolis) and pressure forces that
determine the thickness of the Ekman layer. For planetary bound-
ary layer calculations, it is the'scale of the Ekman layer that
ultimately determines, in conjunction with surface roughness, the
actual values of surface shear and'heat transfer under any given
conditions of surface heating or cooling. The lowest surface
layer can then exhibit Monin-Obukhov similarity•based on the
length .determined from the relative magnitudes of w'T' and u'w1
near the surface.' This behavior is very similar to the behavior
of the classical turbulent boundary layer on a flat plate. In
this .case,, there is a similarity near the wall (the law of the
wall) which is based on the level of shear .near the wall (see
Figure 5.15). The actual magnitude of this shear is not determined
until the whole boundary.layer thickness 6 is known even though
this thickness may be orders of magnitude larger than the region
of the boundary layer where the law of .the wall is valid.
In view of the above discussion, we shall take up in the next
section the calculation of the entire Ekman layer by the methods
discussed in this report. Before going on to this, however, we
should discuss what might be done to improve the accuracy of the
basic model with which the calculations reported here have been
carried out.
-------
11-3
It is clear from the way In which the model was constructed
that the model was carefully tailored to give good results for
the mean and turbulent velocity fields that have'been carefully
measured for certain classical shear flows. Due to the paucity
of experimental data on the behavior of temperature fluctuations
in turbulent shear flows, it was assumed that the model valid for
velocity fields could, for the purposes of evaluating the poten-
tial of the method, be taken-over directly.to the calculation of
temperature-velocity and temperature-temperature correlations.
An indication that this is not true is the inability of the
atmospheric turbulence calculations presented in Section 8 to
2
properly account for the magnitude of T" . Also, the very high
value of critical Richardson number obtained in the'discussion of
superequilibrium theory in Section 7 is an indication that the
modeling may not be quite right. It has been found that a slight
reduction'in'the magnitude of the tendency-towards-isotropy term
u.'T" equation compared to this tendency in the u.'u,' equation
1 1 1C
has a pronounced lowering effect upon the critical Richardson
number found by superequilibrium calculation.
The high value of critical Richardson number obtained by
superequilibrium analysis cannot be taken too seriously by itself,
for this analysis neglects the"effect of diffusion which actually
does play an important role in determining the balance .of -produc-
tion, and loss of a given correlation at a point in actual computa-
tions . All that we can say at this time'is that a further study
of the modeling of the temperature correlation equations needs to
be made. How might this be accomplished? If, indeed, the equations
we have been using contain an adequate description of atmospheric
turbulence, then one should be able to obtain from these equations
the Monin-Obukhov similarity theory. If'this can be done, the
agreement of the computed Monin-Obukhov functions with the"functions
derived from actual measurements might be used to select a more
appropriate model.
-------
My colleague, Dr. W. Stephen Lewellen, has recently
succeeded in deriving the whole of Monin-Obukhov similarity
theory and the Monin-Obukhov functions from the model equations
given in the previous sections [see Ref. 36]. A complete discus-
sion of this development is beyond the scope of this report; the
reader is referred to Ref. 36 for details. It is sufficient to
say here that this new development opens up a very fruitful
approach to the problem of improving the modeling of the temper-
ature correlation equations.
Let us turn now to a discussion of the calculation of the
Ekman layer and the dispersal of passive pollutants in such .a
layer.
-------
12-1
12. INVARIANT MODELING OP THE EKMAN LAYER
In-this section we will discuss the calculation of the •
complete planetary boundary layer and the dispersal of passive
pollutants in such a layer by the technique of invariant
modeling. We will start by writing, as is customary, the
equations of motion for the planetary atmosphere with respect
to a Cartesian coordinate system at rest relative to a point on
the earth's surface. The x axis is taken to be the direction
of constant latitude and positive in an easterly direction, the
y- axis is taken positive in the direction of north, and the z
axis positive in a direction normal to a geopotential surface at
the surface of the earth. In this coordinate system (x, y, z),
let the velocities be (u, v, w). The equations of motion for
velocities small compared to the surface rotational speed of-the
earth are then [see, for example, Ref. 37],
P = - - + 2pft (v simjr - w cos
+ I- T + I- T +|-T (12.2)
8x xx 8y xy 3z xz \ c. i
If - 2i5fl u sln *
~ Dw
P Dt = " " PS + 2p" U COS
-------
12-2
, DT = 9£+ u l£_ + 3_(k 2T\ (12 }
r\ "R4- ^-*--T^v CivlSiv \ -1-*-•->/
PL/0 0 I/ .1 O A . 0 A , \ O A .
J J \ J
In these equations, ft is the angular velocity of the earth, <|>
is the latitude, and
/ 9u. 9u, \ 9uk
T1J =^ + 4j t6UX3^ (12'6)
In order to make contact with the studies we have made in the
previous sections, we will not use these equations directly but
will assume that the Boussinesq approximation that the Prandtl
and Schmidt numbers are one holds, and that variations of the
coefficients of viscosity will not seriously affect any computa-
tion we will make so that (12.1) through (12.6) may be approxi-
mated by
2ft (v sin - w cos
p oX
il = - !- |E - 2B u
(12.8)
(12.9)
-------
12-3
.
Dt
P = - ^ T (12.12)
o
where we recall (see definition (2.10)) that
T = T - T (12.13)
o
These equations are identical to the basic equations used
earlier in this report (see Section 3, Eqs . (3-4), (3.5), (3-7),
(3.8)) to describe the atmospheric motion in the surface layer
except for the Coriolis terms which contain' the earth's rotation
fi . This being the case, we may very easily derive the terms we
must add to our second-order closure model of turbulence to
account for the effect of these terms .
It is clear that the terms that are added to the three mean
momentum equations by the inclusion of Coriolis forces are :
In the equation for u , a term on the ' right-hand side equal to
2ft(v sin (j) - w cos )
In the equation for v , a term on the right-hand side equal to
sin (j)
In the equation for w , a term on the right-hand side equal to
2S7u cos (j>
To obtain the term to be added to the u'u1 equation, it will
be recalled that the equation for u'u' was obtained from the
equation for u' by multiplying this equation by 2u' and
averaging the result. It is clear, then, that the term to be
added to the right-hand side of the u'u' equation is
sin
-------
12-4
Likewise the terms to be added to the right-hand sides of the
v'vr and w'w' equations are, respectively,
-4flu'vf sin
'w' cos
To obtain the term to be .added to the right-hand side of the
equation for u'w', recall that this equation was obtained by
multiplying the equation for u1 by w' and the equation for
wf by u1 and adding the resulting equations. Following this
procedure with the new terms that must be included, we find that
we must include a term on the right-hand side of the u'w1
equation equal to
/ _ 2 2
2ft (v'w' sin 4) - w' cos <|> + uf cos
Likewise for the equation for u' v1 , we have
2ft ( v' sin - v'w' cos - u' sin 4>J
and for the v'w' equation
2ft(-w'u' sin + u'v' cos )
The equation for T is unchanged since no Coriolis terms
appear in this equation. However, in the equations for u'T'
v'T' , and w'T' , there will be terms that result from the
earth's rotation. These terms are:
In the equation for u'T'
2ft(v'T' sin <{) - w'T' cos <)>)
In the equation for v'T'
-2flu'T' sin
In the equation for w'T'
2ftu'T' cos
Since the equation for T is unchanged, there will be no
2
ne\\r terms appearing in the equation for T' as a result of the
-------
12-5
It is evident from an inspection of all the terms that must
be added'to our system of equations to account for the rotation
of the earth that only second-order correlations themselves appear,
Therefore, the inclusion of these terms requires no additional
modeling.
If these additional terms are added to the modeled equations
derived in Section k, we obtain a set of equations suitable for
computing the planetary boundary layer: These equations are
9u dv_ + 9w
9x + 9y + 9z
~- + 2v sin $ - w cos
(P- ?- ?-
(—. r\ C- t\ C-
u + 9u + 9u
x2 9y2 9z2
9
(u'u') - (u'V) - (u'W) (12.15)
- - - sin
Dt p 9
IP- ?- ?-
r\ £. r\t- r\t_
9v9v9v
o" ' n" ~ T
9x2 9y2 9z2
- ^ (v'u') - |y (v'v') - |^ (v'W) (12.16)
cos
- f- (w'u') - f-' (w'v') - I- (w'w1) (12.17)
-------
= ,
Dt ° ax2 3y2
a_
ay
12-6
(12.18)
Dt
sln (j, _ uiwi cos
_ 2u!u'
J
A0 q 2
9u!u
2 - A
AT-lu'u' - 3
V
z—z u'u' - 2v
o „ 2 o
u'u1
(12.19)
Dv' v'
Dt
sin
-------
cos
1
Q
3xj
9w'u!
j 3w'w
3x.
J
3u!w
+ 2 •£-
3z
_ _^_
Du'v'
Dt
(2 2
v1 sin - u' sin <(>
3v
- u ' u ! TT— - v ' u '.
3u
3x.
J
12-7
+ v
0 3x2
o
- 2v
W ' W '
(12.21)
3u!v
3u!u'
J
3x.
J
It
__
0 3x2
_ 2v
(12.22)
-------
12-8
1J U
/ p~ 2 \
= 2n(v'w' sin (J> - w' cos <\> + u1 cos cf>J
r - + - ir
9x
, 9w'u! 9u'u'
AO Q I ^^~i + ^~^-
w
9x.
J
3u]w'
^j
v
_
0 9x2
J
- 2v
(12.23)
DV
_ _
Dt = 2fl(-wfu! sin <(> + u'v' cos ((>)
- v'u
1111
9w
9x
j
T" V'T>
TO
9v'u! 9w'u^.
:t- v^z
ay
3x
J
, iL. ^T . 2v ^f
o ...2 o ,2
(12.24)
-------
12-9
Dt
= 2«(v'T' sin <|> - w'T' cos
^T1 -T-- -•• ^i]
- u'u! -— • T'u'. -1^-
j 3x. j 3x.
A
8ujT> , 3u'T'
3x.
J
3
31
9u'
Ai
u'T1
v ^^r u'T' - 2v
0 3x2
J
u'T
(12.25)
_
2flu'T' sin * - v'uj
-
- T-uj
, I rp
31
^
3v'T'
3x
J
3u!T'
J
A.
v'T'
t rp t
o
v'T' - 2v
v'T
0 X2
At
(12.26)
Dw'T'
Dt
cos
- w'u'
J
- T'u'
a T.T
3x.
J
A n
3tq
- —3— w'T'
1t
+ v
w'T' - 2v
w'T1
0 A2
At
(12.27)
-------
12-10
Dt
_
3X
+ v
,2m, 2
2v
,,2
0
(12.28)
At this point it is appropriate that we present the equations
for the dispersal of a passive pollutant species C . Since the
equation for C , namely, (3-6), will not be changed by the
a _
inclusion of Coriolis forces, the equations for C ,
2
'
and
C "T ' will be those we have previously used. We must, however,
add new terms to the equations for u'C ' , v'C' , and w'C' . The
06 Ob Ok
terms to be added to the right-hand sides of these equations are,
respectively,
2fi(v'C' sin
- w'C' cos
-2nu'C' sin
a
and
2f2u'C' cos (f>
a
With these additional terms added to our previous equations,
namely, (4.24), (4.25), and (4.26), we obtain the following model
equations for the dispersal of a passive pollutant:
DC"
01
Dt
32CQ
'° ^
3x.
8x
ay
_
(v'C'.) -
h (W'CA}
a 6 ex
(12.29)
Du'C'- _
___£ = 2n(v,c, sln
_ w,c, cos
8C
,u, _ _ u,c,
3
3x.
J
3
/8u'ci '"j0^
L^ql3xj ' 3x JJ
0
"u'C1 u'C'
q
H
(12'30)
-------
Dv'Ca
Dt"
sin
u'.v'
3C
__
3x
j a 3x.
J
3x
A<
3v'C'
a
3uJC.
3y"
a
- An q
3y ^
u!C?
a
A-
v
a
DW-C;
Dt
'w'
n'w
uw
3C(
3x~
T' + S _ n i rn i
jLa 3x. TQ a
3w'C' 3u!C'
a J a
9x, "2C^
fe(A3c^|3rujCa
Dt
= -u'.C
8CaT'
12-11
o
^lr I P I
3 v'C
V
a
0 3x2
- 2v
i r.'
a
v'C
0 A2
(12.3D
2''
3w'C
+ v
w'C
0 3x2
- 2v
(12.32)
32C'T' C'T"
+ v —^— - 2v -V
0 9x2 ° X2
(12.33)
-------
12-12
DC'2 3C a 9C'2
5T-- -2uj°i K7 + 1x7
J J
2
1
(12.34)
n — = o -
0 3x2 ° X2
. J c
We may now turn to the problem of computing the planetary
boundary layer or Ekman layer by means of the equations we have
set forth above. The simplest case we may compute in order to
demonstrate the method is that of a planetary boundary layer
such that all derivatives with respect to x and y vanish,
and the atmospheric quantities are functions of z and t alone
For this case, in view of the continuity equation (12.14), we may
take w = 0 .
The equations for this special case of motion in the
planetary boundary layer are then
2-
. - 1- + v(M Bln *> + V 2-! - *. (12.35)
-
(12.36,
cos do - P^ £ (12.37)
3T „ 3 T 3w'T' no
•TT- = v —5- - v- (12.
at O *.„£ aZ
a Z
p _
P = - ^ T (12.39)
o
-------
12-13
au'u'
at
= 4n(u'v! sin - u'w' cos
- 2u'w' ^ + ^—
au'u'
'u1 -
_2V
u'u'
0 x?
(12.40)
sin
av'v1
!tq 9z
^ "'"' -3
9z
-2v
V
(12.41)
3w'w'
at
cos
w'T- + ~
8w'w
A-
w' w
v
0
_2v
ww
(12.42)
3u'v'
at
/ 2 2 \
= 2^ ( v1 sin <(> - v'w' cos - u' sin
-------
12-14
p O \ __«__ ^11
sin - w? cos <)> + uf cos c|))- w'w' -j—
T
f rp «
- 2A2 + A, q
3z \ 2t 3t/
Al,
v
0 X2
At
(12.44)
2Q(-u'W sin <(> + u'V cos $ ) - w'w1 ~ + - v'T'
o Z
2A9 + A, q
^•h J+- I
3v
—i—T j. 3 v'W T ,
. v'w' + v —= 2v^
An O ^_2 O
(12.45)
| rp
-1
^_^_^ _ ___^__ __^__
= 2fl(v'T* sin (j> - w'T' cos ) - u'w1 |- - T'w' -
ou o Z o Z
+ ^r A?. q
3u'T'
3z
- u'T' + v
Alt ° 3z2
(12.46)
-------
12-15
= -2flu'T' sin - v'w' - T'w1
3t oZ 3z
hl*o.q
_
Al ° 3z2 ° A?
•
—m__rT_-_r H T O*
cos 4> - w'w' 3— + ^—
d 2 X
o
+V __2v (12.48)
Alt ° 8z2 ° \l
arid
_
3t ~ ~ 3z 3z
22 2
pi^m | <- rp T ^
+ v £-4 -- 2v ^-5- (12.49)
o ~ ^ ,
9 z .A ,
We will discuss solutions for the complete atmospheric
boundary layer using (12.35) through (12.49) by first considering
the case of a steady neutrally stable atmospheric motion. In
order to obtain a solution, it is necessary for one to specify
the driving forces (1/p )3p/3x and (1/p )3p/3y which, if the
flow is to be steady, must be independent of time and must be
functions of z alone. Let us assume for simplicity, although
-------
12-16
the assumption is not necessary, that these driving forces are
constant and independent of height. In this case, if the surface
of the earth were frictionless, a geostrophic wind would exist
that would be constant' and independent of height. The components
of this geostrophic wind would be
ug ~ 2R sin <{> pQ 9y p f 9y uo
v =
I_ lP_ = 1 IE. - v d? 51)
p 3X p f 9X " O V-L^O-L;
sin
-------
12-17
Let us consider a specific case. We wish to compute the
planetary boundary layer for a case when
Vz> •• p;r If • vg • °
u (z) = - —jr TT^- = u =10 meters/second
° P0f 3y S
We will assume that the surface at z = 0 is smooth and, further,
-4
that 2°, sin = 10 rad/sec. In order to carry out a solution
of our equations, we must specify the length scale Aj_, . We will
assume here that Aj, may be taken equal to 0.15 times the height
at which the wind speed / u2 + v2 becomes equal to the geo-
strophic value / u2 + v2 = |u |. This assumption is very much
o o o
akin to the assumption made for an ordinary boundary layer in
Section 5, namely, A = 0.156 QQ .
-L . y y
A calculation of the character of the Ekman layer on an
absolutely smooth earth under the condition u = 10 m/sec and
S
v . = 0 carried out as described above is shown in Figures 12.1
o _ __
through 12.5. Figures 12.1 and'12.2 show the computed u and v
profiles, respectively. Figure 12.3 depicts the Ekman spiral for
this case. It will be noted that the wind speed first reaches
the geostrophic at an altitude of approximately 150 meters. It
will also be noted, from Figure 12.2, that the total thickness of
the Ekman layer is approximately 1500 meters. For this case, the
surface angle is computed to be 13.8°. Figures 12.4 and 12.5
show the values of all the second-order velocity correlations.
One notes in Figures 12.4 and 12.5 that the distributions of u'u1
and u'w' are not smooth but have a "bump" between a height of
50 and 150 meters. Whether this bump is a physical reality .or due
to some bug in the computer program is unknown at present. As we
shall see, this "bump" did not appear in the calculations carried
out for a rough earth. Of some interest, insofar as the dispersal
of pollutants in the atmosphere is concerned, is the fact that the
root-mean-square values of the vertical and lateral velocity
fluctuations run about 2% of the geostrophic wind velocity for
-------
12-18
0)
«*-
0)
E
N
1600
1400
1200
1000
800
600
400
200
0
0
8
10
12
u, meters/sec
Figure 12.1
Profile of mean velocity u for steady state
planetary boundary layer when u = 10 m/s,
v = 0 and z = 0
g r
-------
12-19
z, meters
I600r
1400
120
-1.0 -.5
0
v, meters/sec
1.5 2.0
Figure 12.2.
Profile of mean velocity v for steady state
planetary boundary layer when UK = 10 m/s,
v = 0 and z = 0
-------
(v-Vg)(m/s)
4
zr=0.0 m
300
-4
u-Ug, meters/sec
Figure 12.3. Ekman spiral for a smooth earth when u = 10 m/s, v =
500
g
0, and z =0
ro
i
ro
o
-------
12-21
a)
N
1600
1400
1200
1000
800
600-
400-
200-
'0 .04 .08 .12 .16 .20
LrV(m2/s2)
Figure 12.4. Distributions of u'u', v'v1, and w'w' for
the planetary boundary layer on a smooth earth
for u = 10 m/s and v =0; z - 0
-------
12-22
.12
.08 .04
0 .04 .08
"u'wMm'Vs2)
.12
Figure 12.5.
Distributions of u'w', v'w', and u'v' for
the planetary boundary layer on a smooth earth
for u = 10 m/s and v =0; z =0
o o ^
-------
12-23
some 200 meters above the immediate surface layer where these
velocities are roughly constant for some 20 meters above the
surface at about 3% of the geostrophic wind.
The angle of the flow at the surface is related to the
amount of surface friction. For this particular calculation,
the nondimensional friction velocity defined by
u*/ug = Vpo(0)Sg
is computed to be 0.02. Both this value and the surface angle
of 13.8° are considerably lower than one usually finds. This is
because we have computed our planetary boundary layer for the
ideal case of an absolutely smooth earth. If we are to obtain
results that are more realistic, we must compute the planetary
boundary layer on a rough earth.
In order to take account of surface roughness in the calcu-
lations, we may proceed as follows. First we assume that the
drag of the surface elements will remove momentum from the lowest
layers of the boundary layer in such a way that for a region
between the true surface z = 0 and the small height z = z^ ,
there is no mean velocity; i.e., we assume u = v = 0 for
® — z —-z ' ^-n this region, although the mean velocity is zero,
there certainly may be small turbulent eddies. Therefore, we
will assume that the viscous boundary condition of uju,* = 0
still applies at the true surface z = 0 . We may further
assume that the scale of the eddies is still proportional to the
distance from this true surface, so that in making the calcula-
tions, we specify that Aj = O.?z as we have done in the past.
t
If these new boundary conditions are applied to our equations,
the flow over a rough surface may be calculated.
In Figures 12.6 through 12.10 the results of such a calcula-
tion are presented for a roughness height z of 0.01 meters.
As in the previous example, the results are presented for the
case when u = 10 m/sec and v = 0 . It will be noted from
S. g
-------
12-24
1600
200
800
0>
-------
12-25
-1.0
z, meters
1600
1200
400
zr =0.01 m
0
1.0
2.0
v , (meters/second)
Figure 12.7,
Profile of mean velocity v for steady state
glanetary boundary layer when ug = 10 m/s and
v_ = 0 for a rough earth, zr = 0.01 meters
-------
zr=O.OI m
(v-vg)(m/s)
u-ug, meters/sec
800
Figure 12.8.
Ekman spiral for a rough earth when u = 10 -m/s, v =0, and Zp = 0.01 m
o g
I
IV)
-------
12-2?
800
600
£ 400-
0>
E
*>
N
200-
0
.08 .16
uV (m2/s2)
.24
Figure 12.9. Distributions of u'u', v'v', and w'w1 for the
planetary boundary layer on a rough earth ;_
u = 10 m/sec and v = 0
zr - 0.01 meters for
-------
12-28
z, meters
800
600
u'w'(m2/s2)
Figure 12.10. Distributions of u'w
IT.T I
and u'v1 for the
planetary boundary layer on a rough surface;
zr = 0.01 meters for ug = 10 m/sec and v~ = 0
-------
12-29
Figure 12.8 that the height at which the wind speed first becomes
equal to the geostrophic value / u2 + v2 is now roughly 250 m.
The height of the total Ekman layer is, according to Figure 12.7,
approximately 2000 m. For this roughness condition, the nondimen
siorial friction velocity was computed to be
UK/U = 0.0281
o
The surface angle was found to be 22°. These results are'more in
line with typical experimental results, and we shall presently
compare these results with those found by other investigators.
Figures 12.9 and 12.10 show the distributions of the second-
order velocity correlations that were computed. It will be noted
that the level of the root-mean-square vertical and lateral
velocity fluctuations in the surface layer where / v'v' - / w'w'
= 0..37 are about l\% of the geostrophic wind speed. At about
100 meters above the surface, these fluctuations have fallen off
to about 3% of the geostrophic wind speed. It will also be noted
in Figure 12.10 that no "bump" appears in the energy or shear
correlations for the rough earth where this "bump" had been found
for the smooth earth calculations shown in Figure 12.5.
We may compare the results of these calculations with those
of other investigators. For this purpose, we choose the theoret-
ical results of Blackadar [Ref. 38], Lettau [Ref. 39], and
Appleby and Ohmstede [Ref. 40]. In order to make the comparison,
it is necessary for us to find the relationship between our
roughness height z and the usual roughness height z defined
by other investigators. The equation for the mean velocity in
the immediate neighborhood of the surface when considering the
classical roughness height z is
/52+ v2 = 1 z_ (12.52)
u* 'u* < z
where K is von Karma'n's constant. This may be used to compute
z from the local mean velocity and the surface friction expressed
-------
12-30
through 11^ , viz.,
z = z (12.53)
o r
If we use the results of our calculations for the case just
described where z =0.01 meters, we find from the mean velocity
p-rofile and the friction velocity that the effective z was
0.0071 meters.
We may now compare our results with those of ReTs . 38, 39,
and 1*0 from which the nondimensional friction velocity
u»// u2 + v2 and surface turning angles may be plotted as
g g
functions of the roughness Rossby number
(12.54)
+ v2
2fl sin ZQ
The results of our calculations are compared in Figures
12.11 and 12.12 with those of Refs. 38, 39, and l|0. It is seen
that our results are in good agreement with those of previous
investigators .
Now let us consider the planetary boundary layer when there
are departures from neutral stability, and it is necessary to
include the temperature equations in the solution. With our
present equations, if there exists a steady source or sink of
heat at the boundary of our flow at z = 0 , then there will
never be a steady state solution of the equations in the sense
we have just described for the neutral atmosphere. One can,
however, find time-dependent solutions which are steady state in
the sense that they are cyclic if heat transfer at the boundary
of the flow is cyclic and the net heat added to the atmosphere
is zero. Such a cyclic heating might be representative of the
diurnal effect of the sun on the earth's surface. To obtain
such solutions, one may start at time t = 0 with initial condi-
tions obtained from a neutral solution of the atmospheric boundary
-------
12-31
"
.I0r
.08
.06
.04
.02
8
log,0R0
•-O Rn=oo
10 12
Figure 12.11. Comparison of calculated friction velocities
with results obtained by Blackadar [38],
Lettau [39], and Appleby & Ohmstede [40]
-------
12-32
I.Or
8
Ref. 40
8 10
log,0 R0
—O RosCO
12
Figure 12.12.
Comparison of calculated surface angles in
radians with results obtained by Blackadar [38],
Lettau [39], and Appleby & Ohmstede [40]
-------
12-33
layer and apply at the surface for t > 0 a cyclical heat transfer.
After several cycles of such heat transfer, the characteristics of
the atmospheric boundary layer at each z will become cyclic in
time and independent of the initial conditions that were chosen.
In this way, representative atmospheric boundary layers under a
very wide spectrum of stability conditions can be found. We have
not, as of this writing, obtained solutions of this type to present
here. Studies of this type will comprise part of our effort under
EPA sponsorship during the coming year.
Before leaving this discussion of atmospheric boundary
layers, we should present the equations for the dispersal of a
passive pollutant in such layers that are equivalent to those
which were solved in Section 10. The thin layer equations appro-
priate for use with (12.35) through (12.49) are the following:
,2TT
9C 9C 9C 9 C 9w'C'
a + - a + v - - \> - 2L n? RR^i
,-._.,.. ^- u^ - — -j- v ft — v ~ — ^ V JL c. • \) \) )
9t 9x 9y ° a 2 <^z
C7 2
3u'C' 9u'C' 9u'C'
a , - a , - a
' C ' sin
-------
12-34
9w'C' 3w'C' 9w'C
cos
- w' w
I w I
92w'C'
w'C' + v —
A-, a o n.,,
1 d Z
w'C
— 5
o ,2
A
9w'C' -i
o_
9~z
(12.58)
9T'C' 9T'C' 9T'C'
3C
9t
-
9x
^
9y
=_wtc, _
a 9z
a
3z
3T'C'
9 T'C '
o Z
T'C '
(12.59)
SC
9T
,2
a
9C
+ u
a
9C
,2
9x
+ v
a
9C
3y
= -2w'C'
a
a 9z
3C
'2
2 2
^r11
9 C
a
- 2v
„
oz
C'
a
0 A2
(12.60)
We have not as yet generated any solutions to these equa-
tions . The development of such solutions is a task that will be
carried out under future funding from EPA.
-------
13-1
13. CONCLUSIONS AND RECOMMENDATIONS
This monograph has presented in considerable detail the work
carried out to date by A.R.A.P. towards the development of an
invariant model of the motion in the atmospheric boundary layer
and the dispersal of pollutants in such a layer. While it is
clear that much remains to be done before a final model, complete
with operational programs, is available for production studies of
atmospheric pollutant dispersal, it is clear that a powerful new
tool for the study of such atmospheric problems is emerging. A
great deal of effort over the past two years has been devoted to
the coding and debugging of the programs required for the exer-
cise of the method. Most of the programs required for such
studies are now complete. It is felt that the exercise of these
programs can lead to a detailed understanding of the results of
many experimental studies of pollutant dispersal, many of which
have been difficult to interpret in terms of the simple diffu-
sivity models that have been used in the past. It'is strongly
recommended that EPA continue to support the study of atmospheric
pollution dispersal by the method of invariant modeling.
-------
14-1
REFERENCES
1. Donaldson, Coleman duP.: "A Computer Study of an Analytical
Model of Boundary Layer Transition," AIAA Journal ]_, 2,
(1969), pp. 272-278.
2. Reynolds, Osborne: "On the Dynamical Theory of Incompress-
ible Viscous Fluids and the Determination of the Criterion,"
Phil. Trans. Royal Society, London, A 186, (189*0, p. 123.
3. Kolmogorov, A.N.: "Equations of Turbulent Motion of an
Incompressible Fluid," Izv. Akad. Nauk, SSSR fiz. VI, No. 1-
2, (1942), pp. 56-58.
4. Prandtl, L. and Wieghardt, K.: "liber ein neuses Formelsystem
fur die ausgebeldte Turbulenz," Nachr. Akad. Wiss. Goettingen
!i, 6 (1945).
5. Chou, P.Y.: "Pressure Flow of a Turbulent Fluid between Two
Infinite Parallel Planes," Quart. Appl. Math.3., 3, (1949),
pp. 198-209.
6. Rotta, J.: "Statistische Theorie nichthomogener Turbulenz,"
Z. Physik 129 (195D, p. 547.
7. Glushko, G.S.: "Turbulent Boundary Layer on a Flat Plate in
an Incompressible Fluid," Bull. Acad. Science, USSR, Mech.
Series No. 4 (1965), pp. 13-23-
8. Bradshaw, P., Ferriss, D.H., and Atwell, N.P.: "Calculations
of Boundary Layer Development Using the Turbulent Energy
Equation," J. Fluid Mech. 2_8, 3 (1967), PP. 593-616.
9. Beckwith, I.E. and Bushnell, D.M.: "Detailed Description and
Results of a Method for Computing Mean and Fluctuating
Quantities in Turbulent Boundary Layers," NASA TN D-4815,
1968.
10. Donaldson, Coleman duP. and Rosenbaum, Harold: "Calculation
of Turbulent Shear Flows through Closure of the Reynolds
Equations by Invariant Modeling," NASA SP-216, 1968, pp. 231-
253.
11. Nee, V.W. and Kovasznay, L.S.G.: "Simple Phenomenological
Theory of Turbulent Shear Flows," Phys. Fluids 12 (1969),
pp. 473-484.
12.' Daly, B.J. and Harlow, F.H.: "Transport Equation in Turbu-
lence," Phys. Fluids 13_, 11 (1970), pp. 26-34.
13- Gawain, T.H. and Pritchell, J.W.: "A Unified Heuristic Model
of Fluid Turbulence," J. Comput. Physics 5_, 3 (1970), pp. 383-
405.
14. Saffman, P.G.: "A Model for Inhomogeneous Turbulent Flow,"
Proc. Royal Soc. London A.317 (1970), pp. 417-433.
-------
14-2
15. Rodi, W. and Spalding, D.B.: "A Two Parameter Model of
Turbulence and its Application to Free Jets," Warme- und
Stoffubertragung 3_ (1970), pp. 85-95.
16. Donaldson, Coleman duP., Sullivan, Roger D., and Rosenbaum,
Harold: "A Theoretical Study of the Generation of Atmospher-
ic Clear Air Turbulence," AIAA Journal 10_, 2 (1972), pp. 162-
170.
17.. Lumley, J.L. and Panofsky, H.A.: The Structure of Atmos-
pheric Turbulence^ Interscience Publishers, New York, 1964.
18. Spiegel, E.A. and Veronis, G.: "On the Boussinesq Approxima-
tion for a Compressible Fluid," Astrophysics Journal 131,
(I960), pp. 442-447.
19. Donaldson, Coleman duP.: "Calculation of Turbulent Shear
Flows for Atmospheric and Vortex Motions," AIAA Journal 10,
1 (1972), pp. 4-12. Dryden Research Lecture, 1971.
20. Donaldson, Coleman duP. and Conrad, Peter W.: "Computations
of the Generation of Turbulence in the Atmospheric Boundary
Layer," presented at RAeS/AIAA/CASI Inter'1 Conf. on Atmos-
pheric Turbulence, London, May 1971.
21. Wygnanski, I. and Fiedler, H.: "Some Measurements in the
Self-Preserving Jet," J. Fluid Mech. 38, 3 (1969), pp. 577-
612. ~~
22. Wygnanski, I. and Fiedler, H.: "The Two-Dimensional Mixing
Region," J. Fluid Mech. 4JL, 2 (1970), pp. 327-361.
23. Gibson, M.M.: "Spectra of Turbulence in a Round Jet,"
J. Fluid Mech. 15_, 2 (1963), pp. 161-173-
24. Donaldson, Coleman duP., Snedeker, Richard S., and Margolis,
David P.: "A Study of Free Jet Impingement. Part 2: Free
Jet Turbulent Structure and Impingement Heat Transfer,"
J. Fluid Mech. 4j>., 3 (1971), pp. 477-512.
25. Tollmien, W.: "Berechnung Turbulenter Ausbrietungsvorgange,"
Zeitschrift fur angewandte Mathematik und Mechanik IV (1926),
p. 468.
26. Prandtl, L.: "The Mechanics of Viscous Fluids," in Aerodyna-
mic Theory (W. Durand, editor), Stanford University Press,
1934, pp.. 34-208.
27. Coles, Donald: "Measurements in the Boundary Layer on a
Smooth Flat Plate in Supersonic Flow. Part 1: The Problem
of the Turbulent Boundary Layer," JPL Report 20-69, 1953.
28. Grant, H.L.: "The Large Eddies of Turbulent Motion,"
J. Fluid Mech. 4^ 2 (1958), pp. 149-190.
29. Donaldson, Coleman duP. and Sullivan, Roger D.: "Decay of
an Isolated Vortex" in Aircraft Wake Turbulence and its
Detection (J. Olsen, A. Goldberg, M. Rogers, editors),
Plenum Press, New York, 1971, pp. 389-411.
-------
14-3
30. Donaldson, Coleman duP. and Hilst, Glenn R.: "An Initial
Test of the Applicability of Invariant Modeling Methods to
Atmospheric Boundary Layer Diffusion," A.R.A.P. Report No.
169, 1971.
31. Rotta, J.C.: "Recent Attempts to Develop a. Generally
Applicable Calculation Method for Turbulent Shear Plow
Layers," AGARD Conference on Turbulent Shear Flows, London,
September 1971. AGARD CP-93, pp. A.l-A.ll.
32. Monin, A.S. and Obukhov, A.M.: "Basic Regularity in Turbu-
lent Mixing in the Surface Layer of the Atmosphere," Trudy
Geophysics Inst., ANSSSR 2_4 (1954), p. 163.
33. Rose, W.G.: "Results of an Attempt to Generate Homogeneous
Turbulent Shear Plows," J. Fluid Mech. 2J5, 1 (1966), pp. 97-
120.
34. Rose, W.G.: "Interaction of Grid Turbulence with a Uniform
Mean Shear," J. Fluid Mech. _4_4, 4 (1970), pp. 767-779.
35. Fuquay, J., Simpson, C.L., and Hinds, W.T.: "Prediction of
Environmental Exposures from Sources near the Ground Based
on Hanford Experimental Data," J. Appl. Meteorol. 3_, 6
(1964), pp. 761-770.
36. Lewellen, W. Stephen and Teske, Milton: "Prediction of the
Monin-Obukhov Similarity Functions from an Invariant Model
of Turbulence," A.R.A.P. Report No. 185, November 1972.
37. Hess, Seymour L.: Introduction to Theoretical Meteorology,
Holt, Rinehart, and Winston, 1959, pp. 161-173.
38. Blackadar, Alfred K.: "The Vertical Distribution of Wind
and Turbulent Exchange in a Neutral Atmosphere," J. Geophys.
Res. 6_7, 8 (1962), pp. 3095-3102.
39. Lettau, H.: "Theoretical Wind Spirals in the Boundary Layer
of a Barotropic Atmosphere," Beitr. Phys. Atmosph. 35• 3-4
(1962), pp. 195-212.
40. Appleby, J.F. and Ohmstede, W.D.: "Numerical Solution of
the Distribution of Wind and Turbulence in the Planetary
Boundary Layer," Proc. Army Sci. Conf. 1_ (1964), pp. 85-99.
41. von Rosenberg, D.U.: Methods for the Numerical Solution of
Partial Differential Equations, American Elsevier Publ. Co.,
Inc., New York, 1969.
42. Carnahan, B., Luther, H.A., and Wilkes, J.O.: Applied Numer-
ical Methods, John Wiley & Sons, New York, 1969, pp. 449-461,
43. Crow, Steven C.: "Viscoelastic Properties of Fine-Grained
Incompressible Turbulence," J. Fluid Mech. 33, 1 (1968),
pp. 1-20.
-------
A-l
APPENDIX A. NUMERICAL INTEGRATION TECHNIQUE
The equations discussed in this report are parabolic partial
differential equations. They are parabolic in the sense that the
time derivative 9/9t is first-order in time and is present in
2 2
every equation with at most a second-order derivative 3 /9z in
a spatial dimension. In this appendix, we will explain the
approximate numerical technique used by A.R.A.P. to obtain reason-
able solutions to these differential equations.
Various numerical techniques are, of course, available to
solve the types of equations we have discussed in this report.
For our purposes, the implicit technique [Refs. 4l and 42] seems
to be best, for a number of reasons. Its main attraction is its
strong stability. The computerized equations should converge to
a solution, even though it may not be the correct one. The
simplicity of our boundary conditions, with all values going to
zero at the earth's surface or to given known values at the upper
surface, makes the application of the implicit technique almost
straightforward. Another reason is computer environment.
A.R.A.P.'s present computer is a 16K core (16 bits/word) Digital
Scientific Corp. META-4 with secondary storage on a CalComp DS-12
disk drive. The relatively slow speed of computation (as demon-
strated below) requires us to use a technique that permits
solution as quickly as possible to the time-dependent problem.
The implicit technique gives us that solution in a time that is
faster (for the number of equations we typically handle) than a
standard iterative technique.
To apply any of these several numerical techniques to a
continuous differential equation, we first replace the differen-
tial equation with a finite difference equation. We represent
the continuum spatial dimension by a series of discrete points.
Three such points in the z direction are represented in Figure
A.. 1 by the points at j-1, j, and j+1. In the same manner, we
replace continuous time with discrete time, as with the time
points n and n+1 , as shown.
-------
A-2
n
» t
n + l
Figure A.I.
The discretized time-space plane for the two-
dimensional implicit solution technique
-------
A-3
We solve the differential equations by starting with a
given set of initial conditions at all z values, and then
follow the solution for increasing time. We form an approxima-
tion to the continuous derivatives present in our equations.
These difference approximations are then substituted into the
differential equations to obtain a difference equation involving
the values of the dependent variables as--a function of the j
and n positions. To solve this equation, we assume that the
solution at tn is known (computed from the last step or a
given initial condition) and use the implicit technique given
below to obtain the solution at time tn . As we stated before,
this technique is unconditionally stable, but it may not be
consistent in that the difference equation may not relax to the
differential equation in the limit of Az, At -»• 0 . To ensure
some success in solving our equation (since existence and unique-
ness of our solutions have never-been proven), we attempt
internally within the computer program to maintain an "optimum"
number of points in the profile in z at any one time step (a
test based on curvature criterion) and an "optimum" step size-in
t based on maximum permitted change per step. This optimization
enables us to compromise between accuracy and time considerations;
as the time step decreases, the accuracy of the results increases,
up to a point, but the total computation time increases accord-
ingly. We know that our solutions are 1% accurate for the various
laminar cases we have checked, and we strive to keep our turbulent
solutions to an accuracy of 5%.
We must then approximate the four possible differentials
3/3t , 3/3z , 32/3z2 , and 3 (f (z) 3/3z)/3z .. They become
n+1 _ n
3t , n+1 ,n
t - t
3z h. + h
T
-------
A-4
2 (un«. _un+l) _ 2 /n+1 _ n+l<
»2,, n+ V j + 1 j / h V j J-l.
LU , _+ z (A>3)
h. + h
f + f f + f
' "in ,,n+i\ J :
8 (f(z) £)* 1
dZ \ dZ /
h+ + h_
(A.4)
where
and u represents the dependent variable. With Eqs. (A.I)
through (A.4) substituted into the differential equation of
interest, we have replaced the equation by a forward-time-centered-
space differencing that is first-order accurate in time and second-
order accurate in space.
To get a feel for this technique, we now use Eq. (10.5) for
C''~ as a very simple application. Since v'C' , w'C' , and C
r r' sr r
are presumably computed elsewhere, Eq. (10.5) contains only one
2
dependent unknown, C' , and may be rewritten as
9C'2 \ » / 30
3t = ~ H + 9?
2 \ 2
i^ \ r. i*
(A.5)
0 X2
-------
A-5
where
9C
H =
_
2v'C'
P
2w'Cf
P
9z
(A.6)
and is known (in general, several equations will be coupled
together to form a vector of unknowns, as shown below). For
the one-space-dimensional problem, we may look at only the z
and t dependence and assume that 3/3y = 0 in Eq. (A.5). The
substitution of Eqs. (A.I) through (A.4) into (A.5) yields the
form
where
vnun+1
Jj-l
Ynun+1
J J
7nun+1
(A.7)
V —
hjh
(A.8)
nil
j .n+1 ,n h . + h
U ~ U i ~"
Aq . , , + Aq . Aq . + Aq . ,
J+l J , J J-l
I ' h+
/I 1 \1 2v
'"^o \ h h / 2
+ - ' X
h
(A.9)
zn =
j
Dn =
n
u
tn+1 - tn
- H
and u
n+1
n+1
u
J
n
Aq.+1
,
and
n+1
Aq
(A.10)
(A.11)
are the unknown values of
at j-1 , j , and j+1 . A similar linear equation may be
written for every other interior point in our discrete space.
Thus, if we have m points in the profile, we will have m
m+2 unknowns. Since both
difference equations (A. 7) with
boundaries are known, U
Q
and
are known values. We,
-------
A-6
therefore, have a sufficient number of equations to solve for
all of the u} at the new time step.
J
We choose to use the general tridiagonal algorithm (called
the Thomas algorithm in Ref. 4l) to systematically solve for
these unknowns. This method requires a double sweep of the z
profile and is actually a straightforward Gaussian elimination.
We begin at j = 0 knowing un and work towards the point
j = m+1 by successively solving the following set of algebraic
equations.
n+l _ n n
-1
J ^ J '
Since u^ is known, the first term in Eq. (A.7) is known for
n
j = 1 and may be combined with D. , enabling us to write
vn - n J
A, - (J .
The sweep is continued until j = m where we downsweep
back to j = 1 using the formula
= Wn+1 - Gn+1u
j j :
Since u"! t is known and all the W1?"1"1 and G1?"1"1 are also
m+1 j j
known, the downsweep give the successive values of u at all
the discrete z points j for the time tn . The procedure
is then repeated for the next time t and so on. For other
boundary conditions, the technique needs some modification at
the end points.
If u. is a vector of-unknowns U. (as, for example, C ,
T^Cf7" , w'C' , and T'C' in Section 10), then X. , Y. , and
-------
A-7
Z. are square matrices and D. is a vector. The inversion in
J J
Eqs. (A.13) and (A.14) is then a matrix inversion. The implicit
technique is seen to be defeating if the size of the matrix (the
number of unknowns) is abnormally high. In our facility, a
matrix size of about 15 would appear to be prohibitive.
A typical running time for our various computer programs
is one hour. If the program is small enough to fit into the 16K
core, its running time would be reduced to 1/2 hour. Comparable
running times on an 1108 Univac and CDC 6600 would be 5 and 2
minutes, respectively. The number of machine operations per
second on the META-4 is approximately 200,000.
Equation (A.5) may also be used to demonstrate the Alter-
nating Direction Implicit (ADI) method used to solve the spatial
pollutant dispersal equations of Section 10. We now have the
additional spatial dimension y with discrete points ..., k-1,
k., k+1, ... as shown in Figure A.2. Our aim is to advance the
plane of solution at tn and'all y and z to the time tn
A straightforward substitution of the finite difference approxi-
mations (substituting y and k where appropriate) into (A.5)
gives the form
In un+1 + In un+1 + n n+1
•Lk-l,juk-l,J xk,J k,j
, Tn n+1 , Tn ,, n ,. ,,-.
+ I. . , u, . , + I, . ,-,u. . ,, = D, . (A.16)
where the I's involve the spacing factors and other known values,
Each equation now contains five unknowns. With the known bound-
ary conditions, there are as many equations as unknowns, but the
difficulty of the problem has increased considerably. An elimina-
tion technique could be employed at this point, but it is far
easier from our point of view to use the ADI method. In this
method, every time step is split and two half-steps are performed,
In the first half-step, one of the spatial dimensions is swept up
-------
A-8
k-l,J
i
k, j + l
^
1
k.j-l
Figure A.2. The discrete spatial plane at time step tn
-------
A-9
and down; in the second half-step, the other dimension is swept.
During the first sweep (say, z) from tn to tn , the y
derivatives are held fixed at their t values. During the second
sweep in y from t to t , the z derivatives are held
fixed at their tn+1/2 values. It can be shown [Ref. 41] that
while this double sweep (alternating direction) uses two sets of
equations similar to Eq. (A.7)5 its full-step results are equivalent
to using Eq. (A.16). Of course, the complexity of the computer
program and its running time are increased significantly.
-------
BIBLIOGRAPHIC DATA
SHEET
1. Report No.
EPA-R4-73-Ol6a
3. Recipient's Accession No.
4. Title and Subtitle
ATMOSPHERIC TURBULENCE AND THE DISPERSAL OF
ATMOSPHERIC POLLUTANTS
5. Report Date Approved
March 1§73
6.
7. Author(s)
Coleman duP.
Donaldson
8- Performing Organization Rept.
N°- 186, Vol. I
9. Performing Organization Name and Address
Aeronautical Resea.rch
50 Washington Road
Princeton, New Jersey
Associates of Princeton,Inc
085^0
10. Project/Task/Work Unit No.
Element A-11009
11. Contract/Grant No.
EPA 68-02-0014
12. Sponsoring Organization Name and Address
U.S. Environmental Protection Agency
Office of Research and Monitoring
Washington, D.C. 20460
13. Type of Report & Period
Covered
Interim
14.
15. Supplementary Notes
16. Abstracts
A detailed derivation of the equations describing the generation of
turbulence and the transport of pollutants in the atmosphere is given.
The equations are closed at the order of the second-order turbulent
correlations through the use of models of the higher-order correla-
tions. The resulting invariant second-order closure model is applied
to three problems: 1) the generation of turbulence in the surfa.ce
layer (computed results are compared with AFCRL experimental data);
2) the transport of passive pollutants in the atmosphere (computed
results are compared with experimental data and suggestions for
improving the model are made); 5) the complete planetary boundary layer
for the case of neutral stability for both smooth and rough surfaces.
The results of computations are compared with the results of previous
investigations.
17. Key Words and Document Analysis. 17o. Descriptors
Turbulence
Turbulence Models
Atmospheric Motion
Atmospheric Turbulence
Atmospheric Transport
Pollutant Dispersal
Atmospheric Boundary Layer
Atmospheric Surface Layer
Ekman Layer
17b. Identifiers/Open-Ended Terms
Second-order Closure
Second-order Modeling
Invariant Modeling
Boundary Layers
Surface Roughness
17e. COSAT1 Fie Id/Group
18. Availability Statement
Unlimited
19. Security Class (This
Report)
UNCLASSIFIED
20. Security Class (This
UNCLASSIFIED
Pace
21. No. of Pages
210
22. Price
FORM NTIS-35 IREV. 3-72)
USCOMM-DC M8S2-P72
-------
INSTRUCTIONS FOR COMPLETII graphic Data Sheet based on COSATI
Guidelines to Format Standards for Scientific ana lecnnicai neports Prepared by or for the Federal Government,
PB-180 600).
1. Report Dumber. Each individually bound report shall carry a unique alphanumeric designation selected by the performing
organization or provided by the sponsoring organization. Use uppercase letters and Arabic numerals only. Examples
FASEB-NS-87 and FAA-RD-68-09.
2. Leave blank.
3. Recipient's Accession Number. . Reserved for use by each report recipient.
4- Title and Subtitle. Title should indicate clearly and briefly the subject coverage of the report, and be displayed promi-
nently. Set subtitle, if used, in smaller type or otherwise subordinate it to main title. When a report is prepared in more
than one volume, repeat the primary title, add volume number and include subtitle tor the specific volume.
5- Report Dote, l-.uch report shall carry a date indicating at least month and year. Indicate the basis on which it was selected
(e.g., date of issue, date of approval, date of preparation.
6. Performing Organization Code. Leave blank.
7. Authors). Give name(s) in conventional order (e.g., John R. Doc, or J.Robert Doe). List author's affiliation if it differs
trom the performing organization.
8. Performing Orgonizotion Report Number. Insert if performing organi/.at ion wishes to assign this number.
9. Performing Organization Name and Address, dive name, street, city, state, and zip code. List no more than two levels of
an organizational hierarchy. Display the name of the organization exactly as i: should appear in Government indexes such
as USGRDR-I.
10. Project/Tosl'/Work Unit Number. Use the project, task and work unit numbers under which the report was prepared.
11. Contract/Grant Number. Insert contract or grant number under which report was prepared.
12. Sponsoring Agency Nome and Address. Include zip code.
13. Type of Report and Period Covered. Indicate interim, final, etc., and, if applicable, dates covered.
14. Sponsoring Agency Code. Leave blank.
15. Supplementary Notes. Knter information not included elsewhere but useful, such as: Prepared in cooperation with . . .
Translation of ... Presented at conference of ... To be published in ... Supersedes . . . Supplements . . .
16. Abstract. Include a brief (200 words or less) factual summary of the most significant information contained in the report.
l! the report contains a significant bibliography or literature survey, mention it here.
17. Key Words and Document Analysis, (a). Descriptors. Select from the Thesaurus of Kngineering and Scientific Terms the
proper authorized terms that identify the major concept of the research and are sufficiently specific and precise to be used
as index entries for cataloging.
(b). identifiers and Open-Ended Terms. Use identifiers for project names, code names, equipment designators, etc. Use
open-ended terms written in descriptor form for those subjects for which no descriptor exists.
(c). COSATI Field/Group. Field and Group assignments are to be taken from the 1965 COSATI Subject Category List.
Since the majority of documents are muhidisciplinary in nature, the primary Field/Group assignment(s) will be the specific
discipline, area of human endeavor, or type of physical object. The application(s) will be cross-referenced with secondary
Field/Group assignments that will follow the primary posting(s).
18. Distribution Statement. Denote releasability to the public or limitation for reasons other than security for example "Re-
lease unlimited". Cite any availability to the public, with address and price.
19 & 20. Security Classification. Do not submit classified reports to the National Technical
21. Number of Poges. Insert the total number of pages, including this one and unnumbered pages, but excluding distribution
list, if any.
22. Price. Insert the price set by the National Technical Information Service or the Government Printing Office, if known.
FORM NTIS-33 (REV. 3-721 USCOMM-DC 149S2-P72
------- |