WATER POLLUTION CONTROL RESEARCH SERIES • 16130DW010/70
I
MATHEMATICAL MODELS
FOR THE PREDICTION OF
TEMPERATURE DISTRIBUTIONS
RESULTING FROM THE DISCHARGE
OF HEATED WATER
INTO LARGE BODIES OF WATER
ENVIRONMENTAL PROTECTION AGENCY • WATER QUALITY OFFICE
-------
WATER POLLUTION CONTROL RESEARCH SERIES
The Water Pollution Control Research Series describes the
results and progress in the control and abatement of pollu-
tion of our Natio&s waters. They provide a central source
of information on the research, development, and demon-
stration activities of the Water Quality Office, Environ-
i ntal Protection Agency, through inhouse research and grants
and contracts with Federal, State, and local agencies, re-
search institutions, and industrial organizations.
Inquiries pertaining to the Water Pollution Control Research
Reports should be directed to the Head, Project Reports
System, Office of Research and Development, Water Quality
Office, Environmental Protection Agency, Washington, D C 20242
-------
MATHEMATICAL MODELS FOR THE PREDICTION OF
TEMPERATURE DISTRIBUTIONS RESULTING FROM
THE DISCHARGE OF HEATED WATER INTO LARGE
BODIES OF WATER
byv
Robert C. Y. Koh
Loh-Nien Fan
Tetra Tech, Inc.
630 N. Rosemead Blvd.
Pasadena, California 9110?
for the
WATER QUALITY OFFICE
ENVIRONMENTAL PROTECTION AGENCY
Program #16130 DWO
Contract #14-12-570
October, 1970
For sale by the Superintendent ol Documents, U.S. Government Printing Office
Washington, D.C., 20402 - Price $1.75
-------
EPA Review Notice
This report has been reviewed by the Water Quality Office,
EPA, and approved for publication. Approval does not signi-
fy that the contents necessarily reflect the views and poii—
cies of the Environmental Protection Agency, nor does mention
of trade names or commercial products constitute endorsement
or recommendation for use.
11
-------
TABLE OF CONTENTS
1. INTRODUCTION
2. CONCLUSIONS AND RECOMMENDATIONS
3. INITIAL MIXING AND SURFACE SPREADING
DUE TO SUBSURFACE DISCHARGE
3. 1 ______________
3. 2
Introduction
Initial Mixing Phase: Multiple Submerged
Buoyant Jets Discharged into an Arbitrarily
Stratified Ambient
3. 2. 1 Formulation
3. 2. 2 Method of Solution and Examples
3. 3 Time Dependent Surface Spreading of a
Buoyant Fluid
3. 3. 1 Two-Dimensional Case
3. 3. 2 Axisymmetric Case
3. 3. 3 Comparison withExperiments
3. 3. 4 Numerical Solutions
Summary and Discussion
FACE HORIZONTAL BUOYANT JETS
Introduction
Formulation and Solutions for the Two-
dimensional Problem
4. 2. 1 Derivation of Equations
4. 2. 2 General Discussion of the Equations
and Solutions
4. 2. 3 Solution of the Equations
4. 2. 4 Matching of Solutions
4. 2. 5 Summary and Discussions
4. 3 Axisymmetric Surface Buoyant Set
4. 4 Example Applications
3. 4
4. SUB
4. 1
4. 2
Page
I
7
15
15
16
17
28
35
35
40
44
47
49
51
51
56
62
67
73
87
92
97
111
ill
-------
Page
5. PASSIVE TURBULENT DIFFUSION FROM A
CONTINUOUS SOURCE IN A STEADY SHEAR
CURRENT OR AN UNSTEADY UNIFORM CURRENT
WITH UNSTEADY SURFACE EXCHANGE
5.1 Introduction 115
5.2 Derivation of Basic Equations 119
5.2.1 The Problem for Excess Heat 119
5. 2. 2 The Problem for a Tracer 122
5. 2. 3 The General Problem 123
5. 3 Steady Release in a Steady Environment 126
5. 3. 1 Formulation 126
5. 3. 1. 1 Source Conditions 126
5. 3. 1. 2 Environmental Character- 127
istics
5. 3. 1. 2a Horizontal Dif- 128
fusion Coefficient K
z
5. 3. 1. 2b Vertical Diffusion 131
Coefficient Ky
5. 3. 2 Method of Moments 137
5. 3. 3 Limiting Solution for Cases with Zero 139
Vertical Transport
5. 3. 4 Dimensionless Equations and Numerical 141
Solutions
5.4 Continuous Release of Heat into a Uniform But 156
Time-Varying Environment
5.4.1 Formulation 156
5.4.2 Method of Moments 157
5.4. 3 Dimensionless Equations and Numerical 161
Solutions
5. 5 Summary and Discussions 164
6. APPLICATION OF THE MATHEMATICAL MODELS TO 167
PRACTICAL PROBLEMS
6.1 Subsurface Discharges 168
6. 2 Surface Discharge 169
REFERENCES 171
APPENDICES A, B, C, D Computer Programs 177
iv
-------
LIST OF FIGURES
Figure No. Title Page No .
3. 1 Definition sketch. 18
3. 2 Jet interference. 18
3. 3 Zone of flow establishment in a submerged jet. 27
3. 4 Predicted trajectories of multiple buoyant jets 31
in a uniform environment.
3. 5 a, b Predicted jet centerline excess temperature of 32-33
multiple buoyant jets in uniform environment.
3. 6 Predicted trajectories of multiple buoyant 34
jets in stratified environment.
3. 7 Definition sketch. 36
3. 8 Definition sketch. 36
3. 9 Definition sketch. 41
3. 10 Definition sketch. 41
3. 11 Comparison of theory with experiments for the 46
surface spreading of buoyant fluid.
3. 12 Growth of a spreading surface pool of buoyant 48
fluid.
4. 1 Definition sketch. 57
4. 2 Entrainment coefficient e as function of Richardson 57
number Ri.
4. 3 Definition sketch. 59
4.4 The zones in a surface horizontal buoyant jet. 61
4. 5a Predicted jet thickness and density deficiency in 77-78
a surface horizontal buoyant jet for k = I/B = 0
(two dimensional case).
V
-------
Figure No. Title Page No .
4. 5b Predicted jet thickness in a surface horizontal 79
buoyant jet for the case when k> kcr+ (two dimen-
sional case).
4. 6 Definition sketch. 80
4.7 Critical relation F = F (s). 88
o ocr
4. 8 Division of parameter space into regions of 91
different flow pattern.
4. 9 Definition sketch. 98
4.10 Predicted jet thickness in a surface horizontal 109
buoyant jet (axisymmetric case).
4.11 Predicted jet density deficiency in a surface 110
horizontal buoyant jet (axisymmetric case).
5.1 Various stages of mixing of submerged cooling 116
water discharges.
(a) Effluent field established at water surface. 116
(b) Effluent field trapped below water surface. 116
5.2. Passive turbulent diffusion cases studied. 117
(a) A steady continuous source in a steady shear 117
current with constant surface heat exchange.
(b) An unsteady continuous source in a uniform 117
unsteady current with time varying surface
heat exchange.
5. 3 Flow configurations in cases with steady releases. 124
(a) Submerged release. 12.4
(b) Surface release. 125
5. 4 Horizontal diffusion coefficient as a function of 129
horizontal scale (from Orlob (1959)).
5. 5 Correlation of Ky with density gradient. 134
5. 6 Dependence of Ky 1 on sea state. 136
vi
-------
Figure No. Title Page No .
5. 7 Profiles of Ky and u used in study. 146
5. 8a Vertical distribution of c (x, y) for PTD cases 149
CN 200, CC 200, CC 100, CC 220 (Ky—profile
uniform).
5. 8b Vertical distribution of c (x, y) for PTD cases 150-151
SN 200, SN 210, SC 100, C 200 (Ky-profile
not uniform, surface release).
5. 8c Vertical distribution of c (x, y) for PTD cases 152
TN 100, TN 120, TN 200, rc 100, TC 200, TC 202
(Ky-profile not uniform, subsurface release).
5. 9a Width of diffusing plume for PTD cases CC 200, 153
CN 200, CC 100, CC 220 (Ky-profile uniform).
5. 9b Width of diffusing plume for PTD cases SN 200, 154
SN 210, SC 100, SC 200 (Ky-profile not uniform,
surface release).
5. 9c Width of diffusing plume for PTD cases TN 100, 155
TN 120, TN 200, TC 100, TC 200, TC 202
(Ky-profile not uniform, subsurface release).
5. 10 Representative sketches in cases of unsteady 159
turbulent diffusion.
5. Il Comparison of PTJJ with UTD to illustrate effect 163
of unsteady current, rate of heat discharged, and
surface heat exchange coefficient.
vii
-------
LIST OF SYMBOLS
For simplicity, symbols of secondary importance which appear only
briefly in the text are omitted from the following list.
Chapter 3 .
CD drag coefficient
3 CD
D 8
D jet diameter
0
E entrainment function
F density deficiency flux
F 1 initial value for density deficiency flux
G temperature deficiency flux
initial value for temperature deficiency flux
K coefficient
L jet spacing
M kinematic momentum flux
momentum flux
M 1 initial value of kinematic momentum flux
Q volume flux
viii
-------
Q 1 initial value of volume flux
discharge
T centerline temperature
Ta ambient temperature
T temperature
T 1 initial value of centerline temperature
a thickness of spreading layer
a 0 initial thickness of spreading layer
b jet characteristic width or radius of spreading layer
b initial radius of spreading layer
d. depth of jet
f buoyancy force
g gravitational acceleration
p pressure
q 0 discharge
s coordinate along jet path
t time
ix
-------
u velocity along jet
u centerline velocity along jet
x horizontal coordinate
y vertical coordinate
a entrainment coefficient
entrainment coefficient for slot jet
a entrainment coefficient for round jet
r
dimensionless radius of spreading layer
eddy viscosity
coordinate normal to s
e angle of jet trajectory
00 initial angle of jet discharge
p centerline density
00 ambient reference density
ambient density
o * density
jet discharge density
x
-------
Chapter 4 .
F
F
0
F 2
F
ocr
K
K
cr+
K
cr-
spreading ratio for slot jet
spreading ratio for round jet
dimensionless time
Xr
1
0
2 T
densimetricFroude number = U /g —
source densimetric Froude number = U / g — h
0 0
0
value of F after an internal hydraulic jump
critical value of densimetric Froude number
surface heat exchange coefficient
upper critical value of K above which jet type
solution exists
lower critical value of K below which the source
is inundated
Reynolds number
Richardson number = -
critical Richardson number
density deficiency
source density deficiency
R
Ri
Ri
cr
T
T
0
xi
-------
U velocity
U 0 source velocity
e entrainment coefficient
e entrainment coefficient for F =
0
g gravitational acceleration
h thickness of surface jet
h initial thickness of surface jet
0
k dimensionless surface heat exchange coefficient
upper critical value of k
k lower critical value of k
cr-
p pressure
q 0 , q discharge
r radial coordinate
1
u horizontal velocity or dimensionless velocity
v vertical velocity
x horizontal coordinate
y vertical coordinate
z vertical coordinate in axisymmetric case
xii
-------
ct, a 1 , a 2 coefficients
dimensionless coefficient, =
free surface elevation
e density deficiency
T shear stress (kinematic)
T shear stress at free surface (kinematic)
S
shear stress at interface (kinematic)
e shear coefficient
o density
ambient density
0
1 jet density
Chapter 5 . (Note: Primed quantities are dimensionless)
A dissipation parameter (K Ac
4/3
AL dissipation parameter (K AL L
Ch specific heat of water
E equilibrium temperature
Fco source strength
H heat content of ambient above a given reference level
a
xiii
-------
H total heat content above a given reference level
K longitudinal diffusion coefficient (horizontal)
K vertical diffusion coefficient
y
K vertical diffusion coefficient at sea surface
yl
K lateral diffusion coefficient (horizontal)
z
surface heat exchange coefficient
K kinematic surface heat exchange coefficient (KE / o
decay coefficient
L source width
0
L plume width
R. Richardson number ( )
dy
T temperature excess
T surface ambient water temperature
a
T surface water temperature
p
c concentration
c zeroth moment of c
0
c 1 first moment of c
c 2 second moment of c
c maximum value of c in z
max
xiv
-------
g gravitational acceleration
h source thickness
0
hb depth of water
t time
u velocity in x- direction
u characteristic velocity
v velocity in y - direction
w horizontal coordinate, in direction of current
x terminal x of interest
K1’ K2’
y - coordinates in defining K - profiles (see Fig. 5. 7)
K3’ K4’
e1’ e y - coordinates in defining u - profiles (see Fig. 5.7)
source level
z horizontal coordinate, transverse to current
coefficients in defining K profile (see Fig. 5. 7)
density gradient
p density
plume width characteristic
xv
-------
source width characteristic
X dimensionless dissipation parameter = Ax
I 2/3
‘(LGZ (o,y)u
xvi
-------
CHAPTER 1 INTRODUCTiON
The rise in the production of electric power has resulted in the attendant
generation of large quantities of waste heat. This waste heat is usually
disposed of either to the atmosphere through cooling towers or ponds or
to adjacent bodies of water. In order to properly manage the vast quantities
of waste heat which will be produced in the future, it is necessary to de-
velop a body of knowledge on the transport behavior and the effects of heat
on the total environment. One important item in this necessary body of
knowledge is the ability to predict the temperature distribution in the en-
vironment given the method of waste heat discharge and the characteristics
of the environment. The present investigation is concerned with the develop-
ment of prediction methods in the case when the waste heat is discharged
into a large body of water. In this case, two limiting schemes can be
envisioned for the method of discharge of heated water. First, one may
employ a multiport diffuser submerged at some depth to promote much
initial dilution such as is done for sewage. Alternatively, the other
extreme would be to “float’ the warm water on the surface, resulting in a
minimum of initial dilution while maximizing the rate of heat toss to the
atmosphere.
In order to properly evaluate the effects of various discharge schemes on
the environment, it is necessary to be able to predict the resulting tem-
perature distributi n given the discharge scheme. This would also provide
a rational basis for the design of the discharge structure.
Of importance in this overall problem of excess temperature prediction
are the following phenomena and their interrelationships:
a) momentum of the discharge . For discharge schemes
employing relatively large efflux velocities, the behavior
of the effluent near the source is strongly influenced by this
momentum and the mixing phenomenon may resemble that
in a jet.
b) buoyancy of the discharge . Since the effluent is warmer
(and hence lighter) than that of the receiving waters, there
I
-------
is a tendency for it to float on top of the ambient cooler
(and hence heavier) water.
c) dispersion due to ambient turbulence . Even in the absense
of momentum and buoyancy, the introduction of any mis -
cible tracer into a body of water would result in the dis-
persion of the tracer due to existing turbulence in the ambient.
d) ambient density stratification . The water in a typical lake,
reservoir or the ocean is often density stratified particularly
in the summer months. The stable stratification has the
profound effect of suppressing vertical turbulence and dis-
persion. In addition, the warm effluent, if discharged at
the surface, tends to float on top of the cooler ambient and
enhance the existing stratification.
e) ambient current structure . The effluent, other than undergo
motions induced by its own momentum and buoyancy, would
also be advected by any ambient currents which may be
there. These currents may change with time and location.
f) solid boundaries . The presence of boundaries (both the shore
and the bottom) also affects the dispersion of the effluent
and the resulting temperature distribution in the local
environment.
g) surface heat exchange . The warm effluent exposed to the
atmosphere would gradually lose heat to the atmosphere,
altering the temperature and density of the water, particularly
when there is wind.
When heated cooling water is released from a power plant through a dis -
charge structure, all the mechanisms discussed above and their inter-
relationships play a role in influencing the resulting temperature distri-
bution in the receiving water. However, different mechanisms would
dominate in different regions of the induced flow field. For example, near
the source of discharge, it can be expected that the momentum and buoyancy
of the effluent would be important in influencing the mixing process. On
the other hand, far from the source, it may be imagined that the ambient
2
-------
currents and turbulence would be dominating factors and the dispersion may
be thought of as more or less passive. In between, all the mechanisms
may contribute to the dispersion process.
The phenomenon of dispersion and mixing of one fluid with another has been
studied by numerous investigators. There is, for example, a body of
knowledge on the dispersion of sewage effluent discharged from outfalls.
These can be applied, with some modifications, to the problem of the
transport behavior of heated water discharged below the surface. Also
there are some laboratory studies on the dispersion of heated water dis -
charged at the surface into a laboratory tank containing cooler water.
These investigations will now be briefly summarized.
a) Dispersion of Sewage Effluent . Present methods of ocean
sewage disposal typically discharges the effluent through a
multiport outfall diffuser submerged at about 200 ft. depth.
The mixing processes undergone by the effluent can be
divided into three separate phases. First, the effluent
undergoes jet diffusion through entrainment of the ambient
water as it rises in the form of a buoyant jet or plume.
Second, on reaching its terminal level of ascent which may
be the surface it spreads out horizontally due to the density
difference or difference in density stratification. Third,
it further diffuses in the prevailing ocean current. The
phenomenon of buoyant jets and plumes has been studied by
many investigators including Abraham (1963), Brooks and
Koh (1965) and Fan (1967). Almost all those investigations
are for the case when the ambient fluid is motionless.
The further passive diffusion of the diluted effluent in the
prevailing current has been studied by Brooks (1960) for
the case of vertical uniformity and constant current velocity
with the lateral dispersion characterized by a power de-
pendence on the plume width. Edinger and Polk (1969)
analyzed the similar problem including vertical variations
3
-------
but all dispersion coefficients and currents were assumed
constant.
The second phase of the horizontal spreading phenomenon
has, to date, not been studied to nearly the same extent.
Some small scale laboratory experiments have been performed
by Sharp (1969), for the case when the spreading is on the
surface of a quiescent laboratory tank.
b) Studies on Thermal Dispersion on the Surface . The growing
concern over thermal pollution has lead to several laboratory
and theoretical investigations of the dispersion of heated
water discharged on the surface of a body of cooler ambient
water. Jen, Wiegel and Mobarek (1966), Hayashi and Shuto
(1967), and Stefan and Schiebe (1968) performed laboratory
experiments where the warm water was discharged from a
finite source horizontally into a quiescent cooler ambient.
Measurements were made for a variety of cases. Wada
(1966) and Hayashi and Shuto further advanced a theory for
this problem. However, it is applicable only for extremely
low discharges and small temperature differences. The
two-dimensional case of the same problem was also in-
vestigated experimentally by Stefan and Schiebe with several
interesting results.
From the above brief summary of previous work on the problem of prediction,
it is clear that no general method exists by means of which the temperature
distribution resulting from waste heat discharge can be predicted. It is
the purpose of the present investigation to advance the status of knowledge
on this problem. It should be pointed out that no general prediction method
is developed in this report. Bather, several simpler prediction models are
developed each applicable under differing circumstances. It is believed that
these models bring us closer to the time when a general model may be
formulated.
4
-------
It will become obvious on reading the subsequent chapters of this report
that the general mixing and dispersion phenomenon which ensues following
the discharge of heated cooling water into a large body of water is highly
complex and thus difficult to analyze. Not only is the hydrodynamical
aspects complicated, the prevailing ambient conditions are usually not
deterministic and can only be statistically described. Moreover, the
interplay of the many mechanisms such as source momentum, buoyancy,
surface heat loss, ambient currents and so on, make the problem difficult
to describe even qualitatively. Therefore, before any attempt is made to
develop a general prediction model, it is necessary to examine the
significance and interrelationships between these mechanisms taken several
at a time. This then is the philosophy of the present investLgation. It is
found that the interrelationship between these mechanisms are such that
the flow field is sometimes entirely different from what may be intuitively
expected.
This report has been divided into several chapters. Chapter 3 deals with
the initial and intermediate phases of mixing in the event the discharge is
made at depth. The problem of a row of equally spaced round buoyant jets
discharging at an arbitrary angle into an arbitrarily density-stratified body
of water is solved in Sec. 3. 2. The unsteady surface spreading of a warm
fluid on top of a cooler ambient is analyzed in Sec. 3. 3.
Chapter 4 deals with the case when the discharge is made at the surface.
The two-dimen ional case is treated in detail while the axisymmetric case
is also examined. The effects of source momentum, source buoyancy,
interfacial shear, surface heat exchange and entrainment are all included.
It is found that the flow field can be entirely different depending on the
relative importance of these mechanisms. In some cases, the flow field
is like a jet while in others, it is like a two-layered stratified flow. Under
certain conditions, the flow field consists of a jet type region near the source
followed by a stratified flow region with an internal hydraulic jump joining
the two. Given the source characteristics and the ambient conditions, the
model developed can predict the flow field and the temperature distribution
and also locate the hydraulic jump, if it occurs.
5
-------
In Chapter 5, two mathematical models have been developed for the
case of passive turbulent diffusion from a continuous source in a uni-
directional current. The vertical dispersion is allowed to be an
arbitrary prescribed function of the vertical coordinate. The horizontal
dispersion is assumed to be proportional to the 4/3 power of the plume
width. In Section 5.3, a model is developed for the case of a steady
release into a steady environment with a shear current (PTD). In
Section 5. 4, a model is developed for the case of a time varying release
into a time varying environment (UTD). In the latter model, the current
is unsteady but uniform.
Most of the models developed in this investigation represent generalizations
of previously existing models. For example, in Chapter 3, the previous
analyses on a single buoyant jet has been generalized to include a row of
jets which interferes and to include an arbitrary ambient density strati-
fication. In Chapter 5, the problem of dispersion from a continuous source
in a current has been generalized to include arbitrary vertical distributions
for current and vertical diffusivity. Previous models have assumed constant
current and constant diffusivities. In these cases, the solutions found are
as expected in the sense that they are qualitatively the same as those pre-
viously found. The model developed in Chapter 4, on the horizontal sur-
face buoyant jet, however, gives results which are qualitatively different
from previous investigations on either the ordinary jet or the submerged
buoyant jet. These results should be verified in the laboratory. Some
laboratory experiments on this phenomenon have been reported by Stefan
and Schiebe (1968), which showed some of the qualitative features found
in this investigation. These should be analyzed in more detail in the light
of the present findings.
6
-------
CHAPTER 2. CONCLUSIONS & RECOMMENDATIONS
En this report, several mathematical models have been developed for
predicting the distribution of excess temperature resulting from the
discharge of heated cooling water from power plants into large bodies
of water. The main conclusions and recommendations are as follows:
a) Initial mixing for subsurface discharge.
1. The flow field and mixing resulting from
a row of equally spaced round buoyant jets
discharging at an arbitrary angle into a
quiescent ambient with arbitrary density
and temperature stratification is formu-
lated. A computer program RBJ (Appendix A)
has been written to obtain the solution given
the jet characteristics and the ambient condi-
tions.
2. The flow field consists of two zones. Near
the source, the individual round jets behave
as if they were single jets. Further away,
they merge together and resemble a two-
dimensional slot jet.
3. The transition from one zone to the other is
taken to be either 1) when the round jet width
is equal to the jet spacing, or 2) when the
entrainment based on round jet analysis and
slot jet analysis are equal.
4. It is found that the two transitions give virtual-
ly identical results except in the small region
between the two transition points.
7
-------
5. Due to the relatively large dilution ratios and
the fact that the temperature excess of the
discharge is usually only 10 to 20°F, a very
small density stratification is sufficient to
prevent the discharge from reaching the sur-
face. In that event, all the temperature ex-
cess is assimilated in the ambient subsurface
water.
It is recommended that
1. A parametric study be performed based on
the model developed (RBJ) to obtain the re-
suiting temperature excess distribution in
a variety of cases.
2. The model be extended to include end effects.
The model developed assumes infinitely many
equally spaced round jets. Practical multi-
port diffusers are of finite length. Thus it
is recommended that the effects of the end
jets be analyzed. It can be expected that the
effects would be most important for short
diffusers. However, the model RBJ should
give conservative results.
3. The model be extended to include an ambient
current. Presence of ambient current would
further contribute to the mixing of the effluent
with the receiving waters. This effect should
be analyzed by incorporating the current into
the model guided by available experiments.
8
-------
4. The model be verified by laboratory investi-
gation. Although there has been much laboratory
investigation on single buoyant jets in linearly
stratified ambient, there is little on multiple
jets in non-linearly stratified water. Experi-
ments should be performed to verify the findings
from this model.
b) Intermediate phase for subsurface discharge.
In the event the buoyant jet reaches the free
surface, a model is developed to obtain
the spreading of the buoyant fluid on the
ambient (Sec. 3. 3).
In the two-dimensional case, the horizontal
extent of the spreading layer is found to grow
(after a brief initial period) at first linearly
with time t gradually becoming proportional
to
3. In the axisymmetric case, the radius of the
spreading layer is found to grow (after a brief
initial period) as at first gradually becoming
proportional to t 1 ” .
It is recommended that
A laboratory investigation be performed to
verify the model and to obtain the coefficients
needed.
9
-------
c) Surface horizontal buoyant jet.
The flow field induced and the dispersion
process in a horizontal buoyant jet discharged
at the surface is investigated in Chapter 4.
The interplay of source momentum, source
buoyancy, inte rfacial shear, entrainment,
and surface heat loss have all been in-
corporated in the model. The two-dimensional
case has been treated in detail.
The flow field is found to possess features
different from that of either an ordinary non-
buoyant jet or a fully submerged buoyant jet.
3. The type of flow field is found to depend on the
relative magnitudes of three parameters:
F , the source densimetric Froude number,
0
k, the dimensionless surface heat exchange
coefficient, and B, the source Reynolds
number.
4. The parameter space can be divided into three
regions such that given F 0 and P for example,
there exist two critical values for k, k > k
cr+ cr-
such that if k > kcr+ the flow field is of jet type.
If k < k , the source is inundated. For
cr-
k < k < k , the flow field consists of a jet
cr- cr+
type region near the source followed by an
internal hydraulic jump. The flow field after
the jump resembles that in a t’wo-laye red
stratified flow.
10
-------
5. The finding in 4. is of importance from
design considerations since; a) the tempera-
ture distribution is dependent on the type of
flow field, and b) inundation of the source
may lead to short circuiting the intake and
discharge of cooling water.
It is recommended that
1. Available laboratory experiments be analyzed
to verify the findings and determine the co-
efficient s.
2. Further laboratory experiments be performed
to supplement those already available.
3. The analogous axisymmetric problem( Sec. 4. 3)
be analyzed in detail to obtain the critical re-
lations in parameter space.
4. The analysis be extended to the three-dimensional
case.
d) Passive diffusion in a current.
I. Two mathematical models have been developed
in Chapter 5 resulting in two computer programs;
PTD (Appendix C) and UTD (Appendix D). In
both models, longitudinal dispersion is ignored.
Lateral dispersion is assumed to follow a 4/3
law. The vertical diffusion coefficient is allowed
to be an arbitrary function of the vertical co-
ordinate.
11
-------
2. Program PTD treats the case of steady
release of contaminant into a steady uni-
directional shear current. Program UTD
treats the case of a time-variable release
into a time-varying uniform current. The
surface heat exchange coefficient, assumed
constant in PTD, is allowed to be time
varying in UTD.
3. Results from the models are as would be
expected. Larger diffusion coefficients
result in larger dispersion.
4. The presense of current shear was found
to enhance the dispersion process.
It is recommended that
The models be applied to a variety of cases
to further obtain the quantitative dependence
of the dispersion process on the parameter.
2. The model be extended to include the effects
of longitudinal diffusion and a time dependent
current shear. Moreover, the current should
be allowed to exist in two directions, so that
the current direction as well as magnitude
are functions of depth.
The models should also be extended to include
the possibility of current reversals and times
of slack (no current). The extension of the
12
-------
models to include those effects can be
achieved by 1) formulating the dispersion
model for an instantaneous source in a
general ambient and 2) superposing the
solutions to obtain the resulting contaminant
distribution.
13
-------
CHAPTER 3 INITIAL MIXING AND SURFACE SPREADING DUE TO
SUBSURFACE DISCHARGE
3. 1 Introduction
The heated cooling water containing the waste heat from a power plant is
often discharged into a neighboring body of water. The discharge may be
made at the surface or submerged at some depth. These may be thought
of as representing two different philosophies of alleviating the thermal
pollution problem. If the discharge is made at the surface then it is
possible to design the discharge structure in such a way that minimal
mixing occurs between the effluent and the ambient water. This promotes
a comparatively high rate of heat dissipation to the atmosphere. On the
other hand, if the discharge is made at depth, then it is possible to promote
much initial mixing reducing the temperature rise in the ambient. In
this case, however, the rate of heat dissipation to the atmosphere is
correspondingly much lower. Of course, it is also possible to design a
surface discharge scheme which results in much initial mixing. The sub-
ject of surface discharge will be treated in Chapter 4 of this report where
it will be shown that proper design to achieve given goals may be more
difficult than it first appears.
In this chapter, we shall investigate some aspects of the flow and mixing
phenomenon when the discharge is made at depth. Since the heated water
is slightly less dense than the ambient water, the efflux would tend to rise
towards the surface. The general mixing and flow field can be conveniently
divided into three phases: an initial phase of jet mixing where the momentum
and buoyancy of the efflux are of importance in governing the flo’ . a final
phase of passive turbulent diffusion where the ambient turbulence and currents
are dominant in the dispersion process; and an intermediate phase joining
the two. The final phase of passive turbulent diffusion will be investigated
in Chapter 5. In this chapter, we shall investigate the initial phase of jet
mixing and the intermediate phase.
15
-------
3.2 Initial Mixing Phase: Multiple Submerged Buoyant Jets Discharged
into an Arbitrarily Stratified Ambient
The mixing phenomenon in buoyant jets and plumes have been studied by
numerous investigators. Morton, Taylor and Turner (1956) applied an
integral method to the problem of a buoyant plume discharged as a point
source into a linearly stratified ambient fluid. Brooks and Koh (1965)
analyzed the two-dimensional buoyant jet problem with application to the
design of a submerged ocean outfall diffuser. Fan (1967) examined the
more general case where the angle of discharge is arbitrary. Abraham
(1963) examined the same problem as Fan but using a slightly different
mechanism of entrainment. All these studies assume 1) similarity in the
flow pattern, 2) density differences are small so that the Boussineq
approximation is valid, 3) an entrainment mechanism which depends only
on the local mean flow, and 4) the ambient fluid is motionless.
For application to the practical situation, it is necessary to generalize
the models developed in several respects. First, all the previous studies
are for a single jet (either axisymmetric or two-dimensional). In a
practical case, there may be many jets spaced by a certain distance from
one another. Thus after an initial period, these jets would merge and
interfere with each other. Second, the ambient density stratification is
most likely not linear. Third, the ambient fluid may not be motionless
in the vicinity of the discharge. Finally, for application to thermal
discharges, it must be realized that the ambient density stratification
and the thermal stratification may be different as, for example, in the
case when salinity variations are also contributing to the density strati-
fic a ti on.
In this chapter, a general integral model is first developed which is
independent of the geometry of the jet (i. e. , equally applicable to two-
dimensional or axisymmetric cases). The ambient density stratification
and temperature stratification are considered as independent and arbi-
trary (may be nonlinear). This model is then specialized to either the
16
-------
two-dimensional case of a slot jet or the axisymmetric case of a round
jet. The case of a row, of equally spaced round jets, is then examined.
Near the source, before jet interference, the individual round jets are
more or less separate and the axisymmetric case applies. At some point,
the individual jets would begin to merge and finally the two-dimensional
case would apply. The transition from one to the other occurs in an
intermediate zone where neither axisymmetry or two-dimensionality
obtain. In the model to be developed, the transition is taken to be sudden.
However, two different transition criteria were used and the results are
found to be virtually the same based on either criteria.
3.2. 1 Formulation
Consider a jet oriented at angle 0 to the horizontal issuing fluid of
density 01 and temperature T 1 into an ambient of density stratification
and temperature stratification T(y). Let be the discharge,
M 1 the momentum flux, F 1 the density deficiency flux, and the tem-
perature deficiency flux at the source. Figure 3. 1 illustrates the general
behavior of such a jet. The bending of the jet path is a result of the fact
that the discharge is buoyant. Define u as the velocity, T the temperature
and the density of the fluid. Since the ambient is motionless, u
is assumed to be along the jet path. Let s be the coordinate along the
jet path, A-plane be the plane perpendicular to the jet path, and the
angle of the jet path with respect to the horizontal. We now define the
volume flux Q, momentum flux M, density deficiency flux F and tem-
perature deficiency flux G as
Q = S u dA, (3. 1)
M = - = 5 u 2 pdA 5 .u dA (3. 2)
0 °A A
F = 5 a p)u (3.3)
Note that M’ is the true momentum flux while M is the kinerratic momentum
flux.
17
-------
y
y y
Figure 3. 1 Definition sketch
Figure 3.2 Jet interference
S
U.
A- plane
Do/\
p 0
18
-------
G = (Ta - T) u dA (3.4)
We note that in the vertical direction, a buoyancy force exists due to the
density difference between the jet fluid and the ambient fluid tending to
bend the jet. This force f is
= g $ a - p) (3. 5)
The conservation equations can now be written in terms of the variables.
The conservation of mass equation is
E (3.6)
ds
where E is the rate of entrainment of ambient fluid. Note that strictly
speaking, since the density is variable we should really have, instead of
Eq. (3.6)
} = Epa (3.7)
However, all density differences are small and we may approximate the
Eq. (3.7) by (3.6)
The conservation of horizontal momentum flux is
d(M cos 0 ) 0 (3.8)
ds
For the vertical momentum, we must include the buoyancy force. Thus
( sin 0) = f (3.9)
The conservation of density deficiency flux equation reads
{ $u(p 0 - p) dA} E(p 0 - (3. 10)
19
-------
where p is a reference density (e. g. p 0 = Equation (3. 10)
can be written
{Su po dA + $u (p dA}=E(p p)
or
dQ dp dF
- + C ds + = E(p (3.11)
Using Eq. (3. 6), we finally have
dp
dF_ aQ 312)
dsds
Similarly, the conservation of temperature deficiency flux equation reads
- { J’u (T - T ) dA} = E(T - Ta) (3. 13)
which reduces, with Eq. (3.6), to
dT
dG_ aQ 314)
ds ds
Equations (3. 6), (3. 8), (3. 9), (3. 12) and (3. 14) constitutes five equations
for the five unknowns Q, M, e, F, G, as functions of s once we can
express E and f in terms of known quantities or these unknowns. To do
this we will make two more assumptions. First, we shall assume similarity
of the shapes of the velocity profile, temperature deficiency profile and
density deficiency profile in the plane A. In particular, it will be assumed
that the profiles are Gaussian. Thus, in the two-dimensional case (slot
jet) we assume
22
* - /b (s) (3. 15)
u (s , ) u(s) e
2 22
-r /X b
- p = a - p (s)] e (3.1 )
20
-------
- /X 2 b 2 (3.17)
Ta T(s, ) = [ T - T(s)] e S
where u(s), p(s), and T(s) are the values along the jet centerline. is
the coordinate normal to s, b(s) is the characteristic jet width, and
is a turbulent Schmidt number for the two-dimensional case. Similarly
in the axisymmetric case, we take
22
- r /b
u (s,r) = u(s) e (3. 18)
2 22
- p (s, r) a - p(s)] e r (3. 19)
2 22
-r/Xb
Ta - T(s,r) [ Ta - T(s)] e r (3. 20)
Secondly, we shall assume that the entrainment function E is proportional
to the jet characteristic velocity u and the jet boundary (2 rrb or ZL) and
the proportionality constant is a. r for round jet and for slot jet).
Substituting these expressions into the definitions for Q, M, F, G and f
(Eqs. (3. 1) through (3. 5)) will give Q, M, F, G and f expressed in
terms of the quantities u, p , T and b. For example, substituting Eq.
(3. 15) into (3.1) gives
Q LSu(se /bdfl
where L is the length of the jet. Thus
Q = Lub$ e d = ubL
Similarly, substituting Eqs. (3. 15) and (3. 16) into Eq. (3.3) gives
2 22 2 2
-il/Xb -i 1 /b
F = U u e S a - p) d
- -( 2 /b 2 ) (1/X 2 +l)
= LU(Pa - p) j e d
-
fz
(x
= LU(Pa P)bJ
1+x
5
2.1
-------
In this fashion, Table 3. 1 may be constructed;
TABLE 3. 1
For Round Jet
For Slot Jet
of Length L
Volume Flux Q
Tub 2
ubL
Momentum Flux M
Density Deficiency Flux
F
uu 2 b 2 /2
x2
r ub 2 (p -p)
a
r
jyr/2 u 2 bL
/ ub(p p)L
j a
S
Temperature Deficiency
Flux G
2
x
nub 2 (Ta T)
1+X
r
[ 2
/TIX
/ S 2 Ub(Ta T)L
\‘ l+X
S
Buoyancy Force f
X b 2 a - p)g
5 bLg(o p)
Entrainment Function E
2rra ub
r
Z X uL
S
This table gives the transformation from the variables u, b, etc.
to the variables Q, M, etc. The inverse transformation is given
in Table 3. 2. Moreover, it is possible to express E and f now in
terms of Q, M, etc. as shown in Table 3. 3
The problem under consideration in this chapter is the mixing pro-
cesses involved for a row of round buoyant jets spaced a distance L
apart. initially, the jets are separate round jets. However, after
a while, they begin to merge and form more nearly a two-dimensional
-------
TABLE 3.2
Pound Jet
Slot Jet
Centerline Velocity u
2M/Q
\Z M/Q
Norminal Half Width b
Q/ J2TiM
Q 2 / [ f LM]
Density Deficiency p—p
2
— F/Q
x
r
Dilution Ratio S
Temperature Deficiency
T-T
a
Q/QJ
+ 2
r G/Q
r
Q/Q 1
2
/ 2 G/Q
V
S
TABLE 3.3
Pound Jet
Slot Jet
E
2 / a JM
2 J LM/Q
f
i+x 2
gQF
/ + 2 QgF
, 2 M
slot jet (see Fig. 3. 2). Thus in the calculations, it is necessary to
provide a criterion whereby the round jet analysis is switched to that
for a slot jet. Two such criteria are proposed. First, we may
assume that transition occurs when the width of the round jet becomes
equal to the jet spacing. This shall be designated transition 1. Re-
ferring to Table 3. 2, this occurs when
Q2 J/2 L or q_= L = O.885L
f rrM 2
(3. 21)
23
-------
Here the jet width” has been taken to be 2 j2 b. Alternatively, we may
assume that transition occurs when the entrainment as calculated by the
round jet theory or the slot jet theory becomes equal. This shall be de-
signated transition 2. Referring to Table 3. 3, we see that this occurs
when
2/ a / = 2 JZa LMIQ
or
a
= L
Since experimental values for a and a.r are 0. 1 6 and 0. 082 respectively,
this is
Q = 0. 16 L = 1. 1 L (3. 22)
0. 082 \/1-r
This occurs somewhat later than the first transition, Of the two transition
criteria, it is felt that the second is more reasonable since, in the first
one, it is necessary to define the “jet width” which is somewhat arbitrary.
It has been found, however, that the two criteria gave virtually the same
results except for the region between transition 1 and 2. Thus the solution
is not sensitive to exactly where the transition is.
It should be noted that the independent variable of integration is s, the
distance along the jet path. However, the ambient conditions a and Ta
are usually only given as functions of y. Thus the following two equations
are needed to allow conversion between s and x, y:
= cos e (3.23)
ds
= sin 8 (3.24)
ds
The system of Eqs. (3. 6), (3. 8), (3. 9), (3. 12), (3. 14), (3. 23) and (3. 24)
constitute seven ordinary differential equations for the seven unknowns,
24
-------
Q, M, G, F, G, x and y as function of s. These equations may be solved
given the initial values of the unknowns at s = 0.
The initial conditions are given by the source conditions, namely, u,
the jet velocity, D 0 , the jet diameter, T 1 , the jet temperature, p 1 , the
jet density, e 0 , the jet discharge angle. However, since the formulation
is in terms of the flux quantities Q, M, F, and G, these jet characteristics
must be converted to initial values in these variables. Moreover, it is
well known that there exists a zone of flow establishment extending a few
jet diameters during which the top hat profiles of velocity, density de-
ficiency and temperature excess change gradually to Gaussian form. In
this formulation, we shall start the integration from the beginning of the
zone of established flow. Thus it is necessary to relate the jet charac-
teristics to the flux quantities at this point. Albertson, et. al (1950), in
their experimental investigations on the round jet found that the zone of
flow establishment extends a distance of 6. 2 jet diameters. Equating the
momentum flux at the beginning and end of the zone of flow establishment,
(assuming that the buoyancy force is negligible in such a short region),
we get
22
rrb U
i-i 22 2 00
u = u 2 rdr =
400 2
0
Thus the initial value for Q is
2 2
Q = irb u = —D u
1 o o 2 o 0
In other words, the volume flux at the beginning of the zone of established
flow is twice that at the source.
By assuming further that the ambient density is uniform in the zone of flow
establishment, we may equate the density deficiency flux at the beginning
and end of this zone to get
25
-------
- p 1 ) = $u(p - p) 2 rdr
x 2
r 2
= 2 ub Uo(Pa P)
1+x
r
Thus the initial value for F is
2
F 1 = D 0 Uo(Pa - p 1 )
However, the centerline density deficiency is
2
- p a - p 1 )
For X = 1.16, p 0.87 a -
Similarly,
G = D 2 u (T -T
1 400 a 1
and the centerline temperature excess is
2
T - T = r (T - T ) = 0.87(T - T
a 2X 2 a 1 a 1
The temperature excess and density deficiency are thus already decreased
to 87% of their values at the jet efflux due to their faster spreading. In
other words, the zone of flow establishment for temperature and density are
shorter than that for velocity. Since we are starting the calculations at
the end of the zone of flow establishment for velocity, the temperature
excess and density deficiency have already undergone some decrease,
namely 13%. Note that if X = 1 (same spreading and length of zone of
flow establishment) then this decrease would be absent. This phenomenon
is shown schematically in Fig. 3. 3.
26
-------
U 0
ZONE OF FLOW
ESTABLISHMENT
Figure 3.3
Zori e of flow establishment in a submerged jet.
-------
3. 2. 2 Method of Solution and Examples
A computer program I BJ has been written in Fortran IV language to solve
the problem formulated in the previous section. A listing of the program is
included as Appendix A. The method of solution is to first obtain the initial
conditions Q 1 , M 1 , F 1 , G 1 , from the given source characteristics. Then
the Eqs. (3. 6), (3. 8), (3. 9), (3. 12), (3. 14), (3. 23) and (3. 24) are integrated
with E and f as given by those for the round jet (column 2, Table 3. 3).
When transition is reached, given by either Eq. (3. 21) or (3. 22), one simply
continues the solution but with E and f as given by those for the slot jet
(column 3, Table 3. 3). The results obtained are then converted from the
variables Q, M, F, G,.. . to the physical variables u, p, T, and W, the
jet width which is takeii to be 2. .j2 b. The conversion is effected by the
relations in Table 3. 2.
The inputs to the program consists of the following:
u jet velocity
D = jet diameter
T 1 = jet temperature
p 1 = jet density (in gm/c. c.
jet discharge angle with respect to the horizontal,
(in degrees)
d. = jet discharge depth
L = jet spacing
a = entrainment coefficient for a round jet (usually
r taken to be 0. 082) *
a = entrainment coefficient for a slot jet (usually taken
S to be 0. 16) *
X = spreading ratio for a round jet (usually taken to be
r 1. 16)
= spreading ratio for a slot jet (usually taken to be 1. 0)*
g = gravitational acceleration
The program is written in such a fashion that all the quantities are dimen—
sional. However, any consistant system of units may be used (FPS, MKS,
or CGS), except that density is always in units of gm/cc. The specification
* See Fan (1967)
28
-------
of g, the gravitational acceleration is utilized as the indicator for the units.
Thus, for example, if g is specified to be 32. 2, then the system of units
which should be used is the FPS system.
In addition to the inputs above, it is also necessary to specify the density
and temperature profiles in the ambient. This is accomplished by speci-
fying a table of depth, temperature, and density. The program linearly
interpolates between specified values to arrive at the values for inter-
mediate points. For example, if the ambient is linearly stratified in
temperature and density, then only two points need to be specified one at the
surface, and one at the jet depth. When the two values coincide, then the
ambient is uniform.
Output quantities from the program consist of x, y, jet width, dilution,
jet temperature, jet density, ambient density, ambient temperature, and
temperature excess. The quantity jet width is taken to be 2A/ b. Dilution
is the ratio of volume flux Q to that at the beginning of the zone of established flow.
Jet temperature, density and temperature excess are all centerline values.
These should be kept in mind when interpreting the results.
Example cases have been solved using the program RBJ and these cases
are summarized in Table 3.4. The solutions are shown graphically in
Figs. 3.4 to 3.6. The effect of the various parameters can be readily
seen from the figures. Figure 3.4 shows that the jet path for various
value of L, the jet spacing and D, the jet diameter in a uniform environ-
ment. It is readily seen that, as expected, decreasing D 0 or increasing
L, implying delayed jet intereference, bends the jet upwards. Figure 3. 5
shows the jet excess temperature as a function of the vertical coordinate.
It is seen that the transition point is not of import except for a short zone
between the two transition points. Otherwise, the predictions of excess
temperature based on the two transition points are virtually identical. It
should be noted that in Fig. 3.5, the scale on y starts at y = 1. The jets
have already travelled horizontally quite a distance before reaching y = 1
(see Figure 3.4). Thus they have already achieved a significant reduction
in T through entrainment. Also note that in most cases in Figure 3.5,
transition occured before y = 1. Since the excess temperature at discharge
is usually only a matter of 10 or 2 0°F, and since the dilution is usually
quite large in a submerged jet, only a very small density stratification is
needed to suppress the jet from reaching the surface. This can be seen
29
-------
in Fig. 3. 6. In these cases there is a temperature difference of but 0. 8°F
over a distance of 60 ft. , yet when the jet temperature was 82. 4°F (jet
efflux excess temperature = 12. 2°F), the jet never even reached a vertical
distance of 60 ft. When the jet temperature is 89. 2°F (jet efflux excess
temperature = 20°F), the jet does reach above the 60 ft. level. However,
it is stopped by the the rmoc line. It can be expected that if the ambient is
stratified, even very slightly, it would not be difficult to design an outfall
diffuser to always keep the discharge submerged.
TABLE 3.4
D u T d.
0 0 1 j
1
0.5
12.5 89.2
2
0.5
12.5 89.2
3
0.25
12.5 89.2
4
1
12.5 89.2
5
1
12.5 89.2
6
1
12.5 89.2
Non-uniform
Ambient
7
0.5 •
12.5 89.2
100
8
0.5
12.5 82.4
100
L
Ambient
Ta (y)
100
5
uniform
77°
100
10
.
uniform
0
77
100
5
uniform
77°
100
5
uniform
77°
100
10
.
uniform
0
77
100
20
.
uniform
0
77
0
20
40
60
5 80-
77°F
0°F
100 ____ 69 .2°F
30
-------
3 2
I 6
4-
-J
U i
I-
Li
0
L, ) Ui
— C.)
z
C,)
-J
C-)
I-
Ui
>
x : HORIZONTAL DISTANCE OF TRAVEL, ft
Predicted trajectories of multiple buoyant
Figure 3. 4
6 : 00, U 0 : 12.5 fps, p 1 :0.99581 gm/cc
Pa:O 9 9 7 O 7 gm/cc, T 1 :89.2°F, Ia: 77°F (UNIFORM)
5
tOO
50
0
0 50 tOO 150
jets in a uniform environmentS
-------
D 0 :lft, U 0 :I2.Sfps
P 1 0.995 18gm/cc
y, VERTICAL DISTANCE OF TRAVEL, FT.
Figure 3. 5a
Prt dicted jet centerline excess ten perature of
multiple buoyant jets in uniform environment.
L: 5’
0
w
I-
4
w
0
LU
I—
C l)
Cl)
LU
C.)
LU
I-
l.0
0.I
=20’
6 =00
Transition I
—— Transition 2
T 1 89.2°F
a O.997O7gm/cc
T 0 !77°F
100 10 I
-------
U-
0
w
w
0 .
U i
F-
Li) Q)
C l )
LU
0
x
LU
I .-
1.0
0.I
I I I
y, VERTICAL DISTANCE OF TRAVEL, FT.
Figure 3. 5b
Predicted jet centerline excess temperature of
8 =o°, U 0 = 12.Sfps
0.99707gm/cc
T 1
89.2°F
p 0
T 0
0.99518gm/cc (uniform)
= 77°F
(uniform)
L=5 , D 0 =0.
L=5’, D 0 =O.25’
100 10 1
multiple buoyant jets in uniform environment.
-------
X, HORIZONTAL DISTANCE OF TRAVEL,ft.
Figure 3. 6
Predicted trajectories of multiple buoyant
100
8
60
‘I -
-J
Lii
I -.
Li.
0
LU
0
z
4
I-
(I )
-J
4
0
L ii
>
40
20
jets in stratified environment.
-------
3. 3 Time Dependent Surface Spreading of a Buoyant Fluid
When the warm efflux discharged at some depth reaches the surface of the
ambient fluid, it may still possess some buoyancy and will thus spread on
the surface. The phenomenon of surface spreading can be likened to that
of a surface horizontal buoyant jet, which is treated in Chapter 4. There
it was found that if surface heat exchange is absent, no steady state solutions
may be found. In this section, the unsteady surface spreading problem will
be investigated. The analysis to be presented is very approximate and
should be viewed as providing only an order of magnitude estimate of the
phenomenon. No detailed flow field will be derived. Only the gross
properties of the spreading pool of buoyant fluid will be obtained. The
analysis further incorporates several coefficients on which no data is
available. Experiments on this phenomenon should be performed in the
future to verify the findings and provide estimates of the coefficients.
In this investigation, no surface heat exchange or entrainment of ambient
fluid will be included. In Chapter 4, it will be found that it is when surface
heat exchange is absent that a steady state solution cannot be found.
Moreover, the spreading layer thickens with time so that after an initial
period, entrainment may also be ignored.
3. 3. 1 Two-Dimensional Case
We assume that at time t 0, a line source of strength 2 q 0 per unit
length injects lighter fluid of density p onto a heavier quiescent ambient
fluid of density p as shown in Fig. 3. 7. We make the following assumptions:
1. no entrainment of ambient fluid occurs;
as the buoyant fluid spreads, the shape of the interface is
similar from one instant to another; and
3. the pressure distribution is hydrostatic.
Under these assumptions, we now examine the motion of the pool as a whole.
Consider a half of the spreading pool at time t as shown in Fig. 3. 8. For
simplicity, the similar shape will be taken to be rectangular. It will become
35
-------
2 q 0
jj 2 b
b
Figure 3. 7 Definition sketch.
__ Jja
b- -
00
Figure 3. 8 Definition sketch.
36
-------
obvious later that taking it to be some other shape will not change the basic
features of the resulting equations but will only change some of the numerical
coefficients.
We now write the equation of motion for the buoyant fluid as a whole. The
time rate of change of the momentum of the center of mass is
d l
- i pa ( b) b
where the symbols are as defined in Fig. 3. 7. The driving force is the
pressure difference induced by the density difference. It can easily be
demonstrated that this is
I Z
g( p) a
where L\p = p 0 - p. The resistive forces of shear and hydrodynamic
drag are respectively
0
and
C p(b a
Now we assume T = b— where e is an effective viscosity coefficient.
The equation of motion of the spreading pool is then
[ abbt] = g( p)a 2 - CDp (b)Z - -b’(b - b) (3.Z5)
For Ap << p, Eq. (3.25) reduces to
1 C (b-b)
abe-] = g’ a 2 - - a(b’) 2 - a b’ (3. 26)
where g’ = - g . Now, for no entrainment, we have the continuity relation
37
-------
ab = qt (3.27)
Hence Eq. (3. 26) can be written
22 2
d db q t CD qt r db 7 \ b(b-b) db
2 - - - { qt - - 11 — 2 g b 2 - 2 T i1 - 2 qt - - (3. 28)
This equation can be solved for b as a function of time t for given initial
conditions. The terms in Eq. (3. 28) represent, respectively, the local
inertia, the pressure driving force induced by the density difference, the
hydrodynamic drag and the shear.
It can be expected that after a brief initial period (for which this analysis
d db
is probably not valid), the inertia term - - (q t —) probably becomes
negligible. This initial period is followed by one when the hydrodynamic
drag is balanced by the driving pressure with the shear of secondary im-
portance. Then the equation is
2 C
g’a — -a(b) 2 (3.29)
Note that this implies
7
dti = (3.30)
g’a CD
The left hand side is simply a densimetric Froude number. In studies
of density currents (such as turbidity currents and cold fronts), it has been
found that this Froude number is constant and equal to approximately 2
corresponding to CD = 0. 5, a very reasonable number. Moreover, this
equation gives rise to the solution
b [ Zg’ q 0 ]hI3t (3. 31)
where upon
38
-------
2/3
qt q q 0
a = — — = 1/3 = 1/3 = constant (3. 32)
(Zgq) (Zg’)
For large time, it can be expected that the shear term dominates as the
resistive mechanism so that the pertinent equation is approximately
2 €bdb
g a = - - j-j -- (3. 33)
With ab = q 0 t, we get
b 4 - - = q 0 3 t 3 (3. 34)
with solution
3 1/5
r 5pg q 0 —
b = L 8€ i (3.35)
and
I
r- 5 pgq 0 -,
a = L 8€ J q 0 t (3.36)
In summary, b, the horizontal extent of the spreading pool, increases
(after a brief initial period) linearly with distance while the thickness is
essentially constant. When the spreading has proceeded far enough for
viscous effects to dominate, the spreading rate decreases to while
the thickness increases hut slowly (11/5). As time t , the thickness
would tend to infinity.
We now note that if the similarity shape is not a rectangle but is some
other shape, the same dependence on time will be found. Only the nu-
merical factors in the proportionality constants will be different.
39
-------
3. 3. 2 Axisymmetric Case
We now derive the equation for the axisymmetric case. Basically this is
the same physical phenomenon and hence the same assumptions will be made.
Thus, instead of a line source as in the two-dimensional case, a small
circular source is taken to be emitting buoyant fluid onto a heavier motion-
less ambient as shown in Fig. 3. 9. Under the same assumptions as in the
previous section, we consider a slice of the spreading pool at time t as
shown in Fig. 3. 10. Instead of assuming a rectangular cross-sectional shape,
we shall assume an illipsoidal shape. The driving pressure force on the
section shown in Fig. 3. 10 can be found by integrating the pressure, thus
= -p) dA
where p is the pressure in the buoyant surface fluid and p that in the
ambient and A is the projected area over which the pressure acts (see
Fig. 3. 10).
Since pressure at the interface must be single-valued, it follows that
p = pga-Ogx 0
-------
Qo
Figure 3.9
b—
Figure 3.10
Definition sketch
Ut
Definition sketch
2b 0
p 0
iii b
p 0
bd8
A a
0
p 0
41
-------
For << 1, we have
p 0
2
= ( p)gb dB (3. 37)
The time rate of increase of momentum of the slice is
d I TTab 2 db
16
The shear and hydrodynamic drag forces are respectively
(b 2 - b 2)
o db
a dt
and
CD % ( - )de
so that the equation of motion is
2 2
d I TTab 2 db 4 p 2 ab db 2 ( b -b 0 db
16 i -J = ga b - CDp4( ) - a (3.38)
The continuity equation is
a b 2 = [ - Qt (3. 39)
Equations (3. 38) and (3. 39) can be solved for given parameters and given
initial conditions.
Again, if we assume that after a brief initial period, the local inertia
is negligible, so that the driving force is balanced by the hydrodynamic
drag, we get
4g’ a 2 b = CD (b’) 2 (3. 40)
where g’ =
42
-------
or
(3. 41)
is again a densimetric Froude number. Since the continuity
[ -}Qt, we have
3C g’ [ f]Qt
2
a. , we get
b
if at time
( b’) 2 - 4
ag’ - 3 CD
The left hand side
equation is a b 2 =
2
A I Zi-i
L3
Letting 3C
D
a Jg ’Q Jt (3.42)
t = O,b=b,therefore the solution to Eq. (3. 42) is simply
b = ‘b + t (3.43)
Thus, after a brief initial period, the diameter of the spreading pooi grows
3/4
as t . The thickness a is
Zn
a = _____
____ [ ]Q
______ 0 (34)
-) 4a Jg ’Q
°
0 3
which first increases and then decreases with time.
We next examine the solution for large time when the dominating resistive
force is the shear. In this case, the Eq. (3. 38) reduces approximately to
j_a 2 b = ( - ) -- - - - (3.45)
43
-------
3
3 2rr
gQ T- 1
Using the continuity Eq. (3. 39) and letting cx = 0 , we get
3 7db
a 1 t = b
with solution
b (2aj)8t,1 (3.46)
The thickness a is then
1 2rr 1
LTJQOt _____
a = 2 = = constant (3.47)
b (Zct 1 ) 4
Thus, after a brief initial period, the diameter of the spreading pool
would grow as gradually decreasing to t 1 while the thickness first
increases and then decreases tending to a constant value
(ZcL 1 ) 4
3. 3. 3 Comparison with Experiments
The theory presented in the previous section on the axisymmetric time
dependent surface spreading can be compared with the experiments by
Sharp (1969). Sharp reported on the growth of the radius of the spreading
pool resulting from the surfacing of buoyant jets discharged at the bottom
of a laboratory tank. These results are summarized in Table 3. 5.
44
-------
TABLE 3.5
1/5 I
b(g ) tjgj
1 s
Qo Q
0
1.5 1
2.3 2
3.4 4
4.3 6
5.2 8
5.4 10
8.8 20
11.9 30
14.4 40
16.5 50
For b << b, Eq. (3. 43) (viscous effects negligible) can be written
b = (4 )2 (g’Q) 4
or putting it in Sharp t s notation
1 3/4
b(g I)h/5 4a 2 I t(g) 3 48
Q 2 /S — 3 ‘t. Q1/5 I
Figure 3. 11 shows Sharp t s data together with Eq. (3. 48) with a = 3/4.
It should be noted that for small t and b, the influence of b becomes
0
more important and one expects higher measured values of b than given
by Eq. (3. 48). Since Sharp did not detail the initial values b 0 for his
various experiments, no estimate can be given as to its influence. It is
clear however, that the agreement between Eq. (3. 48) and Sharp’s data is
reasonable.
45
-------
$ (91)3/5
, 1/5
wo,
Figure 3. 11
Comparison of theory with exp rtrnents tor tflt
surface spreading of buoyant fluid.
I A
— N 0
I0
I 10
100
-------
3. 3.4 Numerical Solutions
The equations for surface spreading derived in the previous sections can
be solved numerically including all the terms. For example, in the
axisymmetric case, the Eqs. (3. 38) and (3. 39) can be combined and by
further normalizing the variables by letting t* = tIt 0 , = b/b, where
to ( g ) it can be written
d 2 - + C 1 -( - -j - --
— — t dt’ 3 — D dt- dt ’
3b
8 o
where C D = CD, K = ZuQ 0 t 0
This equation has been solved for the cases C 1 D = 0. 5, and K = 0,
-3 -z
10 , 10 and are reproduced here as Fig. 3. 12.
47
-------
I
I
t, DIMENSIONLESS
Figure 3. 12
Growth of a spreading surface pool of buoyant
102
C O.5
K:O
K
K
K0 1
I0
I0
TIME
f1uid
-------
3.4 Summary and Discussion
In this chapter, the initial and intermediate phases of mixing due to the
subsurface discharge of warm cooling water is investigated. A mathe-
matical model is developed to describe the initial phase of mixing. In
this model, a row of round jets equally spaced at some distance apart is
allowed to emit the warm water at depth into an ambient fluid which may be
stratified in an arbitrary manner. The model is based on an integral
approach similar to t he ones used by previous investigators. The phe-
nomenon of jet interference is incorporated in the model. A computer
program RBJ based on this model is presented in Appendix A and can be
used to obtain the solution in any practical case. It is found that very little
stratification is enough to suppress the effluent from reaching the surface.
The intermediate phase of surface spreading of a buoyant fluid on top
of a heavier ambient is discussed in Section 3.3. The analysis is very
approximate and is performed primarily to obtain the time and length
scales of the problem. The results of the analysis can be used to pro-
vide a link between the buoyant jet portion of the phenomenon discussed
in Section 3.2 and the passive turbulent dispersion portion discussed
in Chapter 5. It should be pointed out that application of the analysis
requires the knowledge of several coefficients which can only be obtained
by experiments. These should be done in the future. Based on typical
values for the parameters, it can be inferred from the analysis that
the time scale of this intermediate phase is on the order of minutes.
49
-------
CHAPTER 4 SURFACE HORIZONTAL BUOYANT JETS
4. 1 Introduction
In this chapter, we shall investigate the behavior of warm (and hence
buoyant) jets discharged horizontally at the water surface. The receiving
water is assumed stationary and uniform in density. The two-dimensional
case of a slot jet is analysed in detail. The axisymmetric case is also examined.
The phenomena of entrainment, source buoyancy and momentum, inter-
facial shear and surface heat exchange and all included in the model.
A number of interesting results are found. In particular, it will be seen
that the behavior of such jets can be drastically different from either an
ordinary nonbuoyant jet or a fully submerged buoyant jet.
The behavior of ordinary nonbuoyant jets has been studied quite exten-
sively. For example, Albertson, et al. (1950) have performed a series
of laboratory experiments on both slot jets and round jets. The detailed
results will not be repeated here since they are well known. It is found
that the flow field can be conveniently divided into two zones: the zone of
flow establishment near the source where the finite size of the source is
important followed by the zone of established flow where only the source
momentum flux is of importance. In the zone of established flow, the
velocity distributions were found to be similar from one station to another
with the shape well approximated by a Gaussian profile. It was also found
that the momentum flux stays constant with distance downstream while
the mass flux increases with distance downstream due to entrainment of
ambient fluid.
Investigations into the behavior of submerged buoyant jets and plumes
have also been studied quite extensively being stimulated by practical pro-
blems in engineering and meteorology, and more recently by the advent of
multiport submerged sewage outfall diffusers. Such studies have typically
employed an integral approach assuming similarity of velocity and buoyancy
profiles and assuming a certain entrainment mechanism. The primary
effect of the jet buoyancy is in supplying an added force so that the flux
of vertical momentum is no longer constant as was the case in nonbuoyant
jets but is related to the buoyancy force. T ius, in the case of a horizontal
51
-------
buoyant jet, the jet path is deduced to bend upwards due to this buoyancy.
These studies include Abraham (1963), Brooks and Koh (1965) and Fan
(1967), among many others.
An analysis of the problem of a row of submerged buoyant jets discharging
into an ambient fluid with an arbitrary density stratification is presented
in Chapter 3 of this report.
In the problem to be considered in this chapter, the buoyant source is
situated at the surface and is discharging horizontally. For a source
with sufficiently strong initial momentum, it is expected that near the
source the behavior might resemble that encountered in ordinary sub-
merged jets. However, further away, as the momentum diffuses through
jet mixing, the influence of the buoyancy would manifest itself in modifying
the horizontal momentum making this problem fundamentally different
from the submerged buoyant jet analysed heretofore. Another way of
visualizing the difference is as follows. If a source of pure buoyancy
exists at some depth, the resulting flow field is primarily vertical towards
the water surface. However, if a source of buoyancy exists at the water
surface, the resulting flow must be horizontal. The driving force hori-
zontally is due to the buoyancy which modifies the pressure distribution
which in turn provides the driving force for the horizontal spreading.
Another important point of departure between ordinary submerged buoyant
jets and the surface buoyant jet considered herein is that in the present
case, one needs to include in the formulation the effects of interfacial
shear. Some distance away from the source, after the momentum has
diffused, the buoyant fluid tends to simply float on the denser ambient.
Shear forces at the interface then play an important role in the dynamics
of the flow. This mechanism is not of import in ordinary submerged
jets since in that case, the flow belongs to the class of free turbulent flows
with typically Gaussian velocity profiles.
In the following sections of this chapter, the two-dimensional case of a
warm jet discharging horizontally at the surface of an infinite body of
52
-------
water will be investigated in detail. The interplay of source buoyancy,
source momentum, entrainment, interfacial shear and surface heat
exchange will be analysed. It will be demonstrated that if the mechanism
of surface heat exchange is absent such as in the case when the density
difference is induced by salinity, no steady state solution exists. The
source will be inundated by the discharge and the depth of inundation will
increase with time. When the mechanism of surface heat exchange is
included in the formulation, it is found that a steady state solution always
exists. It is further found that the general flow field may possess quite
different features depending on the relative importance of the various
parameters. For example, if the surface heat exchange coefficient K
is sufficiently large, then the flow field may resemble that in an ordinary
jet. On the other hand, for small values of K, the jet may not be able
to persist resulting in an internal hydraulic jump followed by a zone where
the flow is essentially that of a two layered stratified flow. The location
of the hydraulic jump is dependent not only on the source conditions but
also on downstream conditions which in this case of an infinite fluid is
replaced by the surface heat exchange and interfacial shear coefficients.
Under certain conditions, the source may be inundated and the zone of
stratified flow extends all the way to the source. However, with surface
heat exchange, a steady state still exists. The depth of inundation is then
governed by downstream conditions and is quite independent of the source
momentum.
The general flow field can thus be divided into several zones within each
different mechanisms dominate, keeping in mind that under certain conditions,
not all the zones may be present. At the source, the source momentum
may dominate and the flow field is like a jet. However, the buoyancy
reduces the entrainment rate so that the rate of increase of the jet thick-
ness decreases. Far from the source, interfacial shear becomes more
important and the flow field becomes controlled by downstream constraints
such as by a tailgate in a laboratory tank. In between, the flow field may
go through an internal hydraulic jump. In the case with heat exchange at
the water surface, this mechanism plays the role of downstream control
and plays a part in determining the location of the internal hydraulic jump.
53
-------
Although only the two-dimensional case is reported in detail here, the
general qualitative features of the flow field should remain valid for other
geometries. For example, in the axisymmetric case, one would expect
the possibility of a circular internal hydraulic jump. The equations and
some solutions for this case are also included in this chapter although
it has not been carried to the same amount of detail.
Previous investigations of related problems also include Wada (1966),
Lean and Whillock (1965), Hayashi and Shuto (1967), Jen, Wiegal and
Mobarek (1966) and Stefan and Schiebe (1968). Wada (1966) and Hayashi
and Shuto (1967) investigated theoretically the temperature distribution
when warm water was discharged from a rectangular outlet at the surface
of a semiinfinite motionless ambient. The inertia of the fluid was ignored
and the temperature distribution is the result only of the dispersion and
advection. Basically, the flow pattern was first obtained by ignoring the
density differences. Then the temperature distribution was deduced by
using the known flow pattern. Thus the analysis can only be applied for
very small temperature differences. Also entrainment was ignored thus
the analysis becomes less accurate for large Froude numbers. Laboratory
experiments were also performed by Hayashi and Shuto. The temperature
determined experimentally was found to be consistantly lower than that
predicted indicating the effect of entrainment. Lean and Whillock (1965)
and Stefan and Schiebe (1968) performed experiments on the two-dimensional
surface jet problem. Their results are consistant with the findings in the
present investigation. However, insufficient details were reported to allow
detailed quantitative comparison. Jen, Wiegal aid Mobarek (1966) and Stefan
and Schiebe (1968) performed experiments on the three-dimensional surface
jet. Jen, et at. dealed primarily with the case when the source densimetric
Froude number is relatively large. They found that the jet excess tem-
perature first decreases due to jet mixing followed by a region where it
decreased at a faster rate. Stefan and Schiebe (1968) reported on similar
experiments for smaller values of the source densimetric Froude numbers.
Detailed measurements were reported. However no analysis of the data
was included.
54
-------
This chapter has been divided into several sections. Section 4.2 treats
the two-dimensional case in detail. In Sec. 4. 2. 1, the assumptions and
the resulting equations are derived. These equations turn out to be a set
of nonlinear differential equations with four parameters. Some general
properties of these equations are discussed in Sec. 4. 2. 2 where it will
be shown that the character of the solutions are strongly influenced by the
relative magnitudes of the parameters in the system. In particular, it
will be seen that for some combination of parameters, a continuous solution
does not exist. In Sec. 4. 2. 2, an approximate relation between the para-
meters will be derived which specifies the region in parameter space
where a continuous solution can be obtained. Section 4. 2. 3 examines the
solution to this system of equations for various values of the parameters.
In the event a continuous solution is not possible, it will be shown that an
internal hydraulic jump may be developed. The flow field before the jump,
the jump conditions, and the flow field after the jump will be derived and
discussed. The problem of matching the solutions at the jump is discussed
in Sec. 4. 2.4. It will be found that for certain combinations of parameters,
no jump can be found to match the solutions indicating that the source is
actually inundated in these cases. A nomograph will be presented in Sec.
4. 2.4 which divides the parameter space into three regions: a) the region
of jet-type solution, b) the region where the solution is characterized by
the presence of an internal hydraulic jump and c) the region where the source
is inundated. Section 4. 2. 5 summarizes the findings in the previous sections
and describes the procedure of finding the solution for given parameters
by using a computer program SBJ2 listed in Appendix B.
In Sec. 4. 3, the analogous axisymmetric problem is investigated. The
equations are derived and some simplified cases solved. The general
features of the solution also depends on the relative magnitudes of the
parameters. However, the division of parameter space is much more
involved and has not been examined in detail. The detailed study of the
axisymmetric case paralleling that done for the two-dimensional case
should be performed in the future.
55
-------
4. 2 Formulation and Solutions for the Two-dimensional Problem
Consider a two-dimensional surface source of buoyant fluid at the origin
as shown in Fig.4.l. Let the density of the ambient fluid, assumed infinite
in extent and motionless, be Also, let the source be characterized by
a discharge velocity U, source dimension h and discharge density p 1 <
• Assume that the source densimetric Froude number
U 2
F
0 / _____
g( ih
/ °
is sufficiently large so that near the source, it can be expected that the
phenomenon is similar to that of an ordinary two-dimensional jet. Thus,
except for a zone of flow establishment, one would expect that the velocity
distribution is very nearly Gaussian. Entrainment of the ambient fluid
would occur along the jet and the jet dimension would grow with distance x.
Let U(x), (x), h(x) be the mean velocity, density and thickness of the jet
at x. Laboratory experiments by Ellison and Turner (1959) indicate that,
unlike an ordinary nonbuoyant jet which expands linearly with x, the buoyancy
of the efflux tends to decrease the entrainment rate. In particular, it was
found that the entrainment coefficient e is a monotonic decreasing function
of the local Richardson number, Ri, defined as where
• g( 0 -p ) -
R i h/U . e - d
Co x
Their experimental finding is reproduced here as Fig.4.2, and can be well
approximated by the formula
r l.75
e = 0.075 ________ - 1 0 < Ri 0.85
Ri -
0. 85
0 otherwise
It is noted that entrainment ceases for Richardson number exceeding a certain
critical value Ri. Thus the buoyant jet does not expand linearly with x,
56
-------
y
h )
p 0
Figure 4.1
Definition sketch
Ri
Figure 4.2
Entrainment coefficient e as function of Richardson
number Ri.
x
e
0 0.2 0.4 0.6 0.8
57
-------
but rather takes a shape as indicated in Fig.4.l. At some distance from the
source, when the local Richardson number reaches the critical value, the
jet ceases to expand and the phenomenon resembles a two-layered stratified
flow. From that point on (and probably some time before), the flow can no
longer be classified as a free turbulent flow and the mechanism of inter-
facial shear should play a part in determining the flow pattern. It is seen
that for the maintanence of positive flow, the interface must possess a
slight positive slope in order to overcome the interfacial shear (see Fig.
4.3). This further suggests that at some x = X the interface may meet the
free surface, leading to the observation that no steady state solution may
exist for this problem.
We now note that as time goes on, X must increase to accommodate the
continuing efflux. Thus the maximum thickness of the buoyant fluid would
also increase. When this thickness exceeds that which can be provided
by the jet through entrainment, an internal hydraulic jump must occur.
As X increases further, the jump would occur sooner until a point is
reached when the source is inundated. From that point on, the source
momentum drops out of the picture entirely.
From the above discussion, it is seen that the phenomenon of a horizontal
buoyant jet discharged at the surface of a quiscent fluid of infinite extent
may possess features very similar to open channel flow. Near the source,
we may encounter jet type flow analogous to supercritical flow in open
channels while far away, the phenomenon is similar to subcritical flow in
open channels. The subcritical flow region can, therefore, be expected
to be influenced by downstream constraints. For example, in the event
a laboratory experiment is performed on this phenomenon, the conditions
at the downstream end of the tank or flume can strongly influence the
flow field.
In this investigation, we are concerned primarily with the case when the
buoyancy of the efflux is due to heat. In this case, there is the added
mechanism of surface heat exchange between the water and the atmosphere.
This mechanism now takes the role of imposing the downstream constraints.
58
-------
y
ENTRAINMENT
p 0
x
SHEAR
Figure 4. 3 Definition sketch.
-------
From the above discussions, it is seen that the flow field induced by a
surface horizontal warm jet can, in general, be divided into four zones
as shown schematically in Fig.4.4. Zone I is the zone of flow establish-
ment. Zone II is the supercritical region where the flow is basically a
jet with decreasing entrainment rate. Zone Ill is an internal hydraulic
jump while Zone IV is the subcritical region where interfacial shear and
surface heat exchange play dominant roles.
It should be remarked that not all these zones may be present in any given
situation. In particular, as will be shown later, Zones III and IV may be
absent if K, the surface heat exchange coefficient, is sufficiently large.
In that case, the flow field is similar to that in an ordinary jet. This is
reasonable physically since if K is very large, the buoyancy would be lost
to the atmosphere quickly and the resulting jet is virtually not buoyant.
On the other hand, if K is sufficiently small, it will be seen that the source
may be inundated and Zones I, II and Ill may be absent. Thus, given all
the other parameters, there exists two critical values of K; Kcr+ and
K with K > K , such that if K > K , Zones III and IV are absent
cr- cr+ cr- cr+
and if K < K , Zones I, II, and III are absent. For K between K and
K+. all the zones are present. The model to be developed in the following
sections will predict whether all the zones are present and also locate the
hydraulic jump when it occurs.
Since the mixing processes involved in the various zones are quite different,
the above discussion is of importance in design considerations. For ex-
ample, if it is desirable to achieve initial jet mixing so that the temperature
drops quickly, then the discharge structure must be designed so that the
source is not inundated. On the other hand, if maximum surface heat loss
is desired then it would be desirable to achieve inundation. However, the
depth of inundation should not be too large so as to cause the discharge to
interfere with the intake of the cooling water.
The formulation and the solutions to be presented in the following is pri-
marily concerned with Zones II and IV with a brief discussion of Zone III.
60
-------
U
The zones in a surface horizontal buoyant jet.
Figure 4.4
-------
The matching of the zones and the conditions under which some zones
are absent will also be discussed. However, no discussion will be given
on Zone I, the zone of flow establishment, since it is similar to the
corresponding zone in an ordinary jet which is well documented in the
available literature.
4. 2. 1 Derivation of Equations
In this section, the governing equations will be derived for the two-
dimensional surface horizontal buoyant jet. In the formulation presented
in this section, it will be assumed that
a) the velocity and density deficiency profiles in the vertical
direction are similar in shape
b) a steady state solution exists
c) Boussineq approximation: density differences are only
important in modifying the gravity term.
d) the flow is primarily horizontal (boundary layer assumption)
It should be pointed out that the similarity profiles to be used in Zones II
and IV may be different.
Under these assumptions, the equations of motion are as follows:
Continuity:
- + - = 0 (4.1)
Morn e n tu m:
+ - —(uv) = - - - - + (4.2)
p x
0 = -- - og (4.3)
62
-------
Cons erv3tion of density deficiency:
+ (4.4)
y y
where x, y are the horizontal and vertical coordinates
U, V are the velocity components in the horizontal and vertical
directions
p is the pressure
is the kinematic shear stress = shear stress/density
is the density
p is the density of the ambient fluid
D is the diffusivity
We shall now utilize the as surnpti on of similarity and further specify
that
u(x,y) U(x) f( T 1 )
p(x,y) - T(x) f ( - )
where y (x) is the free surface elevation and h(x) a characteristic
thickness of the spreading layer. The function f( ) specifying the shape
of the similarity profile will be left arbitrary. Examples may be f( )
(ij- 1
e (near the source) or f( ) (far from the source).
O I:!> 1
We now integrate the Eqs. (4. 1) through (4. 4) from y = - to y = Ti In-.
tegration of the continuity equation gives
a (Uh) U - v(Ti) +
wh e r e
r 0
a. J f( ) d.
Now the kinematic free surface boundary condition is (y- ) O on y
63
-------
i.e.
_U- - + v(y ) 0
Thu s,
(Uh) v(- )
dx -
It is generally assumed that v(- ) (the entrainment velocity), is related to
U, the characteristic velocity, by an entrainment coefficient e which can
be a function of the Richardson number (Ellison and Turner 1959)
e = e(Ri) , Ri = g Th
p 0 U
TI -rn s
= - U (4. 5)
dx a
In a similar fashion, the other Eos. (4.2), (4.3) and (4.4) can be integrated. The
only term requiring some explanation is the pressure integral. From
Eq. (4.3),
r
p(x,y) =-g c(x, ) d
TI
Letting c = p 0 - 9(x, y), and using the similarity profile
e(x, ) = T(x)
we obtain,
y-r 1
- = p(y - ) - T f( hr) = p(y- ) -Th
TI 0
from which
64
-------
d
1 - — (Tb) f( )d( - Tb f( ) d
g x - - 0 dx dx h dx h J
0
ai
As y - -ce, we expect — -, 0 since there is no motion horizontally.
Thus
a d (Tb)
dx p dx
0
We now calculate
y-r1
ri -i rY h+i Ll
(Th) 1 • f( )dC+Thf( h h dx ] ] dy
_j _ x
- - 0
- (Th)lJrd T$ 1 f()ci + ThIf
Y ) dy + Th j h dx h
h h
- - - -
After interchanging the order of integration for the double integral, and
carrying out one of the integrals, we get
‘ Th)’ (Tb ) 2
1 r - dy - -- (Tb 2 ) i (f(c)d + - a
gJ x dx 0
-
T
Note that — << 1 so that the second term on the right can be neglected
p
0
compared with the first.
Returning now to Eq. (4.2), and integrating with respect to y from - to Ti
and using the kinematic free surface boundary condition, we get
- (U 2 h) a 1 (Th 2 ) + (T 5 - T ) a 2
dx
where
a 1 = [ f( )d/Sf 2 ()d ; a 2 1
- f 2 ( d
•1
-
65
-------
and -. are the shear on the free surface and the interface respectively.
Equation (4. 4), when integrated with respect to y from - to ri gives
- —(ThtJ) = a 2 D H: a 2 H
where H is proportional to the rate of heat loss at the free surface.
To summarize, under the assumptions made, the equations governing
the surface spreading are
C on inn ity
—(Uh) = (4 6)
Iviornen turn:
(U h) = ai (Th 2 ) + a 2 s - (4. 7)
Heat balance:
- --(ThU) - H (4. 8)
The Eqs. (4. 6), (4.7) and (4. 8) constitutes a set of three ordinary
differential equations for the unknowns U, T, and h subject to the given
conditionsU U , T T , h=h atx=O. Theyarederivedfromthe
0 0 0
Navier-Stokes equations by making the assumptions stated in the beginning
of this section.
Before these equations can be solved, it is necessary to specify the functions
H, ‘ and . - as functions of the other unknowns. We shall neglect
the shear at the free surface.
66
-------
For ., we shall assume
1
U
T. =
1 h
where is an effective viscosity coefficient. In the analyses to be pre-
sented in the following sections, € is taken to be constant. However,
the general features of the solutions, as well as the method of solution,
are equally applicable when c is not constant but depends on, say, U.
The critical relations to be derived in the following sections will be different
when is not constant. However, the procedure for finding them will be
similar.
For the quantity H, we shall assume
H=-KT
where K is an effective heat exchange coefficient. In this investigation
K will also be assumed constant. Again the general features of the solution
and the method of solution are also applicable when K is not constant
but depends on, say, the temperature.
4. Z. 2 General Discussion of the Equations and Solutions
Before analyzing quantitatively the flow field in the several zones, it is
instructive to examine some of the properties of Eqs. (4. 6), (4.7) and
(4. 8).
It will be shown that for some combinations of input parameters, a
continuous solution cannot be found. In this section, the region in para-
meter space where a continuous solution can be obtained will be delineated.
In the event the parameters fall outside the region, then the solution is
either non-existent, or discontinuous or not source governed. These cases
will be discussed in detail in later sections.
67
-------
Normalizing the variables by the source conditions U, h, T,
U h T x
u-- — fl—, h - - = -h—, T = X =
0 0 0 0
-U 2 h U K
0 o o 2
and defining F =— - , R = , k =
iT U 0
equations (4. 6), (4.7) and (4.8) become, after dropping the s,
d (uh )
= eu (4. 6 a)
dx
d 2 1 d 2 lu
dx (u h) = - 2 F (Th ) - - (4.7a)
—(uhT) = - kT (4.8a)
dx
For application to practical situations, it is necessary to obtain the
numerical values of the ‘s based on assumed form of f (c). For ex-
ample, iff( ) fi - l
-------
Equations (4. 6a), (4. 7a) and (4. 8a) constitute three nonlinear ordinary
differential equations to be solved for the three unknowns u, h, and T
subject to the conditions u h = T 1 at x 0. The system of equations
depends also on three parameters F , R and k. It will be shown later
0
that the character of the solutions are strongly dependent on the relative
magnitudes of these three paramters. Thus it is important tofirst in—
quire into the typical orders of magnitude of these parameters.
Typical values of the parameters K and are very small. One expects
K to be of order l0 or l0 ft/sec while to be of order 10 to lO_2
ft 2 /sec. Thus for U , h of order unity, k 0 (l0 ), and R 0 (i0 ).
0 0
F 0 , however, can vary over a wide range. Since T is small, F can
be expected to be relatively large. We are, therefore, primarily interested
in the case k very small and R large.
It can be readily shown that Eqs. (4. 6a), (4. 7a) and (4. 8a) may be put into
the following form:
dT kT T
-a-- -----h -e - (4.6b)
e hT k hT 1 1
Ze+ZF Z ZF 3Rhu
dh - ou OU (4.7b)
dx 1 hT
0 U
du u dh
= [ e - (4. 8b)
It can be readily seen that the denominator in Eq. (4. 7b) is simply - 1
where F is the local Froude number which is the inverse of the Richardson
number. It now becomes obvious that for the existence of a continuous
solution, it is necessary that if F -, 1, the numerator in 4. 7b must also tend to
zero. It will be seen in the following discussion that this necessary condition
is not satisfied except possibly fortuitously. On the other hand, for some
combinations of the parameters F , k, R, the local Froude number never
0
approaches one so that a continuous solution does exist.
It is convenient to rewrite Eq. (4. 7b) in the form
e 1 k 11
dh = { - - - +
(4.9)
69
-------
We further note that
ThF = F u
0
Hence
dT dh dF 2 du
T+h F u
therefore,
dF 2 du dT dh
F u T h
Thus, F would increase or decrease according to whether
z j_ -
u T h
is positive or negative. From (4. 8a) we have
T’ h’ u’_ k
T + h + h - - uh
and from Eq. (4. 6a) we have
u’ - e h’
uh h
therefore,
F’ — 31 e h’ 1 k
F - Lh h + uh
Since h is always positive, F would increase or decrease according to
whether
3(e - h’) + -
is positive or negative. In particular, for F to increase, we need
70
-------
> 3(h! - e)
Using Eq. (4. 9), this condition becomes
k 6F 1
— > 3e +
u 2F+1 Rhu
Thus,
dF . 6F 1
> 0 if — > 3e + ZF+1 (4. lOa)
and dF . I c
= 0 if — = (4. lOb)
ux U
and dF . Ic
— < 0 if — (4. lOc)
dx u
With the help of Eqs. (4. 9) and (4. 10), we can now discuss the influences
of the parameters Ic, R, and F on the characteristics of the solution.
As was mentioned in the beginning of this section, we shall be concerned
only with the case when F is fairly large, while Ic and hR are both very
small, since this is the case of practical interest. Note that h, u, F,
and e are all positive quantities for a valid solution.
If k = hR = 0, then Eq. (4. hOc) is satisfied at x 0 and F would decrease
and asymptotically approach the critical value Fcr at which e 0 (see
Fig. 4.2.
If k = 0 and R 0, then condition (4. hOc) is always satisfied. F would
continue to decrease past F and reach F = 1 at which point -,
cr dx
leading to the non-existence of a continuous solution.
If = 0 but k 0, then we may have Eq. (4. lOc) satisfied at x 0.
Thus F decreases which implies e decreases. When it decreases
sufficiently so that k = 3eu, then F would increase again. Thus, F never
decreases below the value F since if it does, e would be zero and F
cr
would increase since (4. hOc) would be satisfied.
71
-------
Finally, we shall consider the most interesting case physically when both
and k are not zero. Suppose first that k is very large. In particular,
suppose k> 3 e + . Then (4. lOa) is always satisfied and we would have
F increasing all the time. When F becomes larger and larger, the solution
becomes more and more nearly that of an ordinary jet. On the other hand,
if k is very small, we would expect F to decrease initially. In that case,
u decreases, h increases and hu increases with x so that (4. lOa) becomes
more likely to be satisfied as x becomes larger. For sufficiently small
k, however, F would decreases to F = 1/Ri at which point e = 0 and
cr cr
hu becomes constant. The right hand side of condition (4. 10) then becomes,
6F
cr 1
ZF -i-i R(hu)
cr
It is clear that this may still be larger than - . The question then becomes
whether condition (4. lOc) is satisfied all the way to F 1. If so, we expect
no continuous solution. Since both k and are small in practice a good
estimate can be derived for the critical value of k such that below that,
no continuous solution exists as follows: since e = 0 for F < F , the
cr
critical value for k must be such that
2
k Rh 1
where h 1 is the value of h at the critical point. For k and small,
this value of h 1 can be assumed to be approximately equal to the asymptotic
value of h as x -. in the solution for k = - = 0. These solutions are
exhibited in Fig. 4. 5 in Sec. 4. 2. 3. When this is done, it is found that the
critical relation may be approximate by
-0. 655
(kR) 2.9 [ F] (4.11)
Thus for
(kR) < 2.9 [ F] 0 655 (4. lla)
we would encounter a discontinuous solution, while for
72
-------
(kR) > 2.9 [ F] 0 655 (4. lib)
we would have a jet type solution. However it should be pointed out that
even though we have called it a jet type solution, the flow field may appear
quite different from that in an ordinary jet. In fact, even if condition (4. llb)
is satisfied, F, the local Froude number may first decrease and then
increase. The ordinary jet is characterized by a local Froude number of
infinity. Thus the jet region may behave somewhat differently until F
becomes quite large.
From the above discussion, it is seen that the flow field in a horizontal
buoyant jet at the surface can be very different from that in either an
ordinary non-buoyant jet or a submerged buoyant jet. For the case when
no heat loss occurs at the water surface, no steady state solution is possible.
The source will, sooner or later, be inundated. For non-zero heat ex-
change, a steady state can be found as will be demonstrated in the following
sections. Moreover, given all the other parameters, a critical value of
k, the surface exchange coefficient exists such that for values of I c larger
than this,jet type solution may be found while for k less than this critical
value, the solution may be discontinuous. For Ic and very small, this
critical relation is given approximately by Eq. (4. 11). In practical situations,
k is expected to be very small and the condition of (4. ha) is likely to be
satisfied. In the following sections, this case will be considered in detail,
since it is the case of practical interest. It is also the case which results
in the most complicated flow pattern.
4. 2. 3 Solution of the Equations
In the beginning of Sec. 4. 2, it was deduced from physical reasoning that
the flow field induced by a surface horizontal warm jet can be divided into
four zones as shown in Fig.4.4. In the general discussion in Sec. 4. 2. 2,
it was found that a jet type solution may be expected if condition (4. hlb)
is satisfied. In that case, Zone II, the jet region extends all the way to
73
-------
infinity. If condition (4. ha) is satisfied, we expect Zone II to extend
at most up to some distance from the source. It will be shown later in
Sec. 4. 2. 4 that there exists another critical relation between the parameters
such that if satisfied, Zone II is absent entirely and Zone IV extends all
the way to the source, inundating it. In this section, we shall examine the
flow field in the Zones II, III and IV separately.
( A) Zone II
Since k and are usually very small, we shall first investigate the case
when both are zero. In that case, Eqs. (4. 6), (4. 7) and (4. 8) become
=
dx a.
- --(UhT) = 0
dx
- (U 2 h) =
We shall consider U, h, T as mean quantities over Lhe vertical and
specify f( ) = -1 < so that a = I and = - . It can be
readily seen that choice of f( j to be some other profile will only change
the coefficients a. and 1 slightly without affecting the essence of the
solution.
Let U , h , T be the source conditions at x = 0 and define dimensionless
o o 0
quantities u, h , T , x* as
U = U/U 0
h* = h/h
° (4. 12)
= T/T
= x/h
0
74
-------
Substituting into the equations and dropping s, the problem then becomes
(uh) = eu (4. 13)
- -(uhT) = 0 (4.14)
- --(T h 2 ) (4. 15)
U 2 p
F ° ° (4.16)
o gTh
00
with u h T 1 at x = 0 (4. 17)
From Ellison and Turner (.‘95 , (see Fig.4.2), the quantity e is a function
of the local Richardson nu . ber w ich is the reciprocal of the local Frou de
number
1 Th
R i = = - -
I .
U
C,
It can be approximated by the function
I
I ‘. -
e — ----- 1 , R
) cr
c(Ri) Rlcr (4. 18)
o , Ri >Pi
cr
where e is the value of e at Ri 0, Ri is the cri ica1 Richardson
o cr
number beyond which e 0, and n is an exponent. From the data, we
may deduce the following values for the parameters.
e 0.075
0
Ri 0. 85 (4. 19)
cr
n 7/4
75
-------
Equations (4. 13), (4. 14) and (4. 15), with e given by Eqs. (4. 18) and (4. 19),
subject to the condition of Eq. (4. 17) constitutes an initial value problem
with only one parameter F. The solutions have been affected by using a
fourth order Runge—Kutta numerical scheme and are shown in Fig. 4. Sa.
for a variety of F’s.
In the case k and are not zero, the Eqs. (4.6), (4. 7) and (4. 8) can also
be solved if condition (4. lib) is satisfied giving rise to a jet type solution.
A computer program SBJ2 written in Fortran IV language is included in
Appendix B with which this solution may be obtained. A few examples are
shown in Fig. 4. Sb. In the event k and are not zero but condition (4. 1 ib)
is not satisfied then it can be expected that an internal hydraulic jump
would occur either away from or inundating the source. The method of
finding the solution is to use the program SBJZ twice, first to solve the
jet portion to the point of the jump and then to continue the solution by
re—initializing the program SBJ2 with the parameters just after the jump.
The method of matching the solutions and locating the jump will be discussed
in Sec. 4. 2. 5 after we have investigated the flow fields in Zones III and IV.
The case when the source is inundated will also be discussed in Sec. 4. 2. 5.
It should be remarked that if k and are such that condition (4. lib) is
not satisfied, the solutions obtained herein for k = 0 may be used to match
with the flow in Zone IV to within a good approximation.
( B) Zone III
From the general discussion, it is seen that an internal hydraulic jump
is a possibility in the development of a surface layer. The conditions
across the jump will now be obtained. Consider an abrupt internal
hydra Lc jurno .vii upstream conJ ons u h 1 and do\vnstrearn
conditions 02. u 2 h 2 as Shown in Fi . 4. 6.
Conservation of mass requires
-f = o 2 u 2 h 2 (4.20)
76
-------
x l NORMALIZED DISTANCE
Figure 4.5a
Predicted jet thickness and density deficiency in
a surface horizontal buoyant jet for k = I/B = 0
(two dimensional case).
0 20 40 60 80 100 120
V I)
U)
w
z
0
I
I-.
0
I ii
N
-J
4
140
0
2
4
6
8
I0
12
14
16
• 0075
I/RO
k:O
-------
C-)
z
U i
C -)
Li
Ui
U)
z
Ui
O
Ui
N
-j
4
z
F-
1.0
0.8
0.6
0.4
0.2
0
Co
N
0 20 40 60 80 100 120 140
x, NORMALIZED
DISTANCE
Figure 4.5a Gontinued.
-------
, NORMALIZED DISTANCE
Figure 4.5b
Predicted jet thickness in a surface horizontal
buoyant jet for the case when k k cr+ (two dimen-
I III I I III I I II I-
I I III
Cl)
Cl)
w
z
C-,
I
0
F 0 — C X )
F 0 : 12
R:U
100
I0
F 0 : II, k:O.I
R :IOO
F 0 :IO, k :O.OI
RIOO
0.1 I 10 100 1000
sional case)
-------
h 2 , u 2 , p 2
Figure 4.é
Definition sketch.
p 0
-------
where E the voiume of fluid entrained in the jump.
Conservation of momentum requires
- 01 u 1 2 h 1 + p 2 u h 2 r( 1 - p 2 ) dy
where the integral extends over the vertical sides of the control surface
dotted in Fig. 4. 6. We note that
02 h 2 p 1 h 1 + p(/\h) 2
where ( h) 2 is the jump of the interface.
The integral r( 1 - p 2 ) dy is thus
2
c gh 2 o 1 gh 1 2 o 1 gh 1 + c 2 gh 2 ( h)
2 - L 2 - 2
Since
02 h 2 - p 1 h 1
( h) 2 =
00
the momentum conservation equation gives
p 2 u 2 2 h 2 - p 1 u 1 2 h 1 [ p 1 h 1 2 _o h + g 2 2 2 2
2 2 200 L°2 h 2 h 1 (4.21)
In addition to the conservation equations for mass and momentum, there
is also the conservation equation for buoyancy. Thus
(p - p ) u h - p 2 ) u 2 h 2 (4. 22)
1 11 0
The three conservation Eqs. (4. 20), (4. 21) and (4. 22) now allow the jump
conditions to he determined. For our purposes here, we will assume
E = 0, i. e. , no significant entrainment occurs in the jump. Invoking the
Boussi.nesq approximation, we get
u 1 h 1 = u 2 h 2
81
-------
and
p
The momentum equation, (4. 21) then gives
2
II = ± [ i + 8F - 11 where F 1 23)
U 2 h 1 2 1 j 1 p-c
g( ° )h
Thus the velocity ratio and the thickness ratio are expressed in terms of
the upstream densimetric Froude number. We note that if F 1 < 1,
h h
1 while for F 1 > 1, > 1. To find if either or both are admissible,
w note that the total head H ‘upstream of the jump is
2
H
while that downstream is
2
U 2
H 2 + 2 (h 2 -
where zero bead has 1 en referred to the free surface upstream of the
jump. The difference ‘ : ad = H - H 2 is therefore
= (l+ )(
Thus, is positive if h?/hj > 1 and negative if h 2 /h 1 < 1. But A being
negative implies a gain in total head in going across the jump which is
clearly impossible. Hence a jump can only occur if F 1 > 1. The jump
considered herein is very simple. A more detailed theory of internal
hydraulic jumps in discrete layered fluids can be found in Yih (1965).
82
-------
( C) Zone 1V
Having found the solutions for Zones II and III in the previous subsections,
we must now consider Zone IV. It should be remarked here that the solution
obtained in Zone 1V would determine where the jump (Zone III) would occur.
In fact, it is possible that Zone II and Ill are completely absent if the source
becomes inundated.
In Zone IV, we have F < 1. Thus the entrainment is zero. The Eqs.
(4. 6), (4.7) and (4. 8) then become
— -(Uh) 0
dx
(UhT) - KT
dx
d (U 2 h) - d Th 2 ) -
dx
0
Normalizing U, h, T with respect to the conditions at the beginning of
Zone IV, U ,h , T and dr opping s as before, we get
0 0 0
(uh) 0 (4. 24)
dx
— --(Thu) - kT (4.25)
dx
d 2 ., 2 Ui
— (u h) = - - -- — Th ) - (4. 26)
dx iLi c lx
0
K
where x 0 corresponds to the beginning of Zone IV, k = -fl—, and
Uh o
R= oo.
€
Although Eqs. (4. 24) and (4.25) are easily integrated, Eq. (4.26) does not allow
closed form soIn cns. To gain some physical insight into the situation,
we note that since in Zone 1Y, we have F < 1, the n-iomenturn flux is
0
very si-nail. Thus we shall first examine analytically the solutions by
neglecting the inertia term. In that: case Eq. (4. 26) becomes
83
-------
I c i — U
2F f LIl J — _i-
0
The Eqs. (4. 23) and (4. 24) readily integrate to
ub = I
- kx
T=e
so that we have
___ - (4.27)
where y R is a ratio of the source Froude number to source Reynolds
number. It can be readily shown that the solution to Eq. (4. 27) is
h 4 e 2 - - + - e’ (4. 28)
Note that if > 1, then for sufficiently large x, becomes negative.
This implies there is no steady state solution. Since
4F
2y o
kkR
this condition is the same as
4 4cq 0 2 c 0
h 0 < gTK (4.29)
If condition (4.29) obtains, no steady state solution exists. This means that
the internal hydraulic jump must occur so that the parameters following
the jump are such that
2
4eq 2 0
h 2 g T 2 K (4. 30)
where the subscript 2 in Eq. (4.30 )has been inserted to stress the fact that
these are the downstream conditions after the jump (and the initial conditions
for Zone lv). If the solution in Zone II (jet region) is such that condition
(4. 30) cannot be met at all by an internal hydraulic jump, then the source
will be inundated to satisfy (4. 30).
84
-------
The considerations given on the previous page lead to the condition given
by Eq. (4.30) which involves the assumption that the inertia or momentum
flux of the flow is negligible compared with the pressure induced forces
and the viscous forces. This assumption may not be adequate under all
conditions. Return now to Eqs. (4. 24) through (4. 26). We first integrate
Eqs. (4.24) and (4.25) to get
uh 1
-kx
T e
Substituting these into Eq. (4. 26), letting kx, s , and rearranging
we get
1 4-
2F h e - s
0 (4.31)
--1
0
This is a first order equation with two parameters F and . Moreover,
since we are discussing Zone IV, F 0 < 1. It should be realized that for a
physically realizable solution, h must be bounded for all finite . We now
note that if 2F s, then Eq. (4.31 ) does not possess any physically
0 dh . dh.
realizable solution since - -- at = 0 is not positive. If — - is negative
1 dh
at some , say at = 0, (i. e. , 2F < s) then will stay negative until
the denominator changes sign since ftie numerator will never change sign
dh . .
as long as < 0. But the denominator changing sign implies it becomes
d . dh dh
zero for some which leads to an infinite . Thus > 0 everywhere
for proper solution. For proper solution, it is necessary (though not
sufficient) that F < . The larger the s, the smaller can be the value
F . Given the numerical value of s, there is then a critical value of F
0 0
say F such that proper solution can exist only if F < F . This
ocr o ocr
value, F corresponds to the case when the energy lost in the internal
ocr
hydraulic jump is minimum. We now proceed to find the critical relation
between F and s.
0
For purposes of discussion, it will be convenient to use the local
densimetric Froude number
85
-------
2 F
U 0
F = F — = — e
oTh h 3
Equation (4. 31) can then be written in its alternate form
h
dh - - S
— 1 (4.32)
F
We note that for proper solution - > 1 and > 0 always. Hence
> s always. Since h T > 0, h is monotically increasing with
However, - may decrease to nearly 1 and then increase again. If this
occurs would become large unless the numerator also approaches
zero. Thus we expect that as F 4 1, -. s. In particular, if F is
nearly 1 we must have s nearly . Thus we have found that one point on
the critical relation
F = F (s)
o ocr
is
1 F ( k)
o cr
To obtain the critical Froude number for other F < 1, we note that given
that we have the correct F 0 for the given s, then integration of Eq. (4. 31)
would continue with F increasing such that at = F becomes almost
1 and from on, F decreases again. We note that we can solve the problem
from backwards in and thus obtain F (s). Let h F
1 dF ocr 1 1
be the values at when - - - = 0 and F 1 = 1 - 5 where S > 0. We
now define new variables h , ri by
= I = +
Eq. (4. 31) then becomes
5 1 4 r
—- h e
dh 1 1
dr 1 3 r
—h e- 1
F 1
86
-------
Moreover, from previous considerations. This equation subject
to = , F 1 = 1 - 6 has been solved and from the solution the critical
re1a ion F cr (s is deduced and plotted in Fig. 4. 7.
It must be pointed out that F is never equal to 1 anywhere in the flow field.
The critical Froude number F cr is just an upper bound for the existance
of a proper solution. In a practical situation, it can be expected that the
internal hydraulic jump would occur in such a location that the critical
relation is nearly satisfied, i. e. , F F - 6 where 6 is a very small
o ocr
quantity.
It is interesting to note that the critical relation derived herein is very
nearly the same as condition (4.30) for F smaller than about 0.04.
The critical relation can be checked by attempting to solve the problem
with several values of F in the neighborhood of F . This was done
o ocr
for s = 1, 2, 5 and 10. The results confirm the critical relation obtained.
Thus, we have found that the solution in Zone IV (after an internal hydraulic
jump) is characterized by a) h is always increasing, b) F, the local
Froude number is always less than 1. Moreover, there exists a critical
relation between the Froude number F 0 at the beginning of Zone IV to
the quantity S = kh where h is the thickness of the buoyant layer at the
beginning of Zone 1 . This critical relation is shown in Fig. 4. 7. For a
proper solution, the value of F must be smaller than F cr
4. 2.4 Matching of Solutions
The solutions obtained in the previous section for Zones II, and IV must
be matched by the conditions in Zone III to give the overall quantitative
description of the phenomenon. We note that as the solution in Zone II
proceeds, we have the following quantities at each step of integration:
h, the jet thickness, T, the density difference, u, the velocity from which
we can obtain the local Froude number F. Let the subscript 1 be applied
87
-------
I
II I I iI
Focr (
F:J-
4s
(Eq 4.30)
Li
z
Li
0
U-
C-)
I-
Li
U)
z
Li
0
L i . 0
0.I
0.01
0.001
-I
0.1 1 10
: I/Rk
102 10
Figure 4.7 Critical relation F F (s).
0 ocr
-------
to all these quantities to denote the fact that they are upstream of the in-
ternal hydraulic jump. If the jump occurred at station where the
solution variables are h 1 , u 1 , T 1 and F 1 then from Eq. (4. 23), Sec.
4.2.3, the variables just downstream of the jump would assume the values
- h 1 r
h 2 L + 8F 1 - 1
2 ul
u 2 =
Jl+8F 1 -l
8F 1
F 2 - r
/l -
These quantities (immediately after the jump) must satisfy the critical
relation (Fig. 4. 7) given by
F 2 F (s)
o cr
___ U 2
where s andF =F
= ‘2 2 o T 2 h 2
Note that, in general, it is expected that as x increases, h 2 decreases
while F 2 increases. Thus s increases while F decreases. The
o cr
condition F 2 < F , if satisfied at x = x 1 , would therefore also be
ocr
satisfied for x x 1 . Thus the analyses given in the previous sections do
not give a unique location for the internal hydraulic jump. In order to
specify the location of the jump, it may be postulated that it will occur at
the farthest possible location from the source. This would correspond
to the smallest possible jump. However, this is only an assumption and
must be verified by experiments before it can be applied to a practical
problem.
It should also be pointed out that there may be no location where the con-
dition F 2 < F is satisfied. In such a case, the implication is that
ocr
89
-------
F 2 is too large or F 0 cr too small (or s = too large). The physical
interpretation of this situation is that the jump should have occurred even
before the source. In other words, the shear is too large and k too small
for the available source momentum to push the jump away. One expects
in this case that the source will be inundated.
The critical condition for which this occurs can be deduced if we note that
this critical state corresponds to the case when the jump occurs just at
the source. F 2 at the source can be obtained by the relation
8F
1:- 0
-
\l+8F -l
L 0
Thu s
F., = F ( s)
ocr
where
S =
and
h h [ [ l+8F - 11]
2 o o
Thus given F 0 , these relations together with the relation in Fig. 4. 7
allows the determination of the critical condition for inundation. This
has been done and is shown in Fig. 4. 8.
At the other extreme, it is possible that the condition F < F is
2 ocr
always satisfied but F = F is not satisfied. This would be the case
2 ocr
when k is large compared with represented approximately by condition
4. lib). As was discussed in Sec. 4.2.2, in this case the local Froude
number would actually increase with x. Thus F 2 would be a decreasing
function of x and h 2 would be an increasing function of x. This implies
s is a decreasing function of x and hence F 0 cr would increase with x.
90
-------
F 0 , DENSIMETRIC FROUDE NUMBER
Division of parameter space into regions of
different flow pattern
1
Jet Type
kR
I0
1
0.I
0.0I
Internal
hydroulic
jump
Inundation
I
I0
102
Figure 4.8
-------
Thus if F 2 < F cr initially, it would remain so all the time for this case.
Physically this case corresponds to one when an internal hydraulic jump
is unnecessary. The flow is internally supercritical all the time and when
F becomes large, the flow field resembles an ordinary submerged jet.
However, it is envisioned that this situation would not be likely to obtain
in any physical case since k is usually very small.
Condition (4. 11) is also shown in Fig. 4.8. With this figure, it is now
possible to determine, before hand, from just the source conditions and
the environmental conditions, whether the solution is of jet type, includes
an internal hydraulic jump or the source is inundated.
4. 2. 5 Summary and Discussion
In the previous sections, the dispersion of heat resulting from the hori-
zontal discharge of a two-dimensional warm jet at the surface into a
quiescent cooler ambient is investigated. The effects of source momentum,
source buoyancy, entrainment, surface heat exchange, and interfacial shear are
all included. It is found that, unlike the case of submerged buoyant jets or sur-
face nonbuoyant jets, the source characteristics are not the only parameters
governing the flow. Downstream conditions can play an important role
in influencing the entire flow field possibly all the way to the source,
inundating the orifice. In this investigation, the case of an infinite
ambient fluid is examined and it is found that the surface heat exchange
mechanism can replace the necessary downstream conditions. In particular,
it is found that the relative magnitudes of the source Froude number F,
source Reynolds number R, and the dimensionless heat exchange coefficient
k play an important role not only in the detailed quantitative description
of the flow field but also in determining the type of flow field. For example,
referring to Fig. 4. 8, it is found that for given F, if kR is larger than
the critical value given approximately by Eq. (4. 11) (topline in Fig. 4.8),
then the solution is of jet type. On the other hand, if kR is smaller than
the critical value given by the lower line in Fig. 4. 8, then the flow field
is not like a jet at all. In fact, the source is inundated. For kR between these
two critical values, then the flow field consists of a jet type region near the
source and a two-layered stratified flow region farther from the source with an
internal hydraulic jump between these two regions.
92
-------
The analysis presented herein allows the determination of the flow field
given the source characteristics and the ambient heat exchange. A com-
puter program is given in Appendix B to solve the problem numerically.
First, the values k, R and F are determined from the given source
0
conditions and heat exchange coefficient. Figure 4. 8 should then be used
to determine the type of solution to be expected. If it is found that kR is
larger than the upper critical value, the program will give the solution.
If it is found that an internal hydraulic jump should occur, then it is
necessary to first run the program until the local Froude number be-
comes nearly unity. Then the output of the program is examined to deter-
mine the location where the condition
F = F (s)
2 ocr
is just satisfied. This is then the location of the jump. The subsequent
flow field may now be obtained by a second run of the program with F =F
02
k k = k h 2 and R = u 1 h 1 R where u 1 , h 1 , h 2 are the values of u,
u 2 u 1 h 1 ’
h and h 2 at the location of the jump according to the first solution. Finally,
if it is found that kR is less than the lower critical value, then the source
will be inundated. Although the flow field in this case cannot strictly be
analyzed using the present technique especially near the source, still,
it is possible to obtain an approximate solution using this method by
requiring that inundation occurs to the extent such that h 2 , the dimensional
depth of inundation satisfies the condition
2
F 2 =
g —i h 3
2
Kh 2
where q is the unit discharge. With this new value of F = F , k =
0 2 q
R = , the program will give the solution.
It must be remarked here that the investigation described in the previous
sections are based on several assumptions. For example, the shear law
93
-------
adopted here is given by the relation
U
¶ =
where € is a constant, and cannot be justified rigorously. It seems
reasonable to assume that the shear should be an increasing function of
u and a decreasing function of h. The relation chosen clearly satisfies
these requirements.
The region where ‘r may be of import is in Zone IV where it may in-
fluence the solution in such a way as to lead to either an internal hydraulic
jump or inundation of the source. In Zone IV there is no entrainment so
that uh constant and thus the shear law chosen is equivalent to
2
= constant u
This is equivalent to saying that the skin friction coefficient is a constant.
The skin friction coefficient in pipes and channels are generally found to
be a function of Reynolds number and surface roughness. Since there is
no entrainment, the Reynolds number in Zone IV is constant. Thus the
shear law adopted seems to be a reasonable choice.
The numerical value of to be chosen in a specific case is difficult to
assess since there is little data on the subject, although there are some
(e. g. Lofquist (1960)). On the other hand, there is an abundance of data
on the shear coefficient in pipes and channels (e. g. , see Schlicting 1960,
Chapter 20). In addition, it has been found by Keulegan (1944) that in
the laminar case, the interfacial shear between two fluids for the case
when the upper and lower fluids are both of infinite extent can be given by
¶ = 0.196
where u is the relative freestream velocity between the two fluids, p
is the density, v the kinematic viscosity and x the distance downstream.
94
-------
The corresponding shear law for the laminar boundary layer on a flat
plate according to the well known solution by Blasius is
T 0.332 pu 2 (?)
Thus the shear at an interface is approximately one-half that at a solid
surface. Until adequate data on interfacial shear are available, it may be
proposed that the shear coefficient be taken to be half the corresponding
value for the case of a fluid-solid surface. If that is done, it is found
-3 -22
that the shear coefficient is of order 10 or 10 ft /sec.
The surface heat exchange law chosen in this investigation is
H = - KT
The rate of heat transfer is taken to be proportional to the temperature
excess above the ambient which is assumed to be at the equilibrium
temperature. This is an often used approximation to a very complicated
phenomenon. Typical values of the measured rate of heat exchange, for
example, is Lake Hefner and Lake Colorado City indicate that K is not a
constant but depends on the equilibrium temperature and wind speed.
This is certainly not surprising. Typical values of K are of the order
1O ft/sec or l0 ft/sec (Edington and Geyer, 1955).
In order to obtain the quantitative description of the flow field, it is
necessary not only to have the source characteristics but also the inter-
facial shear coefficient e, the heat exchange coefficient K and the entrain-
ment coefficient e as a function of the local Richardson number. The
formulation in this chapter assumes that the coefficients €, K are con stants
and the coefficient e(Ri) is as given by the experimental findings of Ellison
and Turner (1959). As data becomes more plentiful it may be found
that and K are not constants and e(Ri) is not as given by Ellison and
Turner. For example, c may be a function of u, h while K may be a
95
-------
function of T. In that case, the formulation can be readily modified to
incorporate them. It is believed however, that the basic features of the
findings are valid and once adequate data is established to fully define the
coefficients, this model will provide a good quantitative description of
the flow field.
96
-------
4. 3 Axisymmetric Surface Buoyant Jet
In this section, we shall investigate the axisymmetric analog of the surface
buoyant jet. A schematic diagram of this phenomenon is shown in Fig. 4.9
It can be expected that since the only difference is one of geometry, the
general physical nature of the phenomenon is the same as the two-dimensional
case. Thus for K, the surface heat exchange coefficient large, we expect
a jet type solution. For smaller K, there would be an internal hydraulic
jump. For K sufficiently small, the source may be inundated. For K = 0,
no steady state solution may exist.
The analysis to be presented in the following has not been carried out
to the same detail as in the two-dimensional case. The critical relations
between the parameters which divide the parameter space have not been
derived for this case. These critical relations are more involved because
there are now four parameters F, R, E and k to be defined later.
Moreover, they are all independent whereas in the two-dimensional case
E was absent and the parameters k and R were found to occur approximately
as a group kR for small k and l/R. Thus the parameter space for the
axisymmetric case is basically four-dimensional, and the critical relations
are two three-dimensional surfaces in the parameter space. These
critical relations, however, can be found and should be done in the future.
In the following, the governing equations will be derived and solutions will
be found for some special cases.
We shall make the same basic assumptions as in the two-dimensional case.
The equations of motion are then as follows:
Continuity:
+ + - = 0 (4.33)
r
97
-------
z
Figure 4.9
Definition sketch.
r
00
98
-------
Mo mentum:
u - + w- -- = -- -- - + - (4.34)
p r
= - pg (4.35)
u 2 + w E = (4.36)
2
The assumption of similarity implies
u(r,z) U(r) f( Z ) (4.37)
p - p = T(r) f( fl ) (4.38)
where fl(r) is the free surface profile and h(r) is the characteristic depth
of the spreading layer. Integration of Eq. (4. 35) with respect to z from -
to using (4. 38) gives, as in the two-dimensional case,
rh
p - P g(z-fl) + g Th f( ) d
0
so that
- —i-- = - g - j- - -s-- T h f ( fl) - _ ( fl) (4.39)
We now integrate the Eqs. (4. 33), (4. 34) and (4. 36) with respect to z
from - to . Equation (4. 33) gives
(Urh) (r°f(C)dC)+ w( ) - = w
The Sf boundary condition is (z- ) = 0, which implies -U + w = 0
on z = . Therefore,
! (Uhr) = w where
99
-------
Introducing the entrainment coefficient e, we have
= (4.40)
rdr
Integration of Eq. (4. 34) yields
2
2 d z = — -- -
)+- ---(wu) +
L r r p r dz + - T.
- -
which gives
± (U 2 hr) i°f2( )dC= (Th 2 ) 0 0 _2
dr j f(C)dC + g
p 0 zp
- 0 -
+ TS - T.
Letting
PU
l/rf 2 ( )dC = a 2 , - Cf( )dç I j f (C)dC a 1 ,
0 .J
0
- - -
we have
d 2
r dr = a -a—- (Th ) + a 2 (TS - T ) (4. 41)
Equation (4. 36) can be written
- -(8w) +
+ r
where 9 = p 0 - p
Inte grating with respect to z from - to r gives
d r’ 1 dzr
rJ 8udz =
- -
100
-------
Now at r is the surface loss which we will assume to be, as
before, -KT. Hence
- — (TUhr) -Ka 2 T (4.42)
Equations (4. 40), (4. 41) and (4. 42) are the three equations for U, T and
h and are entirely analogous to Eqs. (4. 6), (4. 7) and (4. 8) of the last
section for the two-dimensional case. It can be expected that the character
of the solutions are similar. Just as in the two-dimensional case, the
relative magnitudes of the parameters would influence the character of
the solutions. Rather than discussing the general problem, we shall
first examine some special cases. As in the two-dimensional case, we
shall assume = 0.
Case 1) T = 0
For the case when the surface jet consists of the same fluid as the ambient,
then T = 0 and Eq. (4. 42) is absent. The entrainment coefficient may
be taken constant and i. is negligible. Then the equations are
d e
—(tjhr) = —Ur
dr
- -- (U 2 hr) = 0
It can be readily shown that
e
h —r
a
and 2
u h
to o 1
a —
\ e Jr
Thus the jet boundary grows linearly with r while the surface velocity
decreases as hr.
It should be remarked that if K is very large, then this should represent
an approximate solution to the problem for large r.
101
-------
Case 2) Surface Plume; e = K = 0
We next consider the case when the initial momentum at the source is
very small. This can only be the case if the velocity is very small.
Let U be the efflux velocity, h be the thickness at the source radius
0 0
r, then we are assuming
2rrr h U
000 0
rhU 2 0
Under the assumption that e = 0 and K = 0, the equations become
d
—(Uhr) 0
dr
- -(TUhr) = 0
dr
- --(Th 2 )
dr a 1
Just as in the two-dimensional case, we assume that Ti = where e
is an effective viscosity coefficient. We further take
o - < -1
f( ) { _ 0
so that a. 1 = - . The equations then reduce to
- --(Uhr) - --(UhrT) = 0
d(Th 2 ) - _2 p U
dr - g
With the conditions U = U, h = h, T = T at r = r, we get
Uhr = Uh r , T=T
000 0
102
-------
and
Z p Uhr
ZhT 000
odr - - g h 2 r
which integrates to
pUhr
h = h 4 - 0000
gT r
0 0
4
h gT
we note that at r = r = r exp 1 0 0 , h 0. Thus, this solution
max o i. p U h r J
is at best a quasi- steady solution. 0 0 0 0 If we accept this quasi-
steady solution as a time dependent solution, we find that the amount of
fluid contained in the spreading layer is
r r
max max h
V hZirr dr h — Z r dr
o h
0
0 o
But
h ( r r
— £n
h r I r
o max
so that
1
r 4pl r
V Zith r (Ln ° Ln (—) xdx
0 max. r / J L xJ
max r
0
r
max
For r >> r
max o
V hr (Ln
omax \ r
max
We observe that
p Uhr r
0 000 max
h =
o L gT
1 03
-------
As the pool spreads or as r increases the thickness h increases.
max 0
As r - so does h . However we note that h increases extremely
max o o
slowly with r . For example, let
max
r 2. 3 €p U h r
max r o oool 4
h A(log ) where A =
o r 0 gT J
then for
r =10r,h =A
max 0 0
while for
16
r = 10 r , h = ZA
max o o
16 12 .
If r is but one loot, 10 r is 2 x 10 miles which is eighty million
times the circumference of the earth Thus for all practical purposes,
the thickness can be assumed constant at say 1. 5A. In practice, the
uncertainty in the value of the coefficient € and the many assumptions
involved in the formulation certainty justifies this replacement.
Case 3) Surface Plume, e = 0,K 0
We next examine the case where K 0. The governing equations are
= 0
dr
- -(UhrT) = - KTr
dr
2 €p
- °
dr - - g h
We again let U = U , h = h , T = T at r = r . Then the first two
0 0 0 0
equations give
104
-------
Uhr =TJhr Q
000 0
K(r 2 -1 )
T T e
0
K/2Q
Substituting into the third equation and assuming e ° 1, we get
K 2
d . -- - -r 2epQ 1
drt ) gT 2
o hr
Let
Zep Q
K — 00
- ZQ’ “ gT
2
then _ _r -i cr - _L
drie 2
hr
with solution
2
4 2kr 2 1 2k rr -kr 2 dr
h =e h e -2y e — .
0 rj
r
0
We note that the function 1(r) is a monotonically increasing
function of r. Thus for given h , r , and k, there may be a value of r
o 02
- Zkr
such that I is larger than h 4 e ° /2y. In that case h 4 becomes
negative and the solution needs interpretation; one expects inundation to a
new value of h
0
We now note that by letting x = (rlr) 2 , then
2
i = Se = [ E 1 (kr 2 ) - E 1 (kr 2 )]
1 05
-------
As r- ,
I( ) = - E 1 (kr 2 )
where E 1 is the exponential integral which is tabulated.
-Zkr 2
Thus for legitimate solution, we need h 4 e ° yE 1 (kr 2), i.e.
2
Zkr
h e E 1 (kr 2 )]
Zkr
Thus for h > [ ye ° E 1 (kr 2)1 , we have
0 L 0
2 r
-Zkr 2
4 4 o r -kr dr
h = h e -Zy I e —
o r
r
0
2kr
For h < [ Ye 0 Ei(kr 2 )1 , we expect the source to be inundated to
Zkr
the level [ ye ° E (kr 2)1
L 1 0 J
Case 4) Discussion of the General Case
From these special cases, it is seen that the features of the phenomenon
is analogous to the two-dimensional case. Rather than considering further
special cases, we shall discuss the general case. The insight gained in
examining the two-dimensional case and in the special cases treated in
this section will be utilized.
We recall the governing equations
± - (Uhr) = (4.40)
rdr a.
1 06
-------
- -(U 2 hr) a 1 —(Th 2 ) - a 2 (T.) (4.41)
- - — (TlJhr) = - KI 2 T (4. 42)
Let the conditions at the source be: at r = r , U U , h = h , T = T
0 0 0 0
We now normalize the variables by these characteristic values. Thus
define u = U/U; T = T/T; h = h/h; r r/r. We get, dropping
id
— —(uhr) = Eu (4.43)
r dr
- -(u2hr ) - - - - —(Th 2 ) - (4.44)
ld
— - —(Tuhr) = - kT (4.45)
where
er -U 2 Uh Kar
E-—--- - F - ° R- °° k- 2o
- ha’ o - aTh ‘ - ar ‘ Uh
0 100 Zo 00
Thus the system depends on four parameters E, F, R, and k. The
parameter space is therefore four-dimensional. The critical relations
delineating the type of flow field are therefore three-dimensional surfaces
in the four-dimensional space.
The critical relations would divide the parameter space into three regions
such that given all other quantities, there are two critical values of k,
k and k . For k > k , a continuous solution may be expected.
cr+ cr- cr+
This solution would become, for large r, nearly the same as the one discussed
in Case 1. For k < k , one expects the source to be inundated and the
cr-
solution should resemble that discussed in Case 3. For k between these
critical values, a circular hydraulic jump would be encountered. The de-
tailed solutions would be more complicated than the two-dimensional case
1 07
-------
because the location of the jump is not only influenced by the four parameters
but also by the radial coordinate r.
It should be pointed out that although the assumptions made in deriving
these equations in the axisymmetric case, are basically the same as in
the two-dimensional case, there is a lack of experimental data on the
entrainment and shear coefficients. In the two-dimensional case, there
has at least been some experiments. Thus any results obtained herein
must be viewed with great caution. The system of equations 4. 43 through
4. 45 can be solved numerically for k> kr+ Some example solutions
are presented in Figures 4. 10 and 4. 11.
108
-------
U)
(I ,
LU
z
C-)
x
I.-
LU
— N
0 —
0
0
z
.
Figure 4. 10
r, NORMALIZED RADIAL COORDINATE
Predicted jet thickness in a surface horizontal
buoyant jet (axisyrnmetric case).
102
I0
E I .0
hR = 0.0
F = 100.
I0
-------
r, NORMALIZED RADIAL COORDINATE
Figure 4. 11 Predicted jet density deficiency in a surface
horizontal buoyant jet (axisymmetric case).
10_ I
I0 2
>-
C-)
z
w
0
U-
w
0
I-
01)
z
L i i
0
0
L i i
N
-J
0
z
I0
110
-------
4. 4 Example Applications
We conclude this chapter by presenting some examples on the two-
dimensional problem.
Example 1
Given h
0
U
0
= ift.
= 10_i ft/sec.
i0 ft/sec
= 10 ft /sec
K = 10 ft/sec
U 2
0
F 10,
0 gTh
00
k = = 10 ,
hR Uh = io 2
Thus kR = 1. Examination of Fig. 4. 8 shows that we are
of parameter space where the solution is of jet type. Thus
SBJ2 would simply give the solution with input F = 10, k =
e 0 = 0. 075. This particular solution was obtained and the
is shown plotted in Fig. 4. Sb.
Example 2
in the region
the program
lO , 1/P = 10 ,
quantity h/h 0
loft.
= 0.1 ft/sec
T
o ft/sec
g- - =
0
2
hence,
Given h
0
U
0
2
= 10 it /sec
K = 10 ft/sec
111
-------
then,
k = = = u h 0 = 10; F 0 = g h 10
Thus kR = l0 * From Fig. 4. 8, we expect inundation. The inundation
would occur until a depth h 2 , such that the following condition is just
satisfied.
(Uh) 2 < F
T ocr Kh
0 3 2
(gj—) h 2
Thus
1 ________
=
-------
/ -3
and hR = — = 10
q
Example 3
Given h = lOft
0
U = 0. 1 ft/sec
0
T
0 -4 2
g— = 10 ft/sec
p 0
= 10 ft /sec
K = a x 10 ft/sec
This case differs from examp’e 2 only in that the value of K is doubled.
We first calculate
U 2
= T = 10, k = -= 2x10 4 , 1/B = 1O
g— h
p o
0
Thus kR = o. 2. Referring to Fig. 4. 8, we see that an internal hydraulic
jump would develop. The program SBJ2 is thus first used with F 0 = 10,
k 2 x 10 , and 1/B 10 . Portions of the output is tabulated as follows.
Also the quantity s = is calculated and inserted as the last column.
a
x/h h /h F s
2o a
1.0 4.03 0.175 1.24
6.0 4.08 0.262 1.22
7.0 4.07 0.278 1.23
8.0 4.07 0.294 1.23
The values of F 2 and s when referred to Fig. 4. 7 shows that the jump
would occur at about x/h = 7 or about 70 feet from the source. At that
point, the excess temperature is about 80% of that at the source and the
113
-------
flow rate about Z5% more than that at the source. Thus we deduce that
after the jump, F 2 0. 278, h 2 = 41 ft, U 2 = 1.25 0. 03 ft/sec. To
obtain the flow field after the jump, we may use SBJ2 again with F 0 = 0. 278,
k 2 6.7 x 10 , hR = = 8 x 10 where now the output
quantities are normalized to the values just after the jump and x is now
measured from the jump location.
114
-------
CHAPTER 5 PASSIVE TURBULENT DIFFUSION FROM A CONTINUOUS
SOUR CE IN A STEADY SHEAR CURRENT OR AN UNSTEADY
UNIFORM CURRENT WITH UNSTEADY SURFACE EXCHANGE
5. 1 Introduction
The mixing process undergone by the cooling water discharged from a power
plant can be divided into three stages, as discussed in Chapter 3. In
the case when the discharge is from a submerged outfall, these three
stages are (Fig. 5. 1):
1. An initial mixing stage governed by the momentum and
buoyancy of the discharge.
2. An intermediate stage of dynamic spreading governed by
the buoyancy and the density stratification.
3. A final stage of essentially passive turbulent diffusion.
In this chapter we shall treat the third stage of passive turbulent dispersion
when the dilution and dispersion are primarily governed by ambient tur-
bulence and currents. This phenomenon is very similar to the case of
the mixing of sewage effluent from a submerged ocean outfall.
Two mathematical models will be established in this chapter for the cal-
culation of the distribution of excess temperature (or dilution of the effluent)
due to the effects of ambient turbulence and current, and surface heat ex-
change. The first model is for the case of the steady state passive turbulent
diffusion from a continuous source in a steady shear current as shown in
Fig. 5. Za. The second model is for an unsteady case where the current
and surface exchange coefficients are time -varying as shown in Fig. 5. 2b.
However, in the second model the current is taken to be uniform with
depth.
Note that the location of the source can be below the water surface such as
in the case of an effluent field trapped by ambient density stratification.
115
-------
Dynamic spreading Passive turbulent diffusion
u
Buoyant jet
a) Effluent field established at water surface.
I
Dynamic spreading Passive turbulent diffusion
Buoyant jet mixing
, I 1 1 / T 1 I 1 I , I I 1 1 1 / 1
b) Effluent field trapped below water surface
Figure 5. 1 Various stages of mixing of submerged cooling water
discharges. Effluent field established at water surface.
116
-------
Constant heat exchange K Constant
e
y
x
/
y h ,—Reflecting boundary
b / (no exchange)
/ / / / / / / / I / / / I / / I / / / 77
a) A steady continuous source in a steady shear
current with constant surface heat exchange.
y = hb
/ / / / / / / / /
r
/7 -
R eflecti.ng boundary
/ / / 7 / J 1 / / /
Figure 5. Z
b)
An unsteady continuous source in a uniform
unsteady current with time varying surface
heat exchange
Passive turbulent diffusion cases studied.
x
z
u(t)
y
117
-------
If this is the case then the excess temperature over the ambient will be
practically zero unless there is salinity difference. Since the models to be
developed are applicable not only to the prediction of excess temperature
distribution, but also to other water quality indicators, the case of sub-
merged source is also of practical interest.
118
-------
5. 2 Derivation of Basic Equations
The basic relation governing turbulent passive diffusion is a conservation
equation for the diffusant: either heat content or a tracer concentration.
The mixing is assumed to be completely dominated by the ambient current
and turbulence characteristics. Thus, before the problem can be solved,
the environmental conditions must be known.
The models to be developed in this chapter are equally applicable to the
dispersion of heat or any other tracer (such as the concentration of coli-
form bacteria). However, there are some differences and therefore they
will be discussed separately.
5. 2. 1 The Problem for Excess Heat
The equation of conservation of heat content without internal heat source
is:
+u +v—+w—
= - -- (K — ) + (K —k) + — - (K t) (5 ‘)
x x x y y y z Z Z
where H = the total heat content above a given reference
heat content;
t = time;
x, y, z = coordinates in longitudinal, vertical and transverse
directions (see Fig. 5. 2);
u, v, w = velocities in (x, y, z) directions;
KX KY KZ = exchange coefficients (eddy diffusivity plus mole-
cular diffusivity) in (x, y, z) directions.
In a natural body of water, the vertical velocity v is usually small (except
near zones of strong upwelting and sinking flow). In this model, v is
taken to be zero. In addition, it is a reasonable assumption that the flow
119
-------
is predominantly in one direction, say x. Thus w is taken to be zero.
Equation (5. 1) then becomes
+ u = —(K-- ---) + —(K - --- -) + (K (5. 2)
The mechanism for surface heat exchange can be expressed as (Edinger
and Geyer, 1965)
- K = - KE(T - E) (5. 3)
at y = 0, i. e., the water surface; where
KE = surface heat exchange coefficient;
T = surface water temperature;
E = equilibrium temperature.
At the bottom, the exchange of heat may be taken to be zero, i. e.
K 0 (5.4)
y y
at y = hb, i. e., the bottom.
The equation governing heat exchange processes in the ambient without
any waste heat addition is simply:
a a a
+ u— = —(K ) (5.5)
y y y
where H = heat content of the ambient water.
a
Note that in Eq. (5. 5) inhomogeneities in the horizontal directions are
neglected. In general this is a good approximation for a body of water
which is not too vast. The boundary conditions at the surface and bottom are
similar to Eqs. (5. 3) and (5. 4), i. e.
120
-------
K a KE(Ta - E) at y 0 (5. 6)
and
= 0 aty=hb (5.7)
where Ta is the surface temperature of ambient water.
To obtain the relation for excess heat distribution, we simply subtract
Eq. (5. 5) from (5. 2) to obtain
. H
— + u— = — (K —) + —(K —) + —(K —) (5. 8)
x x y y y z
where H = H - H is the excess heat content due to waste heat addition.
t a
The surface and bottom boundary conditions are then:
K = KE(T - Ta) at y = 0 (5. 9)
and
= 0 atyhb (5 . 0)
Note that KE and E are assumed to be the same in Eqs. (5. 3) and
(5. 6). For the temperature range encountered in the passive turbulent
diffusion stage, this is believed to be an adequate assumption.
Defining the excess temperature to be T, and assuming a constant specific
heat Ch over the temperature range of interest, Eqs. (5. 8) to (5. 10)
become
+ u I = - --(K - ) + — -(K I) + — --(K - ‘) (5. 11)
X X X X y y y z z z
K
K = E TK T atyO (5.12)
y y pCh e
121
-------
and
K - = 0 aty=hb (5.13)
where Ke KE/(pCh)
p = density of water
5. 2. 2 The Problem for a Tracer
We now consider the corresponding problem where the dispersing sub-
stance is a tracer which is otherwise not present in the ambient water.
The conservation relation is:
+ U + v + w = - - (K - -) + (K - ) + - -- (K c ) -K c (5. 14)
x x x y y y z z z d
where c = the concentration of the tracer
Kd = the decay coefficient of the tracer or the die-off rate.
(Clearly, Kd = 0 for conservative tracers such as dye.
Again, we assume v = w 0. Equation (5. 14) becomes
+ u - - = - -(K - -) + — --(K - -) + — -(K - -) - K c (5. 15)
x x x x y y y z z z d
For tracers, such as dye, salt or radioactivity, there is no surface or
bottom exchange; hence
= 0 at r = 0 and r = hb (5. 16)
The basic equations and boundary conditions are very similar between the
cases for excess heat and for a tracer. In fact, both problems can be
included in a single more general mathematical model. This is discussed
in the next section.
122
-------
5. 2. 3 The General Problem
The general equation governing either excess temperature or a tracer
substance (based upon Eqs. (5. 11) and (5. 15)) can be written
+ u - - = — --(K - -) + -(K -) + -(K - -) - Kdc (5. 17)
The general boundary conditions can be written
K—Kc aty =O
y y e
(5. 18)
= 0 aty=hb
In these equations, c is either the excess temperature or the tracer con-
centration. For a tracer without surface exchange, then Ke = 0; while
for the case of excess temperature or a non-decaying tracer Kd = 0.
Note that the problem is not yet posed completely. To fully define the
problem, the initial condition, and the source condition must be specified
along with the environmental, characteristics such as u, K , K and K
x y z
These will be discussed as we treat the two models separately in the
following sections.
123
-------
source thickness
h ________________________________
source _____ /
level O ‘ ‘
f
current
y
y = hb
/ 1/ / / I / / / / / / / / / I / /7/
Elevation
z
source
width L=
0 X
Plan
Figure 5. 3 Flow configurations in cases with steady releases.
a) Submerged release.
I 24
-------
source 4 y=o
y,- x
thickness h /Z
Current
y
y = hb
/ /1 / / / /1 / 1/ / / / / / / / /
Elevation
z
Source
Width L= - -
Plan
Figure 5. 3 Flow configurations in cases with steady releases.
b) Surface release.
125
-------
5. 3 Steady Release in a Steady Environment
5. 3. 1 Formulation
For the case of steady release of waste heat or a tracer substance into
a steady environment, the governing equation (5. 17) becomes:
- - -(K - -) + — -(K - -) + - -(K - -) - K .c
u -
(5. 19)
in general, the term representing the longitudinal transport —(K - ) is
small in comparison with the transverse transport term _(K -). Thus
the longitudinal transport is neglected. Equation (5. 19) can then be written
as:
u = - -(K - ) + — -(K - -) - K c
y y y ? z z z d
The boundary conditions are of course, still given by Eq. (5. 18).
5. 3. 1. 1 Source Conditions
(5. 20)
The source will be taken to be located at x =
thickness h and width L as shown in Fig.
release, y y = 0 and the thickness is h/2
The distribution of c at the source is taken
c(o, y, z) = c(o, y) exp { -
0, at a depth y = y with
5. 3a. For the case of surface
as shown in Fig. 5. 3b.
to be:
where c(o, y) = c(o, y 0 , o) (Note: c(x, y) = c(x, y, o).)
Equation (5. 21) defines the source distribution to be Gaussian in the
z-direction, and resembles an ellipse in the y-direction. This is the
assumed distribution used in the examples in developing the model. The
z 1
a 21 2(y _ y 0 ) 1 2 1 J (5.21)
{ h
0
126
-------
actual distribution, if known, should replace Eq. (5. 21) in a practical
application.
5. 3. 1. 2 Environmental Characteristics
The environmental conditions relevant to the passive turbulent diffusion
phase as formulated above include the ambient turbulence characteristics,
represented by values of the eddy diffusivities K and K; the ambient
current structure, represented by u(y); and the surface exchange coefficient
KE. The current is a directly measurable quantity. Given the general
locale of the discharge, field data on u(y) can be gathered and used in the
prediction model. The other quantities, namely K, K and KE are more
difficult to measure. In general, no direct measurements are made and
the values are usually inferred from other observable phenomena.
It is not the intention of this study to develop a detailed method of estimating
KE accurately. Studies on the heat transfer between the water environ-
ment and the atmosphere have been initiated and are being extended by
other investigators. Strictly speaking consideration of the heat exchange
processes at the water surface is very complicated. Factors which are
of importance include solar radiation, back radiation, conduction, con-
vection and evaporation. The introduction of the coefficient KE and the
equilibrium temperature (see Edinger and Geyer, 1965), lumps the effects
of all these mechanisms of heat transfer together. At present, this appears
to be the most practical method available. As better relations become
available, it should be possible to modify the models developed accor-
dingly.
The eddy diffusivities K (horizontal) and K (vertical) are important in
determining the dispersion of the heated effluent. Like the coefficient KE,
these diffusivities are also empirical coefficients which depend on more
basic phenomena such as the turbulence structure in the fluid medium.
In turn, it can be expected that the turbulence structure depends on the
input of energy from the atmosphere through wind and waves, the density
127
-------
stratification or stability of the fluid medium and the current shear which
can supply the energy in generating turbulence. These diffusivities will
now be briefly discussed in the following subsections.
5. 3. 1. 2a Horizontal Diffusion Coefficient K
z
In a large body of water, the horizontal diffusion coefficients are generally
governed by the ‘ 4/3 power law T , i. e. , the horizontal diffusion coefficient
is proportional to the 4/3 power of the length scale of the diffusing patch
or plume:
Kz = ALL ’ (5. 22)
2/3 2/3
where AL as a dissipation parameter (cm /sec or ft /sec.
L is the width of the plume (usually taken to be 4c 2 being
the standard deviation of the concentration distribution).
Equation (5. 22) can be written in terms of
K = (5.23)
It should be noted that use of Eq. (5. 23) results in a nonlinear governing
equation.
In the ocean, numerous field experiments have been performed to estimate
K. This is summarized in Fig.5. 4.It is seen that the value of AL is
in the neighborhood of 10 - 10 ft /sec. Thus the value of A is
in the neighborhood of 10 - 6 x io2 ft2/3/sec.
It should be pointed out that no effect of shear currents were removed in
the field experiments so that direct use of the data requires some caution.
It is believed that since the effects of shear is explicitly taken into account
in the present model, the lower value of A = 10 ft2/3/sec might be more
appropriate.
128
-------
10 ’ 1 — —— ________
/
/ I,
4 / ,
/ I, -
G
0.000 15L /3 _\\
I,
o p
V
/
• M U fJERSON I 0
D}
/
“.. 5 C
1’
N
/
/
0
X 4
JO K, 0.001 “1 \
//. /
/ ./ ,
L’M ITS CF DATA
Z / OLSON AND CHIYE
/ /• } RIFT CARDS
— NRDL
I L. ?
U. 10 DEEP LAYER
“c: ,crv1’’r IN
w EXPERIMENT ’ A ,/’/’ £c L COON
f ,o,, /
o 1968 Id— I ‘J LE OUTFALL FIELD
o a,
CQ’LCURRENT-CJ OSS P .R
JO- /
z I’D CRLDB, S 0.00327, I” MESH
o K : 0005 L ’ _ • J /, //
• CRLOB, S OCCC5S. I MESH
0 SVERDRLP \
D -
I L . -‘ a/DYE () PRQUCMAN J
U. “0
/ TcEA1AS Q MUNK
— ,/ I PEARSON
0 4
0 I K, 000 / H NZAVI4 U MARY BY
/ 00
/ , / co , STOMMEL E R SON-
Z . •MCON,ET AL I
I- ‘ xJ4:// /1
7 L ’ 4 ” / / 0 EIHARLEMAN
í a OVON ARX I
o
I
,o_ 3 / x GUMNERSON I
+ HI AKA )
I V . /
/ PARKER - 1 561
r ’
Kr I I 1 I
0; 0 102 IO Q4 IO 0’ ‘o
HORIZONTAL SCALE OF DIFFUSION PHENOMENON , L-ft.
Figure 5. 4 Horizontal diffusion coefficient as a function of
horizontal scale (from Orlob (1959)).
129
A -2-890
-------
It can also be observed that the data indicates A to be larger when the
scale is smaller. This is reflected in somewhat larger K values
for small L in Fig. 5.4 than those corresponding to A = lO ft /sec.
For our purposes, it may be assumed that A = iO ft2/3/sec. When
more field data is available the model may be readily adapted to take
advantage of them.
The experimental data summarized in Fig. 5. 4 are all for the case when
the diffusing pool is on the ocean surface. A few experiments have also
been performed for the case when the diffusing pooi is at depth (Parr 1936;
Riley, 1951; Ozmidov, 1965; Munk, Ewing and Revelle, 1949; Kolesnikov,
Panteleyev, and Pisarev, 1964; Snyder, 1967; and Schuert, 1969). The
results generally indicate somewhat smaller values for the diffusivity.
It is generally recognized that the presence of a stable density gradient
damps out vertical turbulent fluctuations and hence vertical turbulent
transport. However, conflicting views exist for the effect of stability
on horizontal transport. The main difficulty is the lack of adequate field
data.
On the one hand, Parr (1936), Riley (1951), and Ozmidov (1965) suggested
that the horizontal diffusion coefficient increases with stability. Moreover,
they attempted to verify the hypothesis: Parr by using data on the distri-
bution of Atlantic Ocean waters flowing into the Caribbean Sea, while Riley
by analyzing salinity and temperature distributions in the ocean. On the
other hand, Munk, Ewing and Revelle (1949), Kolesnikov, Panteleyev,
and Pisarev (1964), Snyder (1967), and Schuert (1969) suggest the opposite,
i. e., the horizontal transport decreases with stability. Munk, et al. found
that in Bikini Lagoon at 50 meter depth, the value of horizontal diffusion
coefficient was only one-third of that near the surface. Kolesnikov found
by direct measurements, AL to be 0. 01 cmZ ’3/sec. at the surface and
0. 0046 cm2/3/sec. at 500 meter depth. Snyder found that at 9 foot depth,
the value of AL dropped to one-quarter of the value at the surface.
Schuert found that at 300 meter depth, the value of AL is about one order
of magnitude smaller than in surface waters.
130
-------
It can be seen from the above discussion that the dependence of K on
stability is still controversial and unsettled. The diffusion model to be
developed herein, however, can be modified to incorporate a K as a
function of y, the vertical coordinate once a reliable relation is established.
5. 3. 1. Zb Vertical Diffusion Coefficient
In contrast with the relative abundance of data for horizontal diffusion
on the ocean surface, there is a scarcity of data for vertical diffusion.
Evaluation of vertical diffusion coefficients have typically been implicit,
e. g., based on the temperature and salinity distribution and their time
and space variations. The matter is further complicated by the wide
-a
spread in the measured values, (from as low as 4 x 10 to as high as
200 cm 2 /sec). Moreover, no obvious relations were available between
the vertical diffusion coefficient and other readily measured parameters.
The presence of density stratification tends to supress vertical exchange.
Therefore, one expects the vertical diffusivity to be a monotonic de-
creasing function of density stratification. The presence of shear tends
to be destabilizing and increases vertical exchange. It can be expected
that for similar flows the vertical diffusivity should be a monotonic
non-increasing function of the Richardson number defined as
R - p dy
i — (du)Z
dy
Numerous proposed relations between K and R. are summarized in
Table 5.1 . Unfortunately, these cannot be checked and the constants ( )
cannot be readily determined due to a scarcity of data on the shear.
Rather than relating K to R. which is the physically more logical approach,
it is also possible to attempt a correlation of K with
p dy
0
the density gradient alone. Strictly speaking, one would not expect a one-
to-one relationship to exist between K and c.
131
-------
TABLE 5. 1
Summary of Formulas on Correlation of Vertical
Diffusion Coefficient K with Pichardson T s Number
y
R. (or Density Gradient e)
e., the neutral case : proportionalility constant
* As given by Okubo (1962)
** As given by Bowden (1962)
*** The formulas presented in the translated version are apparently erroneous.
132
NOTE: K K at B. = 0, i.
yO• y 1
varies from case to case
Rossby and Montgomery
(1935)*
K =K 0(1 +
• y 1
Rossby and Montgomery
(1935)’
K = K 0 (1 +
1 1
Holzman (1943)*
K K 0 (1 - I) .
Yamamoto (1959)*
K = K 0 (1 - p lf2 B. 4
Mamayev (1958)*
K = K e B1
y yO
MunkandAnderson
(1948)**
K =K (1 +$
y Y° 1
3. 33 based upon data by Jacobsen (1913)
and Taylor (1931)
Harremoes (1968)
K = 5 x 10 3 x 2 ’3 cm 2 /sec
note: C in m’; approximate experimental
range 5x10 9
-------
All readily available data on K where e is simultaneously measured
are collected and plotted as shown in Fig. 5. 5. It can be observed that
almost all data fall within a factor of 10 of the empirical relation
-4
K = (K in cm sec; in m ) (5. 24)
y y
-7 -1 -2 -1
4x10 m E l0 m
It should be noted that Fig. 5. 5 contains only those available data where
both K and e are available.
y
The relation K = 10 /c can be deduced from the definition of the vertical
y
diffusion coefficient provided some assumptions are made. We assume that
diffusion occurs due to turbulence so that molecular diffusion may be ignored.
With the assumption that the density variations are small in the ocean, the
equation describing the variation of potential density is then
+ u + v + w — -- (K + — - (K - ) + (K P)
Z X X X y y y z z z
where u, v, w are the mean currents in the x, y, z directions respectively.
Since horizontal variation of p are usually much smaller than vertical
variations, we assume = = 0; also, we will assume v = 0. Then
the equation becomes
- = - -(K
t y y y
The term is usually very small except for the near surface waters
which may undergo some diurnal changes. Thus we assume steady state.
Then
- -(K - ) = 0
y y y
or
133
-------
ri J I
102
0
}
X
KOLESNIKOV (1961)
X
A
HARREMOS
0
(1967)
JACOBSEN
0
K
FOX WORTHY
(DEFANT, 1961)
U
0
XX
(1968)
FOXWORTHY
£
PATCH
0
(1968)
FOX WORTHY
K
U
m
N
E
U
w
U
Li.
C-)
z
0
U)
L i .
0
-J
4
C-)
L i i
>
0
PLUME
K
U
0
0
(1968)
0
K
POINT SOURCE
U
0
A
K
a
0
a
a
A
a
0.1
A
0 •
S
0
0
.
a
A
0.01
I0
‘0 - I
£ •-I/p dp/dy ,DENSITY.GRADIENT,m’
lO-z
IO..I
Figure 5. 5
Correlation of Ky with density gradient.
-------
K = constant
y y
Hence
K - Constant
This is exactly the form of the relation between Kand € found depicted
in Fig. 5. 5.
It is proposed that unless independent field data is available, Eq. (5. 24)
be used to estimate K for application of the present model. It is realized
that this is approximate at best. However, it is, at present, the most
rational method available. Future studies and measurements may alter
this relation. The region of applicability of this relation is l0 <
-2 -1
10 m
In the surface mixed layer of the ocean, the density gradient is often zero.
The empirical relation is certainly invalid since it implies an infinite
K . In this case, the vertical transport is governed primarily by the vertical
turbulence created by waves and wind. Relations between the vertical
diffusion coefficient K in the mixed layer and the surface wave char-
y
acteristics have been proposed by Golubeva (1963) and Isayeva and Isayev
(1963). Their relations can be summarized by
H 2
K = 0.02 w
yl T
w
where K = vertical diffusivity aL the surface
yl
H = wave height
w
T = wave period
w
Thus given the sea state, K 1 can be estimated. Fig. 5. 6 shows the
relation between K and sea state.
yl
135
-------
4’
SEA STATE
Figure 5. 6
Dependence of Ky 1 on sea state
U
U)
(‘J
E
()
200
100
0
0
2
3 4 5
A-Z- 1643
136
-------
In summary, data on vertical diffusivity are scarce. Although logically,
K is expected to depend on } ichardson number, for practical purposes,
the empirical relation K 10 1€ is proposed, subject to modification
as better data become available.
In general, K has its maximum value in the surface layer: in the open
ocean K at the surface varies between 10 - 200 cm /sec. ; in coastal
2 2
areas, 10 - 50 cm /sec.; in lakes, 10 cm /sec. Below the surface
mixed layer (or epiliminion) K drops to its minimum in the thermocline
2
(of the order of 1 cm /sec. in the open ocean; in lakes, K may drop to
2
as low as 0. 05 cm /sec. ). Below the the rmocline, K may increase
again. Some typical values and K - profiles in lakes and reservoirs
have been determined by Orlob and Selna (1970).
5. 3. 2 Method of Moments
The problem posed in Sec. 5. 3. 1 is complicated and cumbersome to solve
due to the three independent variables and the complexity of the coefficient
functions. This difficulty can be partially overcome by using the method
of moments.
Define moments of the distribution by:
c(x, y) = Sc(x, y, z) dz (5. 25)
c 1 (x, y) = Jzc(x, y, z) dz (5. 26)
c 2 (x, y) = z 2 c(x, y, z) dz (5. 27)
The zeroth moment c is the integrated amount of excess temperature or
tracer in the z-direction. The first moment c 1 is related to the z-
coordinate of the centroid of the c-distribution. In the present model,
137
-------
it is zero because of the symmetry of the distribution in the z-direction.
The second moment c 2 defines the spread in the z-direction. The width
of the effluent field is usua fly taken to be 4a where a 2 = c /c
z z 2 o
Multiplying Eq. (5. 20) by z°, z 2 , and integrating over z, we obtain the
equations governing the moments:
u = — -(K - K c (5.28)
y y y do
u 2 = — --(K — -) - K c + 2K c (5. 29)
y y y d2 zo
An alternate to Equation (5.29) can be written in terms of a as:
2 2 2
u Z = — --(K Z + 2K .i. __a Z + 2K (5.30)
y y y yc 0 y y z
The boundary conditions expressed in terms of the moments corresponding
to Eq. (5. 18) are:
K = K c
y?y eo
aty=0 (5.31)
K
y y e2
or 2
K Z = aty0
y y
and
2
0 = _ = Z = 0 atyh (5.32)
b
The source condition expressed in terms of the moments corresponding to
Eq. (5.21) are:
138
-------
( r Z(Y-y) 2
c 0 (o, y) = c 0 (o, y 0 ) 1 - L h] J 33)
3/4
2 r Z(y-y 0 ) 1 2
c 2 (o, y) c 0 (o, (o, y) { 1 - L 0 J
2 2 r
or = 1’ - L h ] (5.
Thus by the moment method, the number of independent variables is reduced
by one. Since the c-distribution in the z-direction is usually found to be
of Gaussian form both from field and laboratory experiments, the diffusion
process can be adequately described by knowing the zeroth and the second
moments. In fact, if c is exactly Gaussian in the z-direction, then it is
completely specified by its zeroth and second moments. Equations for
higher moments can be formulated in a similar way. These are only
necessary if the c-distribution is distinctly non-Gaussian in the z-direction
In that case, for example, the third moment would indicate the skewness
of the distribution.
5. 3. 3 Limiting Solution for Cases with Zero Vertical Transport
If the effluent field is trapped within a strong thermocline where vertical
transport is small (for example, if K is close to the molecular value),
a good approximation to the solution can be achieved by taking K = 0.
The proper criterion for the validity of this approximation is that the
vertical spreading over the distance of travel is small in comparison with
the vertical dimension of the source, i. e.
/K << h (5.36)
V yu 0 o
where x is the horizontal distance of interest;
is the characteristic velocity; and
h is the source thickness.
0
139
-------
for K = 0. 1 cm 2 /sec., u = 15 cm/sec. (or 0. 5 fps) and x = 1,500
meters, / Kxt/u = 0. 3 m; therefore, for sources thicker than, say
10 meters, it will be sufficient to consider it as a case with zero vertical
transport.
For cases with zero vertical transport, an analytical solution can be ob-
tained as follows. Equation (5. 28), when K = 0, becomes
u = - Kdc (5. 37)
Integrate with respect to x, noting that u is only a function of y:
K
d
c(x, y) = c(o, y) exp [ - —x (5. 38)
Equation (5. 30) becomes:
u = 2K when c 0 (5. 39)
z 0
Applying the 4/3 power law for K:
K =
Eq. (5. 39) becomes
2
u = 2A (5.40)
z
Integrating with respect to x, Eq. (5. 40) gives:
2/3 2 A 2/3
(x, y) - —x + °, y) (5. 41)
For a given set of source conditions, i. e. , c(o, y) and a(o, y), environ-
mental characteristics u(y), A, and the decay coefficient Kd. solutions
are given by Eqs. (5. 38) and (5. 41) if the vertical transport can be neglected.
The maximum concentration or excess temperature (assuming c is Guassian
in z) is given by
140
-------
c(x, y) (o, y)
Cmax(Xj ) = cmax(Oz c(o, y) :(Xi y) (5. 42)
5. 3. 4 Dimensionless Equations and Numerical Solutions
For general cases with non-zero vertical transport, it is necessary to
resort to a numerical method of solution unless the environmental and
source conditions are very simple.
Before the numerical solution is attempted, the governing equations,
source conditions and environmental conditions will first be normalized
by defining dimensionless variables (quantities with primes) as follows:
Coordinates:
x l = X/X (5. 43)
y l = y/ fR xe/u (5. 44)
Velocity: u’ = u/u 0 (5. 45)
Vertical diffusion coefficient: K 1 = K /K (5. 46)
y y yo
Dissipation parameter: X l = Ax / [ 2/3(0, y)u (5. 47)
Exchange coefficient: Ke = KeV’XtIKyoUo (5. 48)
Decay coefficient: K 1 d = Kdxt/u (5. 49)
Zeroth moment: c’ 0 = c/c 0 (o, y 0 ) (5. 50)
Second moment: c 1 2 = c 2 / [ c 0 (o y 0 ) a 2 (o, y )J (5. 51)
141
-------
Lateral spreading: a = Ia(o, y) (5. 52)
Maximum concentration: c i = c Ic (o y ) (5. 53)
max max max ‘ o
&
NOTE: & =
max 0’
z
Here x is the terminal x value of interest;
K is a characteristic vertical diffusion coefficient; and
yo
u 0 is a characteristic current speed
For example, the characteristic values K and u can be taken to be
yo 0
their values at the free surface if the effluent is at the free surface.
The governing equations, (5. 28), (5. 29) and (5. 30), in dimensionless
form becomes:
u’ = — ---(K ’ ‘)- K’ c’ (5.54)
y’ y y d o
= — - 1 (K ’ ) + 2K’ c’ - K’ c’ (5. 55)
y y y z o d 2
2 2 ,
u’ Z = — - (K’ , ) + 2K’ + 2K’ (5. 56)
y y y yc 0 y y z
The normalized boundary conditions corresponding to Eqs. (5. 31) and (5. 32)
become
K’ = K’ c’
y y e 0
aty=0 (5.57)
K’ , = K’ c’
y y e 2
or
142
-------
a
K = 0 y=O
and
a
= 0 at y = h’b (5. 58)
where
h b = hb/ /i xt/u
The normalized source conditions based upon Eqs. (5. 33Yto (5. 35) are:
I Z(y’-y )
c’(o,y’) - J J
Z(y’-y’) Z 3/4
c’ 2 (o, = - [ h’ j (5. 60)
and
a Z(y -y ) Z
= - [ h’ 01 J (5. 61)
h’
for y’ 0 and y 1 - y’ + h’/2
where
= y / /K x Lu , h 1 = h /,JK x Lu
0 0 yot 0 0 0 yot 0
Note that Eqs. (5. 54) to (5. 58) are identical in form to their corresponding
dimensional equations. The main effects of this normalization are: 1) the
region of interest in x is normalized to 0 x 1; a) the source condition
is normalized; and 3) the u- and K—Profiles are normalized.
Thus, the problem of the steady state distribution of excess temperature
(or tracer) resulting from a continuous source in a steady but non-uniform
environment is formulated in dimensionless form. A computer program
143
-------
(PTD) has been written based on the Crank-Nicolson Method and is included
in Appendix C. Given the input conditions h’, ‘ ‘ K- and u-profiles,
K’e Ktd, h’b, the problem can be solved by using the program.
Before discussing the example solutions obtained, we shall first choose
the various parameters and parameter functions to specify the problem.
The following values have been chosen as representative typical values:
A = l0 ft /sec.
x = 10,000 ft.
K = 10 ft. /sec. (surface)
yo
u = 1 ft/sec. (surface)
0
hb = 100 ft.
h = 2Oft.
0
° ‘ y) = 30 ft.
Ke = 10 ft/sec.
From these,
X I = Axt/ [ a 21 (o, y) u] 10
K’ = K /x /K u = i0
e e t yoo
h
= b = 10
b .JK x/u
yo t 0
h
h ’ = 0 = a
0
JK x/u
yo t 0
In the following discussion, the primes will be dropped for simplicity.
144
-------
The parameter functions K(y) and u(y) are chosen to be as shown sche -
matically in Fig. 5.7. The constant parameters K1’ K2’ K3’ K4’
would specify the dimensionless Ky - profile while the constants
U would specify the dimensionless u-profile. It should be pointed out
here that is more or less an arbitrary number. The program PTD
can be run for x from 0 to any value, not necessarily 1. Also, by proper
choice of y the source may be located at the surface ( y 0 = 0) or at any
depth (y > 0).
Guided by the numerical values mentioned in the preceeding paragraphs,
a total of 14 cases has been computed using the program PTD. The
parameters and parameter functions chosen for each case are summarized
in Table 5. 2. As can be seen from the table, two different profiles were
selected for Ky and u as functions of y. The first is when it is constant
with depth. The second is when it takes on a shape judged typical of situations
when the ambient is density stratified. The identification code is designated
by two letters followed by three numbers. The first letter signifies whether
the Ky-profile is constant (C) or not constant (S or T). In case it is not
constant, S stands for the case when the source is at the surface and T the
case when it is submerged (in the thermocline). The second letter signifies
whether the velocity profile u(y) is constant (C) or not (N). The first
number, n , refers to the value of X: X= 1 corresponds to n = 1 and
X = 10 to n 1 = 2, thus, X = 10 ‘ 1 . The second number n 2 refers to
the value of K , by the relation K = 10 2• The third number n refers
e e _n 3
to the value of Kd by the relation Kd = 10
It should be noted that c , the zeroth moment of the distribution is in-
dependent of X. Figures 5. 8 a, b, and c show c 0 (y) plotted versus y
for various values of x. Several different cases are shown on the same
graph to delineate the effect of various parameters. Figure 5. 8a is for
the case when K is constant; Fig. 5. 8b for the case when K is not
constant and the source at the surface while Fig. 5. 8c for K not constant
and source submerged. The effects of Ke and Kd and the current profile
can be observed by comparing the cases in each figure. The effect of
145
-------
1
u-profile
Figure 5. 7
Profiles of Ky and u used in study.
Ky
Kl
K2
K3
Bz
Ky-profile
U
I
146
-------
Ky
U
Ke
ID
0 1
Ill/I/ I) y
5 Ky
10
fl //If/I
V
I//f
I/If
1
10
10
0
0
10
0
0
0
0
0
CC
100
CC
200
CC
220
CN
200
0 1
(f//I/ I )
5
10
H
I//I
V
/1/I
1
10
10
0
0
0
10_i
0
0
0
0
SC
100
SC
200
SN
200
SN
210
0 1
5
10
H
I//I
V
I//I
1
10
1
10
0
0
0
0
1O
0
0
0
l0
0
0
0
TC
100
TC
200
TC
202
TN
100
TN
120
TN
200
Table 5. 2 Summary of PTD cases. (Ky = Vertical diffusivity; U = Current; = Dissipation
parameter; Ke Surface exchange coefficient; Kd = Decay
coefficient) All quantities dimensionless.
-J
-------
K-profile can be seen by comparing Fig. 5. 8a with 5. 8b. As can be seen
from these figures, the effects of the various parameters and parameter
functions are as expected. For example, K tends to decrease c but
primarily at the surface, and the current shear tends to promote some-
what higher dispersion.
The effect of \ on the solution is on the spread of the plume or in the value
of the second moment c 2 . Figures 5. 9 a, b, and c show the solution
= Jc Tc plotted versus x for y = y 0 . It is readily seen that when
X goes from 1 to 10, the spread at x I increases about ten times.
This is not surprising since X is proportional to the horizontal diffusion
coefficient. The effect of shear on can also be observed to promote
a somewhat larger value of as would be expected. Comparison of
Figs. 5.9 a and b with 5.9 c shows that the effect of shear is correspondingly
larger when the source is submerged than when it is at the surface. This
is because the value of the velocity at y = 3 for the TN - runs is only 0. 7
times that for the TC-runs.
From the above discussion it is seen that the model developed did not yield
any profoundly different results than what can be reasonably expected. In
any practical situation, the parameters and parameter functions and the
source conditions may be different from those chosen. The program can be
readily modified to suit those conditions.
It should be reiterated here that the model which resulted in the program
PTD is for the case of a steady discharge into a steady environment. In
practice, the parameters K, u, and the source intensity are most pro-
bably not constant in time. In such a case, the program PTD should not be
used. In the next section of this report, a limited unsteady case will be
treated and discussed.
148
-------
0
C 0 (x,y), NORMALIZED ZEROth MOMENT
Figure 5. 8a
Vertical distribution of c (x, y) for PTD cases
CN 200, CC 200, CC lOO,° CC 220 (Ky-profile
uniform).
X=O.I
/
/
/
/
w
I-
4
z
0
0
0
()
-J
4
0
I-
U i
>
0
U i
N
-J
4
0
z
>1
I
2
3
4
5
6
7
CN200)
(CC200,
CC 100)
(CC220)
0 0.5
149
-------
0
I
2
3
4
5
U i
I-.
z
0
0
C-)
-J
C .)
I-
U i
>
Ui
N
-J
0
z
Figure 5. 8b
C 0 (x, y), NORMALIZED ZEROth
MOMENT
Vertical distribution of C 0 (x, y) for PTD cases SN 200,
SN 210, SC 100, SC 200 (Ky-profile not uniform,
0
0
0.5 1
surface release).
-------
X=I.O X:O.5
X:O.I X=O
C 0 (x,y), NORMALIZED
ZEROth MOMENT
0
I.0
2.0
3.0
4.0
LrI
w
z
0
0
C -)
-J
4
C-)
w
>
0
w
N
-J
4
0
z
5.0
0
0.5
1.0
Figure 5. 8b
Continued.
-------
0
TCIOO, TC200)
C 0 (I,y), NORMALIZED ZEROIh MOMENT
Figure 5. 8c
Vertical distribution of c 0 (x,y) for PTD cases
TN 100, TN 120, TN 200, TC 100, TC 200, TC 202
(Ky..profile not uniform, subsurface release).
L i i
I—
z
0
0
0
()
-J
C.)
L i i
>
0
Lii
N
-J
0
z
I
2
3
4
5
6
TNIOO, TNI2O, TN200)
TC202)
0 0.5 1
AT X=l
1 52
-------
x, NORMALIZED HORIZONTAL DISTANCE
Figure 5. 9a
Width of diffusing plume for PTD cases CC 200,
CN 200, CC 100, CC 220 (Ky-profile uniform).
0
I,
I-.
_ 4
U I
U) b
C N 2
10
I
CC200, CC220
CC 100
0.01 0.1 1
-------
I0
SN200, SN2IO•
x, NORMALiZED HORIZONTAL DISTANCE
Figure 5.9b
Width of diffusing plume for PTD cases SN 200,
SN 210, SC 100, SC 200 (Ky-profile not uniform 1
0
I’
4
SC200
1
0. 1
I
surface release).
-------
X, NORMALIZED HORIZONTAL DISTANCE
Figure 5.9c
Width of diffusing plume for PTD cases TN 100,
TN 120, TN 200, TC 100, TC 200, TC 202
(Ky-profile not uniform, subsurface release).
TN 200
U
4
bN
I0
TN 100
TN 120
.01 .1 I
-------
5.4 Continuous Belease of Heat into a Uniform But Time-Varying
Environment
In a natural water environment, the ambient current and the surface heat
exchange are usually not constant but varies with time. Also, the rate of
excess heat discharge may vary with time. Thus, the steady state problem
formulated in the previous section should be generalized to allow for these
variations in time. The general unsteady problem is very complex and
will not be solved here. In this section, a somewhat simpler unsteady
problem will be formulated and solved. In particular, the current and sur-
face heat exchange are allowed to be time varying, but are assumed to be
uniform in the space coordinates, i. e. , no current shear will be considered.
The rate of excess heat discharge is also allowed to be time varying. It
is clear that this problem is substantially more complicated and cumber-
some from a computational point of view. For example, the time history
of variations of the input functions must be specified before the solution
can be obtained. The method of approach to solve this problem is similar
to the development in the previous section and will be summarized below.
5. 4. 1 Formulation
Neglecting longitudinal mixing as before, Eq. (5. 17) becomes:
CC + u - = - -(K -) + - 2 --(K - -) - K c (5. 62)
x ?y y y z z z d
The corresponding boundary conditions are:
K K c at y = 0 (5. 63)
y y e
K = 0 aty=h (5.64)
y y b
Note that u in Eq. (5. 62) and Ke in Eq. (5. 63) are known functions of
time.
1 56
-------
The source condition at x = 0 is taken to be:
c(o,y,z,t) c(o,y,o,o) F(t) exp - j (5.65)
2 2 (y-y) 2 2
1- [ -----
zoL h j
0
5. 4. 2 Method of Moments
Define moments as before (see Eqs. (5. 25) to (5. 27)) and integrate Eq.
(5. 62) with respect to z, we obtain:
° + u—a. = - (K ) - K c (5. 66)
y y y do
2 + = (K ) - K c + 2K c (5. 67)
y y y d2 zo
The governing equation for the lateral spreading is:
c
Z + u Z = — —(K Z + 2K 1 0 Z + 2K (5. 68)
y y y y c y y z
The boundary conditions for the moments are:
0
K =Kc
y y eo
f atyO (5.69)
K—=Kc
y y e2
or
2
K Z atyO
y y
157
-------
and
2
___ C 2
= = = 0 at y = hb 70)
The source conditions expressed in terms of the moments are
( 2 (y-y 0 ) 1 2
c(o, y, t) = F(t) 1 - L h J (5. 71)
2 r 2 YY )i 2
c 2 (o, y, t) = c(o, y, o)F(t) zo - L h ° J J (5. 72)
2 2 Z(y-y 0 )
°, y) = a - [ h ] J (5. 73)
for y 0 and y - hi 2 y y + hi 2
where F (t) is a prescribed function of time and G (0 y , o).
CO ZO Z 0
Note that , h and y are taken to be time independent. The time
variation of the excess waste heat release is given by the function
The time-varying environmental conditions, Ke and u, are represented
schematically as shown in Fig. 5. 10.
rt
Note that the front of the effluent field is at a value of x given by J udt
as shown in Fig. 5. 10 since the longitudinal mixing is neglected here ?
Thus, the extent of the limit of the region of interest x is related to the
limit of the time of interest t by
tt
x = u(t) dt (5.
The region of solution to be covered is therefore
158
-------
x
tt
7zz
tt
Pt
udt
0
t
Figure 5.10
Representative sketches in cases of unsteady
turbulent diffusion.
xt
K(t)
e
U(t)
tt
tt
F (t)
Co
t
159
-------
tt
dt (5.75)
Define a new independent variable by
= x - $udt (5.76)
Equation (5. 62) in ( , y, z, t) variables become:
= —(K - -) + - Kdc (5. 77)
and the equations for the moments, Eqs. (5. 66) and (5. 67), become:
0 0
- Kdc (5. 78)
c 2
= —(K- --—) + ZKc - KdcZ (5. 79)
Equations (5. 77) to (5. 79) are all independent of the new variable .
However, the source condition is dependent upon ; i. e., Eqs. (5. 71)
to (5. 73) apply at = - J :udt• Therefore, the number o independent
variables has not been reduced by the - transformation although the
governing equations appear to be simpler. Along = constant, we
are following a certain part of the effluent field downstream. In particular,
= 0 represents the front of the effluent field, i. e. , following the very
first part of the release. Negative values represent following later
portions of release. Since longitudinal mixing (in x-direction) is neglected,
there is no exchange between adjacent portions in the x-direction. Thus,
we can treat the problem by investigating each portion of the release as it
travels downstream. For the portion released at time t = t , the source
distributions are c 0 (o, y, t.) and c 2 (o, y, t.). By solving Eqs. (5. 78) and
(5. 79) we obtain the solution for this portion of the release. Note that x
is related to the time variable t by:
160
-------
= (5.80)
By solving a number of cases with different t. from 0 to t , the solution
for the whole region of (x, t) is obtained.
5. 4. 3 Dimensionless Equations and Numerical Solutions
Dimensionless variables are defined as in Sec. 5. 3. 4. In addition to those,
we also define:
= t xe/u (5. 81)
The governing equations, (5. 78), (5. 79) and (5. 80), in dimensionless forms
are:
= (K’ - K’dc’ (5. 82)
12 = (K’ ) + 2K’ c’ - K’ c’ (5. 83)
y y y zo do
and
x’ = ju’ dt’ (5. 84)
The boundary conditions are identical to those given by Eqs. (5. 57) and
(5. 58) except K = K’ (t).
e e
Again, the Crank—Nicolson method was employed in solving this problem
numerically. A Fortran IV program entitled ‘IJTD” was prepared and listed
in Appendix D to handle this particular mathematical model.
The unsteady problem (UTD) requires substantially more input data to
specify the problem than the corresponding steady problem (PTD). In
161
-------
particular, it is necessary to specify the functions F(t), u(t) and Ke(t)
all quantities being dimensionless. Moreover, the program requires a
substantially longer time on the computer particularly if the functions
u(t) and Ke(t) are specified for many values of t. Since the general
basic model is the same as that used in PTD, it can be expected that
similar results should obtain. Only a few cases have been run using UTD.
The results are similar to those of PTD. These results are difficult to
present in such a form as to give ready comparisons with those obtained
from PTD, since the solution depends on the previous history of F,
u, and K. Figure 5. 11 shows one such comparison. The solid line in the
figure is c 0 (x, o) from case PTD-CC-lOO. The points are from using UTD
with the following functions for F , u, and K
co e
t u K F
e Co
o i. 0.0 1
0.2 1. 0 1
0.4 1. 0.1 0.5
0.6 2. 0.1 0.5
0.8 1. 0.2 0.8
1.0 0.5 0 1.0
162
-------
X, NORMALIZED HORIZONTAL DISTANCE
Figure 5. 11 Comparison of PTD with UTD to illustrate effect
of unsteady current, rate of heat discharged, and
surface heat exchange coefficient.
-------
5. 5 Summary and Discussions
In this chapter, two mathematical models have been developed for the
calculation of the distribution of excess temperature due to the effects
of ambient turbulence, current, and surface heat exchange. Both models
assume passive diffusion and ignore longitudinal dispersion. The first
model (PTD) treats the case of a steady release into a steady unidirec-
tional shear current while the second model (UTD) treats the case when
the discharge, the ambient uniform current, and the surface heat
exchange coefficient are time varying. Two computer programs (PTD
and UTD) based on these models are listed in Appendices C and D.
With these programs, the excess temperature distribution can be
determined given the input conditions.
It should be pointed out that this model applies only after the initial
phases of mixing (jet mixing and surface spreading) has subsided and
the buoyancy of the discharge no longer influences the dynamics of the
flow. The initial phases of mixing have been treated in earlier chapters
of this report.
It was not possible to perform a detailed parametric study based on the
models within the scope of this investigation. This should be done in the
future.
Since in practical situations, the conditions are usually unsteady, the
second model (UTD) is likely to be more useful. However, in that model
the ambient current is assumed to be unidirectional and uniform whereas
in practice, shear currents are likely to occur. The model should
therefore be extended to include these effects. Such a model can be
developed by using the concept of superposition. Thus, the case of an
instantaneous release into a general environment should first be solved and
the results superposed. This method has the further advantage of in-
corporating the longitudinal diffusion which is ignored in the present models.
164
-------
It is recommended that this general model of unsteady passive
turbulent diffusion in an arbitrary unsteady environment be
developed in the future.
165
-------
CHAPTER 6 APPLICATION OF THE MATHEMATICAL
MODELS TO PRACTICAL PROBLEMS
In this report, several mathematical models and computer programs
have been developed for the prediction of the excess temperature dis-
tribution in a large body of water resulting from the discharge of heat
such as from power generation plants. Each of these models deals with
a specific portion of the mixing phenomenon.
Although all these models are more general than previously existing ones,
they cannot be regarded as complete. Moreover, no unified model is
available which can calculate the excess temperature distribution from
the beginning phase of discharge through the terminal stage of passive
turbulent dispersion. It is the purpose of this chapter to provide a prac-
tical guideline by which the various models developed herein can be used
in succession to arrive at the solution of practical problems.
Before discussing the practical application of these models, they will
first be briefly summarized. In Chapter 3, two mathematical models
have been developed. The first one (RBJ) solves the problem of mixing
involved in a sub-surface discharge of heated water through a multiport
diffuser from the discharge to the point when either the effluent reaches
the surface of when it reaches its terminal level of ascent. This model
is more general than previously available ones in that in includes jet inter-
ference and an arbitrary ambient density gradient. The second model de-
veloped in Chapter 3 deals with the time dependent surface spreading of
the effluent. The primary purpose of that model is to provide time and
length scales of the phenomenon. Application of this model requires several
numerical coefficients which are as yet not available. These should be ob-
tained by experiments in the future. In Chapter 4, the steady state surface
buoyant jet discharged horizontally is analyzed. The two-dimensional case
of a slot jet is analyzed in detail while the axisymmetric case is also in-
vestigated. The more realistic case of a slot jet of finite length is not
solved and must await future studies. However, the two-dimensional case
167
-------
can be used in certain situations such as if the discharge slot is wide.
In any case, the predictions based on the two- dimensional model should
be conservative in the sense that the predicted temperature excess would
be larger than the actual one, since the model does not include lateral nix-
ing. In Chapter 5, again two models and computer programs are developed.
The first one (PTD) examines the steady passive turbulent dispersion in a
non-uniform current while the second one (UTD) examines the unsteady case.
However, while PTD allows the current to vary with depth, UTD assumes
a uniform, though time varying current. Previous models on this phase of
dispersion assumes constant diffusities and uniform and steady conditions.
It is strongly recommended that before these models and computer programs
are used in a practical situation, the individual using them thoroughly under-
stand the assumptions involved in their derivation and their limitations. This
can be achieved by studying the previous chapters of this report. With this
in mind, the following sections of this chapter are prepared to aid in the
practical applications of these models.
6. 1 Subsurface Discharges
In the event the discharge of cooling water is made at depth, the program RBJ
(Chapter 3) should first be used to obtain the buoyant jet portion of the mixing
phenomenon. The reader is referred to Section 3. 2. 2 for a discussion on the
use of this program. This program terminates the calculation either when
the diluted effluent reaches the water surface or when it reaches its terminal
level of ascent. This can be seen from the output of the program. In either
case, the temperature excess, dilution ratio, and jet width at the end of this
phase are available from the program output.
Having obtained these quantities, the program PTD (Chapter 5) should be
used to continue the calculation. Besides the environmental conditions such
as water depth, the diffusities, and current profiles the program PTD also
requires knowledge of the initial conditions of h, y, and L, the source
thickness, source level, and source width respectively. y should be taken
168
-------
as the terminal depth of ascent from RBJ if the pool is submerged or
zero if the pool is surfaced. It should be noted that since, in PTD,
the excess temperature is just being passively carried along, the flux
of excess temperature, which is proportional to the product u L h C
0 0 max
at the source of PTD,is known from the flux of excess temperature
at the end of RBJ. Here u is the ambient current, C the excess
max
temperature. Thus if we take C to be the same as provided by RBJ,
max
then there is freedom in choosing only one of the two quantities L 0 andh 0 .
Two choices are available. First, L may be chosen to be the total
length of the diffuser plus the jet width. Second, h may be chosen to be half
the jet width at the end of RBJ. It is proposed that separate calculations
be made based on these choices and the worse of the two regarded as a
conservative estimate for design purposes. The reader is referred to
section 5.3. 1.2 for a discussion of the environmental conditions. The
program PTD uses dimensionless quantities in order to minimize the
number of necessary inputs. These are defined in section 5.3.4. Also
included in section 5.3.4 are example solutions which should serve as a
guide on the use of the program. The output of the program include
and C which are dimensionless plume width parameter and dimension-
max
less plume width parameter and dimensionless centerline temperature
excess, each normalized with respect to their values at the source center.
From these and the output from RBJ, the actual temperature excess and
plume width can be readily obtained.
6.2 Surface Discharge
In the event the discharge is made at the surface through a relatively wide
discharge structure, the programs SBJ2 coupled with PTD can be used
to provide an estimate of the excess temperature distribution. Since the
program SBJ2 is based on a two-dimensional slot jet, the prediction would
be more accurate the wider the actual discharge structure. In any event,
the effect of lateral spreading (not included in SBJ2) is to widen the plume
169
-------
which would promote a ‘faster rate of dispersion. Thus the use of SBJZ
constitutes a cons ervative approach.
The calculations based on SBJZ should be carried out to the point of the
internal hydraulic jump and from that point, the program PTD may be
used using the conditions from SBJ2 after the jump as the source conditions
for PTD. The source thickness for PTD may be chosen as the depth of
the flowing layer from SBJ2. The source width L can then be obtained by
a balance of the flux of excess temperature. In the event the source is
inundated, the new source conditions can be obtained in the manner as
described in example 2, section 4.4 This can then be used as the source
condition for PTD. The event that SBJZ would predict a jet all the way is
not likely based on typical values of the relevant paramaters. In the un-
likely event it is the case, then PTD is not needed. SBJZ itself would
probably provide a sufficient estimate of the temperature excess.
It should be reiterated that use of SBJ2 (two-dimensional) for a practical
discharge structure (not two-dimensional) would result in overestimates
of excess temperatures, thus leading to conservative designs. The three-
dimensional problem analogous to SBJ2 should be analysed in the future
to obtain a better prediction tool.
An alternative approach to the use of SBJZ for the initial mixing stage,
especially in the case of narrow discharge structures or in the event
SBJZ predicts a jet all the way, is to use simple submerged jet theory
(e.g. Albertson, Et at (1950) ) to calculate the dilution, based on which
the excess temperature can be obtained. The actual temperature should
be between the predictions based on these two alternatives (i.e. SBJZ-
PTD and simple jet theory.)
170
-------
R EFER ENCES
Abraham, G., ‘Jet Diffusion in Stagnant Ambient Fluid”, Deift Hyd. Lab.,
Pub. No. 29 (1963).
Albertson, M. L. , Dai, Y. B. , Jensen, R. A. , and Rouse, H. , “Diffusion
of Submerged Jets”, Trans. ASCE, 115 (1950).
Bowden, K. F., “Turbulence”, Chapter VI of “The Sea” Vol. I edited
by M. N. Hill, Interscience Publishers, New ‘York (1962).
Brooks, N. H. , “Diffusion of Sewage Effluent in an Ocean Current”, Waste
Disposal in the Marine Environment, Pergamon Press, Inc. , New
York, New York.
Brooks, N. H. and Koh, R. C. Y., “Discharge of Sewage Effluent from a
Line Source into a Stratified Ocean”, XI Congress, Int’l Assoc.
for Hydr. Res. (1965).
Defant, A., “Physical Oceanography”, Vol. I, The MacMillan Co., New
York, New York (1961).
Edinger, J. E. and Geyer, J. C. , “Heat Exchange in the Environment”,
Research Project RP-49, Dept. of Sanitary Engineering and Water
Resources, The John Hopkins University, Baltimore, Maryland,
June (1965).
Edinger, J. E. and Polk, E. M., “Initial Mixing of Thermal Discharges
into a Uniform Current’, Report No. 1, National Center for
Research and Training in the Hydraulic Aspects of Water Pollution
Control, Department of Environmental and Water Resources
Engineering, Vanderbilt University, Nashville, Tenn. Oct. (1969).
171
-------
Ellison, T. I-i. and Turner, 3. S., Turbulent Entrainment in Stratified
Flows”, 1. Fluid Mech., Vol. 6, Pt. 3, Oct. (1959).
Fan, L. N., “Turbulent Buoyant Jets into Stratified or Flowing Ambient
Fluids”, Rept. No. KH-R- .15, W. M. Keck Lab. of Hydr. and
Water Resources, Calif. Inst. of Tech. , Pasadena, California
(1967).
Foxworthy, 3. E., “Eddy Diffusivity and the Four-thirds Law in Near-
shore (Coastal Waters)”, Allan Hancock Foundation, Rept. 68-1,
Univ. of Southern California (1968).
Golubeva, V. N., “The Formation of the Temperature Field in a
Stratified Sea”, Bull, of Acad. of Sci. of the USSR, Geophy. Ser.
(Transl. by F. Goodspeed), No. 5, pp. 4670-4671 (1964).
Harremos, P. , “Diffuser Design for Discharge to a Stratified Water”,
The Danish Isotope Center, Copenhagen, Denmark, 18 pages (1967).
Hayashi, T. and Shuto, N., “Diffusion of Warm Water Jets Discharged
Horizontally at the Water Surface”, Proc. XII Cong. , Int’l Assoc.
for Hydr. Res. (1967).
Isayeva, L. S. and Isayev, I. L., “Determination of Vertical Eddy Diffusion
in the Upper Layer of the Black Sea by a Direct Method”, Issue No. 2,
1963 series, Soviet Oceanography Trans. of the Marine Hydro-
physical Institute, Acad. of Sci. of the USSR, (Transi. by Scripta
Technica, Inc. ) pp. 22-24 (1963).
Jen, Y., Weigel, IL, and Mobarek, 3., “Surface Discharge of Horizontal
Warm Water Jet”, Proc. ASCE, 3. Power Div. , Vol. 92 (1966).
172
-------
Keulegan, G. H., “Laminar Flow at the Interface of Two Liquids’,
3. Res. Nat. Bur. Stands., 32 (1944).
Kolesnikov, A. G. , Ivanova, ? . S. , and Boguslavska, S. G. , “The Effect
of Stability on the Intensity of Vertical Transfer in the Atlantic
Ocean”, Okeanologiya, Vol. 1, (4), English Translation (1961).
Kolesnikov, A. G., Panteleyev, N. A. , and Pisarev, V. D., “Results
of Direct Determination of the Intensity of Deep Turbulent Exchange
in the Atlantic”, DokI. Akad. Nauk USSR, 155, No. 4, (Transl. by
Scripta Technica, Inc. ), pp. 3-6 (1964).
Lean, G. H. and Whillock, A. Z., “The Behavior of a Warm Water Layer
Flowing Over Still Water”, Eleventh International Congress, LAHR,
Leningrad (1965).
Lofquist, K. , “Flow and Stress Near an Interface Between Stratified
Liquids”, Physics of Fluids, 3 (1960).
Morton, B. P. , Taylor, G. I. , and Turner, 3. S. , “Turbulent Gravita-
tional Convection From Maintained and Instantaneous Sources”,
Proc. Roy. Soc. London (A), 234 (1956).
Munk, W. H., Ewing, G. C., and Revelle, R. R., “Diffusion in Bikini
Lagoon”, Trans. Amer. Geophy. Union, 30, (1), pp. 59-66,
February (1949).
Okubo, A., “A Review of Theoretical Models for Turbulent Diffusion in
the Sea”, J. of Ocean. Soc. of Japan, 20th Anniv. Vol. (1962).
Orlob, C. T., “Eddy Diffusion in Homogeneous Turbulence”, 3. Hyd.
Div. Proc. ASCE, September (1959).
173
-------
Orlob, G. T. and Selna, L. G., “Temperature Variations in Deep
Reservoirs”, 3. Hyd. Div. Proc. ASCE, Feb. (1970).
Ozmidov, R. V., “Turbulent Exchange in a Stably Stratified Ocean”, Bull.
Acad. of Sci. of the USSR, Atm. and Oceanic Phys. Ser. 1, (Transi.
by D. and V. Barcilon), pp. 493-497 (1965).
Parr, A. E., “On the Probable Relationship Between Vertical Stability
and Lateral Mixing Processes”, J. Conseil, perman. internat.
explorat. mer., II, No. 3 (1936).
Riley, G., “Parameters of Turbulence in the Sea”, 3. Marine Res. , 10,
No. 3 (1951).
Schlichting, H., “Boundary Layer Theory”, McGraw Hill, New York,
New York (1960).
Schuert, E. A., “Turbulent Diffusion in the Intermediate Waters of the
North Pacific Ocean”, J. Geophysical Research, 75 1970).
Sharp, J. J., “Spread of Buoyant Jets at the Free Surface”, I. Hyd. Div.
Proc. ASCE, May and Sept. (1969).
Snyder, W. H., “A Field Test of the Four-Thirds Law of Horizontal
Diffusion in the Ocean”, M. S. Thesis Dissertation, U. S. Naval
Postgrad. School, June (1967).
Stefan, H. and Schiebe, F. R., “Experimental Study of Warm Water Flow
into Impoundments, Part I, II and III”, Rept. 101, IOZ and 103,
St. Anthony Falls Hydr. Lab., Minneapolis, Minn. (1968).
Wada, A., “A Study on Phenomena of Flow and Thermal Diffusion Caused
by Outfall of Cooling Water”, Proc. 10th Conference on Coastal
Engineering, Tokyo, Japan, Vol. II, Sept. (1966).
174
-------
Yih, C. S., ‘Dynamics of Nonhomogeneous Fluids”, The MacMillan Co.
(1965).
175
-------
AFJ ’hINiJ1X A
The problem discussed in Chapter 3 on the dispersion in a row of buoyant
jets can be solved using the program listed in this appendix. The numerical
integration utilizes a fourth order Runge-Kutta scheme. To facilitate the
use of the program, the following lists are prepared relating the names of
variables used in the text to those used in the program.
Input :
NC
DO
UO
TO
DEN1
THE TA 0
DJ
SPA CJ
D(I=l, NC)
TA(I=l, NC)
DENA(L1, NC)
ALPHAB
A LPHAS
LAMBDR
LAMBDS
GRA VA C
B emarks
number of points for specifying
ambient
diameter of individual jets
velocity of jet discharge
temperature of discharge
density of discharge
angle of discharge
depth of discharge
jet spacing
depth at which ambient specified
ambient temperature
ambient density
In Program
In Text
D
0
u
T
p 1
e
0
d
L
T
a
a
S
x
r
xs
g
gravitational acceleration
177
-------
Output:
In Text In Program Remarks
x x
y Y
afib JET WIDTH
Q/Q 1 DILUTION
T JET TEMP
p JET DENSITY
p AMB DEN
a
T AMB TEMP
a
T- T DELTA T
a
178
-------
C PROGRAM RBJ——ROW BLJIJYANT JET IN A STABLY DENSITY—STRATIFIED ,
C STAGNANT ENVIRONMFNT
0001 DIMENSION TA(50) ,D(50),DENA(50) ,ET(50),FD(50),YT(50)
000? DIMENSION Y( ),YP(6)
0003 REAL LAMBDR,IAMBDS,M
0004 COMMON LAMBDR,LAMBDS,M,H,ALPHAR,ALPHAS, NC,ET,ED,PAT,GRAVAC,YT,IK
1, ICHEK, I Q, SPACJ
0005 204 READ (5,1) NC,00,UO,T0,DEN1,THETAO,DJ,SPACJ
0006 1 FORMAT(1I1O,7E10.5)
0007 IF (00) 2,2,3
0008 2 CALL EXIT
0009 3 READ (c,10) (D(T),TA(T),DENA(I),I=l,NC)
0010 10 FORMAT(3F10.5)
0011 READ (5,11) AIPHAR,ALPHAS,LAMBDR,IAMBDS,GRAVAC
0012 11 FORMATtRE IO.5)
0013 PAI 3.14159265
0014 DC 999 1=1,NC
00 1 999 YT(I)=1)J—D(T)
0016 THETA=THETAO*P4I/180.
0017 ICHEK=0
0018 [ =0
C CHECK PHYSICAL UNITS
0019 IF (GRAVAC—900.) 97,97,98
0020 97 IF (GRAVAC—30.) 101,99,99
C IN FPS UNITS
0021 99 WPITF (6,100) D(],Ufl,T0,DEN1,THETAO,DJ,SPACJ
0022 100 FORMAT (75H1POW BUOYANT JETS IN AN ARBITRARILY DENSITY STRATIFIED
1STAGNANT FNVLRONMENT///SX, I3HJET DIAMETER=,1F6.2,4HFFET,5X,
?23HJET DISCHARGE VELOCITY=,1F6.2,8HFFFT/SEC/SX, 2ÔHJET DISCHARGE T
3EMPERATLJRF, 1F6.2,13HDFGREE FAHPEN, 5X, 22HJFT DISCHARGE DENS 1TY,
41F10.7, 11HGRAM PER M1/4X,18H JET DESCH. ANGLE=, 1F6.2,8H DEGREES/
5X,2OHJET DISCHARGF DEPTH=, 1Fo.2,4HFEET, 5X,L6HJET SPACING C—C ,
61F6.2,4HFEFT)
0023 GO TO 110
C IN MKS UNITS
0024 101. WRITE (6,102) PQ,UO, TO, DENI, THETAO,DJ, SPACJ
0025 10? FORMAT (75HIPOW BUOYANT JETS IN AN ARBITRARILY DENSITY STRATIFIED
1STAGNANT ENVIRflNMENT//15X, 13HJET DLAMETER=,1F6.2,6HMETERS,5X,
223HJ T DISCH.4PGE VFLflCITY ,1F6.2,8HMET./SFC/5X, 2oHJET DISCHARGE I
-------
3EMPERATURE=,1F6.2,13HDEGRFF CENTIG, 5X, 22HJET DISCHARGE DENSITY=,
41F1O.7, IIHGRAM PER ML/4X,18H JET DISCH. ANGLE, 1F6.2,8H OEGREES/
.5X,2OHJET DISCH RGF 0EPTH ,1F6.2,6HMETERS, 5X,I6HJET SPACING C—C=,
6 1F6.2, 6HMETFRS)
0026 GO TO 110
C IN CGS UNITS
0027 98 WRITE (6,103) D0,UO, TO, DEN I, THFTAQ,DJ, SPACJ
0028 103 FORMAT (75H1R.OW BUOYANT JETS IN AN ARBITRARILY DENSITY STRATIFIED
ISTAGNANT FNVIRONMENT///5X, 13HJET DIAMETER=,IF&.2,4H CM.,5X,
2 8 3HJET DISCHARGE VELOCITY=,1F6.2,8H CM.ISECI5X, 26HJET DISCHARGE I
3EMPERATURE=,1F6.2,13HOFGREE CENTIG, 5X, 22HJET DISCHARGE DENSITY ,
41F10.7, IIHGRAM PER ML/4X,18H JET DISCH. ANGIF=, 1F6.2,8H DEGRFES/
55X,2OHJFT DiSCHARGE flEPTH ,1F6.2,4I-4 CM., 5X,L6HJET SPACING C—C=,
8 1F6.2,4H CM.)
0029 110 WRITE (6,111)
0 30 111 FQRMAT (///5X,1HX,IOX,IHY,12X,9HJET WIDTH,ÔX,8HDILUTION,6X,8HJET T
IEMP,4X, ILHJFT DFNSITY,6X,RHAMB DEN ,5X, 8HAMB TEMP,4X,7HDELTA 1)
Oflhl S0.
C TO FIND REFERENCE TEMPERATURE AND DENSITY
0032 I 1
0033 IF (OJ—D(IR) 112,113,114
C 0034 113 TR=TA(I )
0035 DENR=OENA(IR)
0036 GO TO 118
0037 112 IR=IR+1
0038 IF (DJ—fl(IR)) 112,113,117
00.39 114 WRITE (6,120)
0040 120 FORMAT(5X,53H INSUFFICIENT DATA ON AMBIENT DENSITY AND TEMPERATURE
1)
0041 GO TO 204
0042 117 SL=(DJ—D(IP))/(D(IP—1)—D(IP))
0043 TR=TA(IR)+SL*1T4(!R—1)—TA(IR))
0044 DENR=DENA(IR)+SL#(DENA( IR—1)—DENAUR))
C INITIAL CONDITIONS
0045 118 Y(1)=PAI*DO*DO*UO*0.5
0046 M=Y( I )*(JQ*O.5
0047 VOLFJ=Y(1)
0048 H=M*COS(THETA)
0049 Y(2)=M*SIN(THETA)
-------
0050 Y(3)=Y(1)*(DENR—DEN I)/DENR*0.5
0051 Y(4)=Y(l)*(TR TO)/TR*0.5
oc. y(5)=6.?*oo*cos(THE TA)
0053 Y(6) =6.2*00*SIN(THETA)
0054 IQ=0
0055 IP=0
0056 IK=2
0057 SQLAM=(l.+LAMBDR*LAMBDR)/(IAMBDR*IAMBDR)
0058 SQRLAM=SORT( 1.+LAMBDS*LAMBDS)/LAMBDS
C CALCULATION OF DENSITY AND TEMPERATURE GRADIENTS
0059 NC I=NC—1
0060 00 912 I=1,NCI
0061 I1 1+1
0062 DP I=YT(T1)—YT(1)
0063 ET(I) (TA(I1)—TA(I))/(TR*DP1)
0064 912 ED(I) (DENA(I1)—DENA(T))/(DENR*DP1)
C CHOICE CF INTEGRATION STEP
0065 OS1=0 0120.
0066 0S2=DJ/2000.
0067 K=1
0068 IF (DSI—DS2) 301,301,302
0069 301 DS DS1
0070 GO TO 303
0071 302 DS=DS2
C INTEGRATION BY RUNGE—KUTTA METHOD
0072 K=1
0073 303 CALL RUNGS (S,DS,6,Y,YP,L)
0074 304 Y20=Y(2)
0075 CALL RUNGS (S,DS,6,Y,YP,L)
0076 IF (Y(2)*Y20) 20,21,21
0077 20 K=K+1
0078 IF (K—3) 21,22,22
0079 22 IF (ICHFK—1) 204,511,204
0080 2l CONTINUE
C LOOP FflP TRANSITION POINT TWO
0081 IF (ICHEK—2) 5l3,514,204
0082 513 IF (ICHEK—1) 203,206,206
0083 203 TRANW=SPACJ
C ROUND JET SOLUTION
-------
0084 514 IF (Y(6)—OJ) !530,531,53L
0085 531 WRITF (6,532)
0086 532 FORMAT (1OX,2OHTHIS IS FREE SURFACE)
00R7 GO TO 204
0088 530 IF (10) 533,533,206
0089 533 M S0RT(H*H+Y(2)*Y(2))
0090 WIDTH2.*Y(1)/SQRT(PAI*M)
0091 IF (WIDTH—TRANW) 207,206,206
C PRINT SPACING CONTROL
0092 207 SJP=2.*OO
0093 PI=IP*SJP
0094 IF (S—PT) 220,221,221
0095 220 GO TO 304
0096 221 IP=IP+1
009w DENDIF=SQLAM*DFNR*Y(3)/Y(1)
0098 TDIF=SQLAM*TR*Y(4)/Y(1)
0099 DILU=Y(1)IVDIFJ
0100 IF (DENDIF) 401,920,920
0101 sf01 DENDIF=OENDIF*0.5
0102 TDIF=0.5*TD.IF
C TO FIND AMBIENT DENSITY AND TEMPERATURE VALUES
0103 920 IY=2
0106 906 IF (Y(6)—Y1(TY)) 9o0, o1,9O2
0105 90! DENttA=DFNA(IY)
0106 TAA=TA(iY)
0107 IY=IY+1
0108 GO TO 909
0109 900 IY=IY—l
0110 IF (V(6)—Y1(IY)) 900,901,905
0111 905 IYYIY+1
0112 SYY fY(6)—YT(IYV))/(YT(IY)YT(IYY))
0113 TAA=SYY*(TA(IY)—TA(IYY))+TA(IYY)
0114 DENAA SYY*(DENA(IY)—DENA(lYY))+DENA(IYY)
0115 GO TO 909
0116 902 IY=IY41
0117 GO TO 906
OIlS 909 TJ=TAA—TDIF
011 DFNJDENAA—DENDIF
0120 TOIFM=—TOIF
-------
0121
0122
0123
0124
01.25
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
01 38
0139
0140
0141
0142
0143
0144
o 145
0146
0147
01.48
0149
0150
0151
0152
0153
0154
WRITE (6,222) Y(5), Y(6),WIDTH, DILU, TJ, DENJ,DENAA, TAA,TD!FM
222 FORMAT (9G14.7)
GO TO 304
C SLOT JET SOLUITON
C CHEC)( TRANSITION POINT ONE OR TWO
206 IF (Y(6)—DJ) 522,511,511
511 ICHEK=ICHEK+1
IF (ICHEK—2) 512,512,204
C TRANSITION POINT TWO
512 S=SO
Y( 1 )=Y1
Y(2)=Y2
Y(3)=Y3
Y(4)=Y4
v(5)=Y5
Y(6)=Y6
IP=IPC
LKIKC
IQ=0
IY=IYC
L=O
K=KI
WRITE (6,520)
520 FORMAT (lOX, 2OHTRANSITION POINT TWO)
GO TO 303
522 IQ=1
IF (ICHFK—1) 240,241,241
C TRANSITION POINT ONE
240 WRITE (6,122?)
1222 FORMAT (lOX, 2OHTRANSITJON POINT ONE)
C STORE SOLUTIflNS AS INrTTAL CONDITIONS FOR TRANSITION POINT TWO
So=s
Y1=Y( 1)
Y2=Y(2)
Y3=v( 3)
Y4=Y(4)
y5& ’(5)
Y=Y(6)
TRANW=2.*ALPHAS*SPACJ/(PAT* LPHAR)
-------
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
01.69
0170
01.71
01.7?
0173
1PC=IP
KI=K
IKC=iK
TCHFK=ICHEK+1.
I YC=IY
C Ph INT SPACING CONTROL
241 PI=IP*SJP
IF (S—PT) 304,501,501
501 IP=IP+1
M=SQRT(H*H+Y(2)*Y(2)
WIDTH=Y(1)*Y(l)/(SQRT (PAI)*M*SPACJ)*2.
DENDIF=SQRLAM*DENR*Y(3)/Y( 1.)
TDIF=SQRIAM*TR*V(4)/V(1)
DILU V(1 )/V(UFJ
IF (DENDTF) 402, 906,906
402 C0NST 0.5*SQRT(PAI*0.5)
DENDTF=CONST*DEND IF
TDIF=CONST*TDI F
GO TO 906
END
-------
0001 SUBROUTINE DERIVE (S,N,Y,YP)
0002 DIMFNSION Y(6), YP(6)
0003 DIMENSION FT(50),FD(50),YT(50)
0004 REAL LAMBDR,LAMf3DS, M
0005 COMMON IAMBDR,LAMBDS,M,H,AIPHAR,ALPHAS, NC,ET,ED,PAI,GRAVAC,VT,IK
1, ICHFK, IQ,SPACJ
C COMPUTATION OF DENSITY AND TEMPERATURE GRADIENTS AT Y
0006 814 IF (Y(6)—YT(l)) 811,811,812
0007 812 IFIV(6)—YT(NC)) 806,813,813
DOOR 811 EDD=FD(1)
0009 ETTET(1)
0010 GO TO 70
0(111 813 EDD=FD(NC—1)
0012 FTT ET(NC—1)
0013 GO TO 70
0014 806 IF (Y(6)—YT(IK)) 800,801,80?
0015 801 EDD (FD(1K)+ED(1K—1))*O.5
0016 ETT=(ET(IK)+FT(IK—1))*0.5
0017 tK [ K+1
00 18 GO TQ 70
0019 800 IK=II<—1
0020 IF(Y(6)—YT(IK)) 800,801,805
0021 805 EDO=FD(IK)
0022 FTT=IT(IK)
0023 IK=IK+1
002’ GO TO 70
0025 802 IK=!K+1
0026 IF (1K—NC) 814,814,807
0027 807 WRITE (6,808)
0028 808 FORMAT(IOX,25H THIS IS THE FREF SURFACE)
0029 RETURN
0030 70 IF (IQ) 71,71,72
C ROUND JET SOLUTION
0031 71 ENTRAN ?.*ALPHAR*SQPT(2.*PA1*M)
0032 CLAM (t.+LAMBDP*LAMBDR)/2.
0033 GO TO 73
C SLOT JET SOLUTION
0034 72 ENTRAN=2.*SQP.T(2.)*ALPHAS*SPACJ*M/Y(1)
0035 CLAM=SQRT( (1 .+LAMBDS*LAMBDS) /2.)
-------
I-
0036 73 SQRflTM=SQRT(Y(2)*Y(2)+H*H)
0037 YP(1)=ENTRM’
0038 YP(2)=CLAM*GRAVAC*Y(1)*Y(3)/SQROTM
0039 YP(3)=Y(I)*F00*Y(2)/SQR OTM
0040 YP(4)=Y(I)*ETT*V(2)/SQR OTM
0041 YP(5)=H/SQR OTM
0042 YP(6)=Y(2)/SQRUTM
0043 RETURN
0044 END
-------
0001 SUBROUTINE RUNGS (X,H,N,Y,YPRIME, INDEX)
000? DIMENSION Y(7) ,YPRIME(7) ,Z(7),W1(7),W2(7hW3(7),W4(7)
CRUNGS — RIJNGF—KUTTA SOLUTION OF SET OF FIRST ORDER O.D.E. FORTRAN U
C DIMENSIONS MUST BE SET FOR FACH PROGRAM
C X INDEPENDENT V PI4BLE
C H INCREMFNT DELTA X, MAY BE CHANGED IN VALUE
C N NUMBER OF EQUATIONS
C V DEPENDENT VARIABLE BLOCK ONE DIMENSIONAL ARRAY
C YPRIME OFTRIVATIVE BLOCK ONE DIMENSIONAL ARRAY
C THE PROGRAMM R MUST SUPPLY INITIAL VALUFS OF V(1) To Y(N)
C INDEX IS A VARIABLE WHICH SHOULD BE SET TO ZERO BEFORE EACH
C INITIAL ENTRY TO THE SUBROUTINE, I.E., TQ SOLVE A DIFFERENT
C SET OF EQUATIONS OR TO START WITH NEW INITIAL CONDITIONS.
C THE PROGRAMMER MUST WRITE A SUBROUTINE CALLED DERIVE WHICH COM—
C PUTES THE DERIVATIVES AND STORES THEM
C THE ARGUMENT LIST IS SUBROUTINE DEPIVF(X,N,V,YPRIME)
0003 IF (INDEX) c,5, l
0004 1 DO 2 I=1,N
0005 W1( I )=H*YPRI ME(I)
0006 ?Z(I) . Y(T)+(W1(T)*.5)
0007 A=X+H/2.
— 0008 CALL DERIVE(A,N,L,YPRIME)
0009 DO 3 I=1,M
001.0 W2(I)=H*YPRIME(L)
0011 3 Z(I)=V(I)+.5 W?(!)
001? A=X+H/2.
0013 CALL DEPIVE(A,N,Z,YPRIME)
0014 DO 4 I=1,N
W3(1)=H*YPRIME II)
0016 4Z(I)=V(T)+W3(I)
0017 A=X+H
0018 CALL OFRIVE (A,N,Z,YPRIME
0019 DO 7 I=1,N
0020 W4(I)H YPRIME(I)
0021 7 V(I)=Y(1)+(((2.*(W2(I)+W3(I)))+WIU)+W4(I))/6.)
0022 X=X+H
00.23 CALL DERIVE (X,N,Y,YPRIME)
0O21e GE) TO 6
0025 5 CALL DERIVE (X,N,Y,YPRIME)
0026 INDEX=1
0027 6 RETURN
0028 END
-------
APPENDIX B
The problem discussed in Chapter 4 on the two-dimensional surface buoyant
jet can be solved using the program listed in this appendix. The reader
should fully understand the investigation reported in Chapter 4 (Sec. 4. 2)
before attempting to use this program. To facilitate the use of the program,
the following lists are prepared relating the names of variables used in the
text to those used in the program.
Input:
In Text In Program SBJZ Remarks
e E entrainment coefficient
0
k CAY dimensionless surface heat
exchange coefficient
F F2 densimetric Froude number
0
1/B EPS inverse of Reynolds number
XSTOP value of x to stop integration
D control variable for step size
I ] II It I I II
Output :
e E
0
F F2
0
x dimensionless distance
T T dimensionless density deficiency
u U dimensionless velocity
h H dimensionless thickness
2
F = F -p---- FR local densimetric Froude number
oTh
h 2 HZ layer thickness after jump
F FB2 Froude number after jump
2
189
-------
C
C PROGRAM SBJ2—TWO DIMENSIONAL SURFACE HORIFONTAI BUOYANT JET
C ESP IS 1/RE
C
0001 DIMENSION Y(3),YP(3)
p002 COMMON F ,CAY,F2,EPS, XSTOP,DX,OXP
0001 1 READI5,10,END=999) E,CAV,F2,EPS,XSTOP,D1,t)
0004 10 FORMAT(7F 10.6)
00 5 WRITF(6,i000) E,CAY, 2,FPS
00(16 1000 FOKMAT(1HI,1X,2HE,F 10.6,2HK:,F1O.6,3HF2,F10.3,4HEPS,F10.6,f//,
17X,1HX,13X, IHT,13X,LMU,13X, IHH, 13X,2HFR,12X,2HH2, 12X,3HFR2)
00”7 XPRINT=O.0
0008 IF(F2.LE.1.) XSTflP= 8./CAY
0009 IF(c-2.LE.1.) 01=0.1
00.10 !F(F2.IE.1.) 0=1./ 50./CAY
0011 OX=D1*D
0012 DXP=OX—0.00001
0013 X=0.0
0014 Y(1)=1.c
0015
0016 Y(3)=1.0
0017 N=3
0018 1=0
Øfl 9 CALL RtJNGS(X,DX,N,V,YP,L)
0?0 3 CALL RUNGS(X,DX,N,Y,YP,L)
0021 IF(X.LF.XPP!NT) GO 10 3
00?? IF(X.GE.10.*D) DX=O.2*0
0023 IF(X.GF.40.*D) DX=0.5* [ )
0024 Lr(X.GE.100.*D) DX=D
00?5 IF(X.GF.200.*D) DX=2.*t)
00’ IF(X.GE.500.*D) ox=5.*r)
0027 IF(X.GE.1000.*O) DX=10.*D
00 )8 IF(X.GE.D) DXP=10.*OX—O.00001
0029 FR=Y(2)*Y(2)/Y( l)/V(3)*F2
0030 IF(FR.LF.0.) GO TO 1
0011 H2H I=(SQRT(1.+R.*FR)—1.)/?.
0012 FR2=FRIH2H I**3
0O’ 3 i12=H2H 1*Y(3)
0014 XX=X÷0. 0000I
-------
nU: 5 WRITE(6,100)XX,Y(1),V(2),Y(3),FR,H2,FR2
0036 100 FORMAT(IH ,7G14.7)
00 7 IF(F2.LT.l.) GO TO 11
0038 IF(FR.LF.1.) GU TO 1
0019 11 IF(X.GE.XSTOP) GO TO 1
0040 XPP1NT XP 1NT+DXP
0041 GO TO 3
004? 9 CALL FXIT
0043 END
-------
0001 SU WOUTINE RUNGS (X,H,N,Y,YPRIME,INOEX)
0002 DIMENSION Y(3) ,YPRIME(3),Z(3),Wt(3),W2(3),W3(3),W4(3)
CRUNGS — RUNGE—KUTTA SOLUTION OF SET OF FIRST ORDER O.D.E. FORTRAN 11
C DIMENSIONS MUST BE SET FOR EACH PROGRAM
C X INDEPENDENT VARIABLE
C H INCREMENT DELTA X, MAY BE CHANGED IN VALUE
C N NUMBER OF EQUATIONS
C V DEPENDENT VARIABLE BLOCK ONE DIMENSIONAL ARRAY
C YPRIME DERIVATIVE BLOCK ONE DIMENSIONAL ARRAY
C THE PROGRAMMER MUST SUPPLY INITIAL VALUES OF Y(1) TO Y(N)
C INDEX IS A VARIABLE WHICH SHOULD BE SET TO ER0 BEFORE EACH
C INITIAL ENTRY TO THE SUBROUTINE, I.E., TO SOLVE A DIFFERENT
C SET OF EQUATIONS OR TO START WITH NEW INITIAL CONDITIONS.
C THE PROGRAMMFR MUST WRITE A SUBROUTINE CALLED DERiVE WHICH COM—
C PUTES THE DERIVATIVES AND STORES THEM
C THE ARGUMENT LIST IS SUBROUTINE DERIVEIX,N,Y,YPRIME)
0003 IF (INDEX) 5,5,1
0004 1 00 2 I=1,N
0005 W1(I)=H*YPR!ME(I)
0006 2 t(I)=Y(I)+(W1(I)*.5)
0007 A=X+H/2.
0008 CALL DERIVE(A,N,L,YPRIME)
0009 00 3 I=l,N
0010 W2(I)=H*YDPIME(I)
0011 3 L(1)=Y(I)+.5*W2(I)
001? A=X+H/2.
0013 CALL DEPIVE(A,N,Z,YPRIME)
0014 00 4 I=1,N
0015 W3( T )=H*YPRIME( I)
0016 4 L(I)Y(I)+W3(I)
0017 A=X+H
0018 CALL DERIVE (A,N,Z,YPRIME)
0019 00 7 I.=1,N
0020 W4(I )=H*YPRIME( I)
0021 7 Y(I)=Y(I )+I( (2.*(W2(I)+W3( I) ))+Wl( I)+W4( I) )/6.)
0022 X=X+H
0023 CALL DERIVF (X,N,Y,YPRIME)
0024 GO TO 6
0025 5 CALL DERIVE (X,N,V,YPRIME)
0026 INOEX I
0027 6 RETURN
0028 END
-------
0001 SUBROUTINE DERTVE(X,N,V,YP)
0002 DIMENSION Y(3),YP(3)
000 COMMON E,CAY,F2,EPS,XSTOP,DX,DXP
0004 EE=E*FUN(Y(1)*Y(3)/Y(2)IY(2)/F2)
0005 YP(1)=—CAY*Y(I)/Y(2)/Y(3)—EE*Y(1)/Y(31
0006 VP(3)=(—2.*EF—Y(1)*Y(3)*YP(1)/2./F2/Y(2)/Y(2)—EPS/Y(3)/Y(2))/(Y(I)
I*Y(3) /F?/Y(2 )/V(2)—l.)
0007 YP(fl Y(2)/Y(3) *(EE—YP(3)
ooo RETURN
0009 END
0001 FUNCTION FUN(X)
0002 IF(X.GF.0.R5 FUN 0.
0003 IF(X.t.T.0.85) FUN (2./(1.+X/0.8 )I.)**I.lS
0004 RETURN
0005 END
-------
APPENDIX C
The problem discussed in Chapter 5 on the passive turbulent dispersion
from a steady source can be solved using the program listed in this appendix.
The numerical scheme used is based on the Crank-Nicolson method. To
facilitate the use of the program, the following lists are prepared relating
the names of variables used in the text to those used in the program.
Input:
In Text In Program I emarks
NEND A program control number. If
not equal to zero program will
continue, otherwise the program
wilt exit.
NDY number of variations of y-mesh
s c heme s
NDX number of variations of x-step
sizes
NEXP number if variables < 10 NEXP,
it will take it zero
X LAMBDA dimensionless quantity
y YO II
0
h HO ii
0
u UFS
0
y YE1
el
y YE 1
e
YK1
1
YKZ
‘K2
K3 YK3
YK 4
4
R BETAI
F- ,’
BETAZ 1
z
K CKE
CKD
195
-------
Input (continued):
In Text In Program Bemarks
XDY(I) x to change y-mesh scheme
NDYT(I) number of y-mesh changes
NPR(I) number of printout of the
y-mesh lines
DY(I, J) mesh size in y constant for
NYC grids
NYC(I, 3) number of grids that y has
mesh size DY
DX(I) step size in x constant for
NXC steps
NXC(I) number of steps that x has
step size DX
NPX number of x steps for one
printout
Output :
See Input List LAMBDA See Input List
II YO
1 HO IT
YK 1
YKZ
YK3
YK4
BE TA 1
BETAZ I’
UFS
YE1
YE
K KE dimensionless quantity
a SIGMAZ
z
C CMAX ‘I
max
196
-------
C PR1]G AM PIE) —PASSIVE TURBULENT DIFFUSION OF A CONTINUOUS SOURCE
0001 INTEGER 01
0002 REAL LAMBDA
0OO DIMENSION A(100),B(100),C(l00),D(100),F(100),Y(100),rM(100,3),
ISOL(100,3),V0(99),Si(°9),CMAXt99), 0Y19,9), DX(9),NYC(9,9),
2NXC ( 9) , NTAB( 9) , Xl(9) , X2 (9) ,C0(99) ,U(99
0004 DIMENSION NDYT(9) ,YY(100) ,NPX(9) ,XDY(9) ,NPR(9)
0005 COMMON LAMBDA,YU,Hfl,YF1,YF,YK1 ,Y 2,YK3,YK4,BETA1,BETA2,CKE,CKD,
I V, CM, UFS
0006 COMMON /HOID/A,B,C,D,F, SOL,YO, S7,CMAX,CD
0007 01=6
0008 100 READ (5,1)NENO,NOY,N0X,NEXP
0009 1 FcJRMAT(1315)
0010 PA!=3. 1415927
0011 TESTXP=1./I0. *NEXP
0012 IF (NEND) 3,4,3
0013 4 CALL EXIT
0014 3 RFAD(5,?) IAMBDA,YO,HO,UFS, VE1,YE,VK1,VK2,VK3,YK4,BETA1,BETA2,
1 CKE,CKD
0015 2 FORMAT (8E10.5)
0016 READ (5,834) (XDY(I),NDYT(T),NPR(I), I=1,NDY)
0017 834 FORMAT (8(F5.t,13,I?))
0018 83? DO 24 1=1,9
0019 DO 847 J=1,9
0020 DY(I,J)=0.0
0021 847 NYC(I,J)=0
0022 X1(I)=O.0
0023 X2(1)=0.0
0024 DX(I)=0.0
0025 NXC(I)=0.O
0026 24 NTAb(I)=0
0027 8 WRITE (01,6) LAMBE)A,YO,H0,YKI,YK2,YK3,YK4, BETA1,BETA2,UFS,
lYE 1,YE,CKE,CKF)
00?8 6 FORMAT(9H1LAMBOA =,G12.5,//7H V I ] =,G1.2.5,15X,5H 1-40 =,G12.5,//9H
l Y—PROF.,2X,4HVKI,G12.5,2X,4HVK2,G12.5,4HY 3,G12.5,2X,4HYK4 ,
?G1?.5,2X,6HBETA 1,G1?.5,2X,6H8FT42,G12.5//9H U—PROF. ,2X,SHUFS =,
3G12.5,2X,4HYE I,G12. 5,2X,3HYE=,G12.5,//17H SURFACE EXCHANGF,5X,
45H KF =,G12.5,//13H DECAY COEFF.,1OX,SH KU =,012.5/)
0029 N OY1=NDY—l
-------
0010 00 1.10 I 1,MDY1.
0011. NSAV=NDYT(I)
0032 READ (5,5) (DV(J,j), NVC(J,I),J 1,NSAV)
0033 WRITE (OT,IL I) I, XOV(I)
0034 WRITE (01,7) (DY(J,I), NYC(J,I),J 1,NSAV)
0 03! 111 FORMAT (!10,5X,6H X =,Gt2.5)
0036 110 CONTINUF
0037 READ (5,55) (DX(I),NXC(I),NPX(I), 11,NDX)
0038 55 FORMAT (4(F1.0.5,215))
0039 5 FORMAT (4(E10.5, 15,5X))
0040 WRITE (OT,7) (DX(j),NXC(I),I=1,NOX)
0041 7 FORMAT (5(G12.5,15))
C SET UP TABLE
0042 EKF=CKE* DY(1,1)
0043 Y(1)—DV(1,1)
0044 Y(2)=0.
0045 (=2
0046 NSA V=NDYT(1)
0047 DO 10 I=1,NSAV
0048 OELY=OY(I,1I
0049 NUM=NYC(I,1)
0050 00 15 J=1,NUM
0051 K=K+1
0052 YIK)=Y(K— l)+ DELY
0053 15 CONTINUE
0054 10 NTAB(I)=K
0055 NPRTNT=NPRU)+1
0056 NTA B(NSAV)=O
0057 M=K+1
0058 Y(M)=Y(K)+DEIY
0059 M1=M—1
0060 M2=M1—1
0061 WRITE (01,8) (Y(I),I=1,M)
0062 R FORMAT (/1 H TABLE OF INITIAL Y/(8G 13.5))
0063 WRITE (01,9) (NTAB(I),I=t,NSAV)
0064 9 FORMAT (/?H NTAB =,516)
C SET UP SOURCE CONDITIONS
0065 00 11 I=2,M1
0066 AY=ABS(Y(I)—YO)
-------
0067 ALPHAJ=0. SsHfl
00 P 00 13 J=1,3
0069 13 CM(I,J).=0.
0070 F(J—1)=CA Y(I—1,!)
0071 U(I—1).= US U B(y(I))
0072 IF (AY—ALPHAI) 12,11,11
0073 12 SAVI=SQRT(1.—(AYIALpHA1)**2)
0074 CM(I,1)=SQRT(s4v1)
0075 CM( I ,2)=SAVL*CM( 1,1)
0076 CM( I ,3)=CM(I ,2)
0077 00 11 J=1,3
0078 IF (ABS(CM(I,J))—TFSTXp) 872,87?,11
0079 872 CM(I,J)=o.O
0 R0 11 CONTINUE
C SET BOUNDARY CONDITION AT X=0
0031 DO 14 J=1,3
0082
0083 14 CM(M,J)=CM(M2,J)
C LOOP ON NUMBER OF DELTA X—S
0084 X=O.
0085 IDY=2
C PRINT SOURCE CONDITIONS
0086 L=1
0087 00 17 I=2,NPRINT
0038 IF (ABS(CM(I,1))—1.OF--O8) 17,17,16
0039 16C0(L)=CM(I,1)
0090 Y0(L)=Y(j)
0091 SZ(L)=SQRT(CM(r,2)/ CM(!, 1))
009? CMAX(L)=CM(J ,1)/5i(L)
0093 L=L+1
0094 17 CONTINUE
0095 L=L—l
0096 WRITF (OT,90) X
0097 90 FORMAT(1H//3H X,Gl2.5,//,5X,lHy,11x,?Hco, IIX,7HSTGMA Z,7X,4HCMAX,
1/)
0098 WRITE (IJT,491) (YO(1),CO(j),SZ(J),CM AX(I), 11,L)
0099 491 FORMAT (4G13.5)
0100 D I ) 50 NDXL=1,NDX
0101 DELX= DX(NDXL)
-------
NPRINT=NPR( IDY)+1
NSAV=NDYT( iDY)
DO 123 I=1,NSAV
DELY=flV(I ,IDY)
NUMM=NYC(1,!DY )
00 124 J=1,NUMM
K=K+1
YY(K)=YV(K—1)+OELY
124 CONTINUE
123 NT4B(I)=K
NT4B(NSAV)=0
MSAV=M
M=K+l
YY(M)=YY(K)+DEIY
M1=M—1
M2=M1—1
C SET UP PROPER SOLUTION AT DISTANCE N*DX
C LOOP ON YY
IXK=2
DO 126 I=2,M
IK= IKK
DO 127 J=IK,MSAV
I KK=J
IF (ABS(YY(I)—Y(J))—O.00001) 129,129,511
IF (YY(1)—Y(J)) 128,129,127
DO 131 IJ=1,3
SOL( T,IJ)=CM(J,IJ)
GO To 126
128 DO 132 1J=1,3
13? SOL(!,IJ)=(CM(J,IJ)—CM(J—1,IJ) )*(YY(I)—Y(J—1) )/(Y(J)—Y(J—1) )+CM(J—
it, IJ
GO TO 126
NUM=N C (NDXL)
IF (ABS(X—XDY(IDY))—0.00001)
C SET UP NEW Y TABLE
121 VY(1)=—DY(1,IDY)
YY(2)=0.
EKE=CKE*DY(1 ,IDY)
K= 2
121,121,122
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
N) 0118
0119
0120
01 ?1
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
01 34
0135
511
1 29
131
0136
-------
0117
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
o i si
0152
0153
0154
— 01 55
0 1 56
0157
015
0159
0160
0161
0162
01 63
0164
0165
0166
0167
0168
0169
0170
0171
127 CONTINUE
126 CONTiNUE
C RESET 6OUNDARY CONDITIONS
DO 130 1J 1,3
SOL(M,IJ)=SOL(M2,fJ)
130 SQL( 1,IJ)=S0L(1,IJ)—2.*FK [ *SQL(2,iJ)
DO 13 1=l,M
Y( I )=YY( I)
DO 133 IJ 1,3
133 CM( I,1JI=SOL(I,IJ)
00 378 L=2,M1
tJ( I—I )=USUB(Y( I))
378 F(I—1)=CAy( I—i ,!)
IDY IflY+1
C SET UP MATRIX COEFFICIENTS FOR A CONSTANT DELTA X
122 1=1
IDYY WV—i
X1(L)=0.5* OEL X/Dy(1,ID YY)**2
00 20 I=2,M1
J=I—1
IF (I—NTAE3(L)) 21,22,21
21 A(J)=—X1(L)*F(J)
B(J)=X1(L)*(f-(I)+F(J))+tJ(J) I-CKO*DELX
C(J)=—X1 (L)*F( I)
GO TO 20
22 X2(L)=DELX/(flY(L,IDyy)*rjy(Lf , yy)*( 0 y( , IOYY)4-IJV(L÷1,Ioyy)))
+CKD*DELX
C(J) =—X?(L)*F(1)*Dy(L,Ir)yy)
L=L+1
X1(L)=0.5*DFLX/py(L,jDyy)**2
20 CONTINUE
C SET UP BOUNDARY CONDITIONS
8(1) =8(1) —A( 1) *EKE *2.
C(1)=C(1)+4( 1)
A(M?)=A( M2)+C(M2)
C TR1AN( )L4TF 4TRIX
DO 30 J=2,M2
-------
0172 I=M?—J+1
0173 (I)=8(I)—C( 1)*A(I+1)
0174 30A(I)=A(I)/B(I)
C LOOP TN X—COORDINATF
0175 00 51 NTIME1,NUM
0176 X=X+DELX
C LOOP OVFR NUMBER OF FQUATIONS
0177 DO 52 NFQ=1 ,3
C GENERATE NON—HOMOGENEOUS TERMS
0178 1=1
0179 DO 40 i=2,M1
0180 J=1—1
0181 11=1+1
0182 IF (I—NTAB(L)) 41,42,41
0183 41 OtJ)=CM(1,NFO)*U(J)+X1(U*(F(I)*(CM(I1,NFQ)—CM(I,NEO) )—F(J)*(CM( I,
1NEQ)—CM(J,NFO)
0184 GO TO 43
0185 42 D(J)=CM(I,NEO)*U(J)+X2(L)*(DV(1,TDYY)*F( I)*(CM(T1,NEO)—CM(I,NEQ))—
1DY(L+l,IDYV)*FIJ)*(CM(I,NEQ)—CM(J,NEQ)))
0186 L=L+1
0187 43 CONTINUE
0188 GO TO (40,71,72),NEQ
0189 71 D(J)=D(J)+DEIX*(CAYZ(E,CM)*(SOL(I,1)+CM(I,l)))
0190 GO TO 40
0191 72 D(J)=D(J)+DELX*(CAYZ.(i,SO1)*(SOL(I, 1)+CM(I,1)))
0192 40 CONTINUE
0193 D(M2)=D(M2)/B(M2)
0194 00 66 J=7,M2
0195 I=M2—J+1
0196 66 D(I)=(D(T)—C(I )*D(I+1))/B(I)
C COMPUTE SOLUTION VFCTOR
0197 SOL(2,NFQ)=D(1)
0198 DO 67 I=2,M2
0199 IF (ABS(SOL(J,NEQ))—TESTXP) 881,881,67
0200 881 SOL(1,NEQ)=0.
0?01 67 SOL(I+1,NFQ)=D(I)—A(IF*S1)1(I,NFO)
0202 88? SO1(M,NFQ)=SOL(M2,NFQ)
0203 883 CONTINUE
0204 SOL(1,NFQ)=SOL(3,NFQ)—2.*EKF*SO1(2,NEQ)
-------
0205
0206
0207
0208
0209
0210
0211
01?
0213
0214
021 5
0216
0217
021 P
0219
N) 0220
o 0221
L ) 0222
0223
0224
0225
O??6
0227
0.28
0229
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239
0240
52 CONTINUE
00 73 J1,M
73 SOL(J,2)=0.5*(S OL(J,2)÷SOL(J,3)
IF (MOD(NTIrIE,NPX(NE)XL))) 98,54,98
C COMPUTE INTFGRAL OF CO OVER DEPTH
54 N1=3
SUM=0.0
DO 220 I=1,NSAV
N2=NTAB( 1)—i
IF (t—NS4V) 221,222,221
22? N2=M2
221
DO 2.3 J=N1,N2
223 SUM1.=SUMI+SoL(J,1 )*U(J—1)
SUM=SIJM1*OY( 1,Ir)yv)#SUM
220 N1=N2+2
C COMPUTE DESIRFO OUTPUT
L=1
DO 80 I=2,NPRINT
IF (ABS(SflL(I,1))—1.i )E—o5) 80,80,81
81 YO(L) Y(I)
CO(1)=SCL(T , 1)
SZ2=SOL( I ,2)/S01 (1,1)
IF (SZ2) 82,M3,83
82 SZ (L)=—SQRT(—s72)
CMAX(L ) C0(L)/S7 (1)
GO TO 76
83 SZ(L) SQRT(572)
CMAX( L)=C0(L)ISZ(L)
76 L=Lf1
80 CONTINUE
WRITE (01,290) X,StJM
290 FOPMAT(IHI//7H X = ,GI?.5,5X,8H!(CO*U)=,G12.5//SX, lHy,1LX,3H CO,
l8X,7HSTGM4 1, 7X,4HCMAX,/)
LM=L—l
IF (LM—38) 94,94,95
94 L1=1
12=L—1
GO TO 96
-------
0241 95 L1=l
0242 12=38
0243 96 WRITE (OT,91)(Y0(I),C0(I),S?(I),CMAX(T), 1=11,12)
0244 91 FORMAT (4G13.5)
0245 IF (12—IM) 97,9 ,98
0246 97 11=39
0247 1 2 =1—1
0248 WRITE (01,99)
0249 99 FORMAT (tHu/I)
0250 GU TO 96
C SHIFT SOLUTION 10 NEXT X—STEP
0251 98 DO 92 J=1,2
0252 DO 2 I=I,M
0253 92 CM(!,J)=SO1(I,J)
0254 DO 49 I=1,M
0?55 49 CM(I,3)=CM(I,2)
0256 51 CONTINUE
0257 50 CONTINUE
0258 GO TO 100
0259 END
0001 FUNCTION CAY(I,J)
C COMPUTE KY — VERTICAL DIFFUSION COFFFICIENT
C UNIFORM CAY(=BETA2) IF YK4 IS NEGATIVE
0002 REAL LA 8DA
0003 DIMENSION V(100),CM(100,3)
0004 COMMON LAMBDA,Yfl,H0,YF1,YF,YK I,YK2,YK3,YK4,BETAI,BETA2,CKE,CKD,
1Y,CM ,UFS
0005 P=(V( I )+Y(J) )/2.
0006 IF (YK4) 9,9,1
0007 1 IF (P—YK I) 2,2,3
0008 2 CAV=1.
0009 RETURN
0010 3 IF(P—YK2) 4,5,5
0011 4 CAY=(BETA 1*(YK1—P)+(P—YK2))/(YK 1—VK2)
0012 RETURN
0013 IF (P—YK3) 6,7,7
0014 6 CAY BETA1
0015 RETURN
0016 7 IF (P—YK4) 8,9,9
0017 8 CAY=(BETA2*(YK3—P)+BETAI*(P—YK4))I(YK3—VK4)
0018 RETURN
0019 9 CAY=BETA2
0020 RETURN
0021 END
-------
0001 FUNCTII)N (JSUB(V)
0002 PEAL LAMBDA
0003 DIMENSION Z(I00),CM(100,3)
0004 COMMON LAMBDA,Yfl,HO, YE1,VF,YK1,YK2,YK3,YK4,BFTA1,BETA2,CKE,CKD,
I ,CM,UFS
0005 IF (UFS) 10,10,11
0006 10 USUB 0.
0007 RETURN
0 00R 11 IF (V—YE) 6,6,7
0009 6 IF (Y—YF 1) 5,4,4
0010 5 USUB=UFS
0011 RETURN
0012 4 USUB=UFS*((Y—YE)/(YFI—YE))
0013 RETURN
0014 7 USUB=0.
0015 RETURN
0016 END
0001 FUNCTION CAVZ(L,CM)
C COMPUTE K1
0002 REAL LAMBDA
0003 COMMON LAMBDA,Yfl,Hfl,YE 1 ,VF,YK1,YK2,YK3,yK4,BFTA1,BETA2,CKE,CKD,
IV
0004 DIMENSION Y(I00),CM(100,3)
0005 IF (ABS(CM(1,1))—1.OE—08) 1,1,2
0006 1 CAYZO.
0007 RETURN
0008 ? CAYZ=LAMBDA*(ABS(CM(L,2)ICM(L,1)))**O.666666667
0009 RETURN
0010 END
-------
APPENDIX D
The problem discussed in Chapter 5 on the unsteady dispersion from
a continuous source can be solved using the program listed in this
appendix. The numerical scheme used is based on the Crank .-Nicolson
method.
Input:
In Text In Program Bemarks
NEND see list in Appendix C
N DY It H U
NDX number of variation of t-step
sizes
NEXP see list in Appendix C
NT number of times in specifying
u, K , F values
e co
NXPR number of variations of DXPR
sizes
LAMBDA See list in Appendix C
YO
HO
YK1 It
YK2
yK3
yK4
BETA 1
BE TA 2
CKD
XDY(I) t to change y-mesh scheme
NDYT(I) see list in Appendix C
NPR(I)
TI(I) times at which u, K, Fco
spe cified
207
-------
In Text In Program Remarks
u U(I) u at time TI(I)
K CKE(I) K at time TI(I)
e e
F (t) FC(I) F at time TI(I)
CO CO
DXPR(I) spacing in x for each printout,
constant for NXPRC steps
NXPRC(I) see above
DY(I, 3) y-mesh constant for NYC grids
NYC(I, 3) number of grids that y has
mesh size DY
DX(I) step size in t constant for
NXC steps
NXC(I) number of steps that x has
step size DX
Output :
TI See Input List
U
CKE
FC
LAMBDA
YO
HO
YK1
YKZ
YK3
YK4
BETA 1
BE TA 2
KD
XPR x values for printout
TPR t values for printout
Y see output list for Appendix C
CO II
SIGMZ
CMAX
208
-------
C PROGRAM UTi)—--tJNSTEADY TURBULENT DIFFUSION OF A CONTINUOUS SOURCE
0001 INTFGFR Or
0002 REAL LAMBDA
0003 DIMENSION A(lOO),B(I00),C(IOO),D(100),F(100),Y(lO0),CM(j0O,3),
lSOL(IQO,3),Y0(°9),S7(99),CMAX(gc), DV(9,9), DX(9),NYC(9,9),
2NXC (9) ,N14f3( 9) ,Xl(9) ,X2 ( 9) ,c0( 99) ,U(99)
0004 DIMENSION T(99) ,TI(99),CKF(99),FC(99),XI(99),DxPR(9),NXPRC(9),
LTPR( 99) , XPR( 99 ), SOLU( 100 ,3)
0005 DIMENSION NDYT(9),YY(l00),NPX(9),XDy(9),NPR(9)
0006 CflMMON IAMRDA,Y1],HO, YK1,YK7,YK3,yK4,BFTAI,BET A?, CKD,
1 Y,CM
0007 CCMMI1N /HC)LD/A,B,C,D,F, SOL,Y0,Si,CMAX,C0
0008 flT=6
000 ’ ) 100 READ (5,1)NEND,N)y,NDX,NEXP,NT,NxPP
0010 1 FORMAT(10 15)
0011 PAI 3.14l5927
0012 TESTXP 1Jl0.**N [ XP
0013 IF (NEND) 3,4,3
0014 4 CALL EXIT
001% 3 RFAD(%,2) LAM8DA,Yfl,HO, VK1,VK2,YK3,VK4,BEIA1,BETA2,
ICKO, TE
00 16 2 FORMAT (RF1O.5)
0017 REM) (5,834) (XDY(L),NDYT(I),NPR(1), I=1,NDY)
0018 834F ORMAT(8(F5.l,r3,I2))
COiO READ (5,835) (TI (I) ,U(I),CKF (I),FC(I),I=l,Nt)
0020 835 FURMAT(4E10.5)
0021 TF=TI(NT)
0022 XI(l)=0.
0073 JS=2
0074 DO 730 I 2,NT
0025 XI (I )=XI (1—1 )+0. 5*( U( I) +(J( I—I)) *(T 1(I)—TI (I—I ))
0026 730 CONTINUE
0027 WRITE (6,336) ( T,TI( I),U(I ),CKE(I),FC( I) ,XI( I),I=1,NT)
0028 336 FORMAT (//12H TI,U,CKE,FC/(15,5G13.5))
0029 READ (5,31) (OXPR(I),NXPRC(IhI=1,NXPR)
0030 31 FORMAT (4(FI0.5,I10))
0 ” 31 K=2
003? XPR(1)=0.
0033 DO 32 I=1,NXPR
-------
0034 NXPRCC=NXPRC(I)
0035 OXPRC=DXPR(I)
0036 00 33 J=1,NXPRCC
0037 XPR(K)=XPR(K—1)+DxpRC
0038 IF (XPR(K)—XI(NT)) 331,331,332
0039 331 K=K+1
040 33 CONTINUE
0041 32 CONTINUE
0042 332 NXPR1=K—1
0043 832 DO 24 1=1,9
0044 DO 847 J=1,9
0045 OY(I,J)=0.O
0046 847 NYC(I,J)=0
0047 X1(I)=0.0
0048 X2(I)=0.0
0049 DX(I)=O.0
0050 NXC(I)=0.0
0051 24 NTAB(I)=0
0052 58 WRITE (OT,6) LAMBDA,Yfl,H0,YKL,YK2,YK3,YK4,RETA I,BFTA2,
ICKD
0053 6 FQRMAT(9H ILAMBDA =,G12.5,//7H VI) =,(;12.s,15x,5H HO =,G12.5,//9H
1KY—PROF.,2X,4HYK 1=,G 1?.5,2X,4HYK2=,G12.5,4HYK3=,G12.5,2X,4HYK4=,
2G12.5,2X,6HBFTA1,G12.5,2X,6HBETA7 ,G12.5//13H DECAY COEFF., 1OX,SH
3 KD =,G12.51)
0054 NDYL=NDY—1
0055 DO 110 I=l,NDY1
0056 NSAV=NDYT(1)
0057 READ (5,5) (DV(J,1 ), NYC(J,I),J=1,NSAV)
0058 WRITE (flT,111) I, XDV(1)
0059 WRITF (OT,7) (DY(J,1), NYC(J,I),J=1,NSAV)
o oo O Ill FORMAT (I10,SX,bH X ,G12.5)
0061 110 CONTINUE
0062 READ (5,55) (DX(1),NXC(f),I=1,NDX)
0063 55 FORMAT (4(E10.5,110))
0064 5 FORMAT (4(Fl0.5,15,SX))
0065 WRITF ( (11,7) (OX(I),NXC(I),I=l,NDX)
0066 7 FORMAT(5(G1?.5,I10))
C MAIN LOOP OVER TIME OF RELEASES
0067 NT I=NT—1
-------
0068 00 10? NTI=l,NT 1
0069 XINTI=X1(NT I)
0070 DC) 740 1=NTE,NT
0071 740 X1(I)=XI(I)—XINTI
0072 TPP(1)=O.
0073 J=JS
0074 NXPRL1=NXPR I+1
0075 DC 38 I=2,NXPR11
0076 37fF (XPR(I)—XI(J)) 34,36,35
0077 35 J=J+I
0078 IF(J—NT) 37,37,39
0079 36 TPR(fl=TI(J)
OOP O GO TO 38
0081 39 TPR(1)=—1.2345
0082 NXPR2=I—l
0083 GO TO 1002
0084 34 TPR( I )=Tj(J—l)+(Tj(J)—TT (J—]j )*(XpR( I)—XI(J—1))/(XI(J)—XI(J—l)
1—TIC JS—1
0085 38 CONTINUE
00 6 1002 WRITE (6,333) (I,XPR(I),TPR(I),1=1,NXPR2)
0087 333 FORMAT (8H1XPR,TPR/(15,2G15.5))
0088 MT=2
0089 TIN=TI(NTI.)
0C 0 TF 1=TF—TIN
0091 IF (TFL) 100,100,101
0092 101 FCI=FC(NTI)
0093 RKE=CKF(NTI)
C SET UP TABLE
0094 EKER E*DY( 1,1)
0095 V(1) —DY(1,1)
0096 Y(2)=0.
0097 K=2
0098 NSAV=NDYT(l)
0099 DO 10 I=1,NSAV
0100 DFLY=DY(I,I)
0101 NUM=NVC(I,1)
010? D C ) 15 J=1,NUM
0103 K=K+1
0104 Y(K)=V(K—1)+DELY
-------
0105
0106
0107
0108
O1OQ
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
N) 0120
— 0121
0122
01.? 3
0124
0125
01 6
0127
0128
0 1 2
0130
0131
013?
01 3
0134
0135
0136
0137
0138
01 3
15 CONTINUE
10 NTAB(I)=K
t’IPRINT=NPR( 1 )sj
NTAB(NSAV)=()
M=K4 -1
Y(M)=Y(t +flELY
M1=M—1
M2=M1—1
WRITE (01,8) (Y(T),j=1,M)
A FORMAT (/ 19H TABLE OF INITIAL Y/(8G13.5))
WRITE (flT,g) (NTAB(j),I=1,NSAV)
9 FORMAT (/7H NTA =,516)
C SET tIP SOURC CONDITIONS
00 11 I=?,M1
AY=ABS(Y( I )—Y0)
ALPHA 1=0 • 5*HQ
DO 13 J=1,3
13 CM(I,J)=0.
F (I—I) CAY( I—i, I
IF (AY—ALPHAI) 12,11,11
1? SAV1=SQRT(1.—(AY/AIPHAI)**2)
CM(1,1)=SQRT(SAV I)*FCl
CM(I ,2)=SAV I*CM( 1,1)
CM(I ,3)=CM( 1,2)
00 11 J=1,3
IF (ABS(CM(I,J))—TESTxP) 872,872,11
872 CM(i,J)=0.0
ii CONTINUE
C SET BOUNDARY CONI)ITION AT X=O
00 14 J=1,3
CM(1,J)=CM(3,J)—2.*FKF*CM(2,J)
14 CM(M,J)=CM(M?,J)
C LOOP ON NUMBER OF OELTA X—S
x=O.
IDY=2
PRINT SOURCF CONDITIONS
L=1
00 17 I=?,NPPINT
IF (ABS(CM(I ,l) )—1.OE--OR) 17,17,16
-------
C)140 16 CO(L)=CM(I,1)
0141 Y0(L)=Y(I)
014? SZ(U=SQRT(CM(I,2)/CM(I,t))
0143 CMAX( L)=CM( 1,1) /Sl(L)
0144 1=1+1
0145 17 CONTINUE
0146 LL—1
0147 WRITE (01,90) TIN
0148 90 FORM4T(1H// I5HIRELEASE TIME =,G12.5//,5X,1HY,1IX,2HCO,11X,7HSIGMA
IZ,7X,4HCMAX,/)
0149 WRITE (P1,491) (Y0( I) ,C0( I) ,SZi! ) ,CMAX( I), I=1,L)
0150 491 FORMAT (4G13.5)
0151 DO 50 NDXL=l,NDX
0152 DELX= DX(NDXL)
0153 NUM=NXC(NDXL)
0154 IF (Af3S(X-XDV(TOY))—0.0000l) 121,121,1??
C SET UP NEW V TAbLF
0155 121 YY(I)=—DY(1,IDY)
0156 YY(2)=0.
0157 EKE=RKF*DY(1,I OY)
0158 K2
0159 NPRINT=NPR(IDY)i-l
0160 NSAVNDYT(1DY)
0161 DO 123 I=1,NSAV
0162 DELY=DV(I,IOY)
0163 NUMM=NYC(1,IDY)
0164 DO 124 J=1,NUMM
0165
0166 YV(K)=YV(K—1)+DELY
0167 124 CONTINUF
0168 123 NTM3(I) K
01.69 NTA (NSAV)=0
0170 MSAV=
0171 M=K4-1
0172 VY(M)=YY(K)+PE L.Y
0173 M1=M—1
0174 M2=M1—l
C SF1 UP PROPER SOLUTION AT DISTANCF N*DX
C LOOP ON YY
-------
0,75
0176
0177
0178
0179
0180
0181 511
0182 129
0183 131
01R4
0185
0186
0187
01’38
0189
N 0190
Z 0191
0192
fl193
0194
0195
0196
0197
0198
0199
0.00
0201
0702
0203
0204
0205
0206
0407
0208
0209
0210
IKK=2
00 126 !=2,M
IK= IKK
00 127 J=IK,MSAV
IKKJ
IF (ABS(YY(I)—y(J))—O.000 01) 129,129,511
IF (YY(j)—y(J)) 128,129,127
DO 131 IJ=1,3
SOL I I,IJ)=CM4J, I J)
GO 10 126
1?8 00 132 1J1,3
13? SOIl 1,lJ).=(CM(J,jJ)—CM(J—1,1J) )*(YY(I)—Y(J—l))/(y(J)_y(J-.l) )+CM(J—
ii, IJ)
GO TO 126
127 CONTINUE
126 CONTINUE
C RESET BOUNDARY CONDITIONS
00 130 TJ l,3
SOL (M,I J)=SOL(M2, IJ)
130 SOL(1 ,IJ)=S01( , IJ)—2.*FKE*SQL(2,jJ)
00 133 I=1,M
Yl I )=YY( I)
00 133 IJ=1,3
133 CM( 1,1 J)=SOL (I, I J)
DO 378 I=7,M1
378 Fl I—I) =CAY( 1—1,1)
IDY=I OY+1
C SET UP MATRIX COEFFICIFNtS FOR A CONSTANT DELTA X
1.? 1=1
I DYY= IDY—1
X1(t)=O.5*DFIX/Dy(1,jQyy)**2
DO 20 1=2,MI
J=I—1
IF (I—NTAS(L)) 21,22,21
21 A(J) —X1(L)*F(J)
B(J)=X1(L)*(F(!)+F(J))+I.0O +CKD*DELX
C(J)=—Xl(L)*r( I)
GO If) 20
2? X2(L)=DELX/(py(L,IDyy)*fly(L+1,IDyy)*(Dy(L, JDVY)÷DY(L+1,jDyy)))
-------
0211 A(J)=—x2(L)*E(J)*DY(L+l, IDYY)
0212 B(J) X2(L)*(1)Y(1,IDYY)*F( I)+DY(1+l, IDYY)*F(J))+1.00 +CKO*DEIX
0213 C(J)=—X2(L)*F-(I)*DY(1,IDYY)
0214 1=1.41
0215 Xl(L)=O.5*DFLX/ OY(L,!DYY)**2
0216 20 CONTINUE
0217 C( 1 )=C( 1 )+A( 1)
02l 61=8(1)
0219 Al=4(l)
0220 A(M2)=A(M2)+C (M2)
C TRIANGULATE MATRIX
0221 A(M2) 4( M.)/B( M2)
0222 M3=! 2—1
0223 00 30 J=2,M3
0224 I=M2—J+1
0225 B(I)=B(I )—C(1)*A(I#1)
0226 304(1)=4(I)/B(I)
C LOOP IN X—C1)ORDTNATE
LI 0227 00 51 NTIME=1,NUM
0228 XX+DELX
0229 IF(X—TFI) 1001,1901,102
0230 1001 XT=X+TIN
0231 IR=l
0232 109 IF (XT—TI(TR)) 106,108,107
0233 107 IR=IR+1
0?34 IF (JR—NT) 109,109,735
0235 735 IR=IR—1
0236 108 RKE=CKF(IR)
0237 GO TO 103
0?38 106 IF (IR—l) 733,733,732
0239 733 RKE=CKE(1)
0240 GO 10 103
0241 732 RKE=CKF ( IR)—(TI ( IR)—XT) / (TI ( JR )—T! ( JR—i) )*(CKE( IR)—CKE( IR—1) )
024? 103 EKE=P KE*OY(1,IDYY)
C SET UP BOUNDARY CONDITIONS
0243 Btl)=B1—A1*FKE ‘p2.
0244 P( 1).=B(1)—C(1 )*A(2)
0?45 A(1)=A1/B(1)
C LOOP OVFR NUMBER OF EQUATIONS
-------
024f DO 52 NEQ=1,3
C GENERATE NON—HOMf)GENEOIJS TERMS
0247 1=1
0248 DO 40 1=?,Ml
0249 J=T—l
0750 11=1+1
0251 IF (1—NTAB(L)) 41,42,41
0252 41 D(J)=CM(I,N O)*I.O0+Xl(L)*( (I)*(cM(Il,NFQ)_CM(1,NFQ) )—F(J)*(CM(I,
INFO)—CM(J,NEQ)
GO TO 43
02 4 ‘+2
IDV(I+1,TDYY)*F(J)*(CM( I,NEQ)_CM(J,NEQ)))
0255 L=1+1
0256 43 CONTINUE
0757 GO IL ] (40,71,77),rn-Q
0258 71
0259 GO Tn 40
0.60 72
0261 40 CONTINUE
0262 D(M2) D(M2)/8(M ?)
0263 00 66 J=2,M2
0764 1M?—J+1
0265 66 fl(I)=(D( I)—C(I)*D(j+1))/B(I,
C COMPUTE SOLUTION VECTOR
Q266 SO1(2,NFQ)=D(I)
0767 DO 67 I 2,M?
IF (ABS(S01(I,NEQ))—TFSTXP) 881,881,67
0269 881 SO1(I,NFQ)=0.
0270 67 SOL(I+1,NEQ)=D(Ip—A(I)*SOL(J,NFQ)
0271 IF (A8S(SO1(M2+1,NEQ))—rESTXp) 882,R82,883
0272 882 SO1(M?I 1,NEQ) o.o
0273 883 CONTINUE
0274 SOL(M,NFQ)=SCJL ( M2 ,NEQ)
0?75 SOL(1,NEQ)=S OL(3,NfEQ) —2.*EKE*S0 [ (2,NFQ,
0276 5? CONTINUE
0277 DO 73 J=1,M
0778 73 SUL(J,2)=0.5*(ScjL(J,2)4.S01 J,3)
077 IF (X—TPP(MT)) 98,301,101
0280 301 COEF=(X—TPR(Mr))fD [ Lx
-------
0291 XPRR=XPR(MT)
028? TPRR=TPP(MT)÷T IN
0293 !4T=MT+1
0294 DO 731 T 1,M
0295 SULU(T,1)=S OL(T,1)—(S(IL(y,1)—CM(I,1))*CQFF-
0286 731 SOLU(1,2)=SOt ( I,2)—(SOL( I,2)—CM(I,2))*COFF
C COMPUTE INTFGRAL OF CO OVER DEPTH
0287 54 N1=3
0239 SUM=0.O
0299 DO 220 I=1,NSAV
0290 N2 NTAB(I)—l
0291 IF (I—NSAV) 221,222,221
0292 222 N2=M2
0293 2?l SUM1=(SOL(N1—1,1 )+SOL(N2+1,l))/2.
0294 DO 223 J=Nl,N2
0295 223 SUM1=SUM 1+S0I(J,1)
0296 SUM StJM1 *1)y( I, IOYY)+SUM
0297 220 N1=N2+2
C COMPUTE DESIPED OUTPUT
02 8 L=1
0299 DO 80 I ?,NP INT
0300 IF- (A9S(SOLU(I,1))—l.0F—08) 80,90,81
0301 81 YC(L) Y(I)
0302 C0(L)=SOL(J(I,1)
0103 SL2=S OLU(T,2)/SO1U(I,1)
0394 IF (SZ2) 82,83,33
0305 82 Sl(L)=—SQRT(—S12)
031)6 CMAX( L)=C0( L)/SZ(L)
0307 GO TO 76
0308 83 SZ(L)=SORT(Si?)
0309 CMAX( 1) =C0(L)/S7(L)
0310 76 L=L+1
0311 80 CONTINUE
0312 WRITF (UT,290) TTN,FCI,TPRR,XPRR,SUM
0313 290 FORMAT(1H1//18H TIME OF RELEASE =,G12.5,SX,20H AMOUNT OF RELEASE =
1,GI?.5,5X,7H TIME =,G12.5,5X,5H X =,G12.5,5X,//8H i(CO) =,G1.2.5//
2 5X,1HY,IIX,3H C0,8X, 7HSIGMA l,7X,4HCMAX,/)
0314 LM=L—1
0315 If- (LM—38) 94,94,95
-------
0316 94 11=1
0117 L21—1
0318 GO TO
0319 95 11=1
0320
0321
03??
0323
0374
0375
0326
0327
0178
0329
0330
0331
033?
03 3 3
0334
0315
0336
0337
0338
0339
96
L2:38
96 WRITE (Ut,91)(YO(I) ,C0( I),SZ(T),CMAXL1 ),I=L1,L2)
91 FORMAT (4G13.5)
IF (L2—LM) 97,98,98
97 L139
L2=L-1
WRITE (01,99)
99 FORMAT (tHu/I)
GO TO 96
C SHIFT SflLUTION T i) NEXT X—SIEP
98 DO 9? J=l,?
00 92 I=1,M
92 Chl( I,J)=SOL(1,J)
00 49 I=i,M
49 CM(t,3)=CM(T,2)
IF(TPR( 1T) ) 102,51,51
51 CONTINUE
50 CONTINUE
102 JS=JS+1
GO TO 100
END
-------
0001 FUNCTION CAY(I,J)
C COMPUTE KY — VERTICAL DIFFUSION COEFFICIENT
C UNIFORM CAY(=BETA2) IF YK4 IS NEGATIVE
000? REAL LAMRDA
0003 D!MFNSION Y(100),CM(100,3)
0004 COMMON LAMPDA,Y0,HO, YK1,YK2,YK3,YK4,BETAI,8FTA2, CKD,
I Y,CM
0005 P=(Y( I )i-Y(J))/2.
0006 IF (YK4) 9,9,1
0007 1 IF (P—YK1) 2,2,1
0008 2 CAY=1.
0009 RETURN
0010 3 IF(P—YK2) 4,5,5
0011 4 CAY=(BFTA 1*(YK1—P)f(P—YK2))/(YKI—YK2)
0012 RETURN
0013 5 IF (P—YK3) 6,7,7
0014 6 CAY= BETAI
oois RETURN
0016 7 IF (P—YK4) ,9,9
0017 8 CAY (BETA2*(YK3—P)+BETA1*(P—YK4))/(YK3—YK4)
0018 RETURN
0019 9 CAY=RETA2
0020 RETURN
0021 END
000! FUNCTION CAYZ(L,CM)
C COMPUTE K?
0002 REAL LAP BDA
0003 COMMON LAMBDA,YO,Hfl, YK I,YV2,YK3,YK4, BFTA1.,BETA2, CKD
0004 DIMENSION CM(100,3)
0005 IF (ABS(CM(L,1))—l.0E08) 1,1,2
0006 1 CAYZ=0.
0007 RETURN
0008 2 CAY7=LAMBDA*(ABS(CM(L,2)/CM(t ,1)))**O.ô 6 b 6 ôôôôl
0009 RETURN
0010 END
-------
BIBLIOGRAPHIC : Tetra Tech, Inc., C. Y. Koh and Loh-Nien Fan, ACCESSION NO.
“Mathematical Models for the Prediction of Temperature Distri-
butions From the Discharge of Heated Water into Large Bodies
of Water”, FWQA Publication No. 16130 DWØ1O/70.
ABSTRACT : Mathematical models for heated water outfalls were developed KEY WORDS:
for three flow regions. Near the source, the subsurface discharge into
a stratified ambient water issuing from a row of buoyant jets was solved Mathematical mode
with the jet interference included in the analysis. The analysis of the Thermal outfall
flow zone close to and at intermediate distances from a surface buoyant
jet was developed for the two-dimensional and axisymmetric cases. Far
away from the source, a passive dispersion model was solved for a two-
dimensional situation taking into consideration the effects of shear
current and vertical changes in diffusivity.
A significant result from the surface buoyant jet analysis is the
ability to predict the onset and location of an internal hydraulic jump.
Prediction can be made simply from the knowledge of the source Froude
number and a dimensionless surface exchange coefficient.
Parametric computer programs of the above models are also developed
as a part of this study.
This report was submitted in fulfillment of Contract No. 14-12-570
under the sponsorship of the Federal Water Quality Administration.
-------
Accession Number 2J Subject Field & Group
Ø23C SELECTED WATER RESOURCES ABSTRACTS
INPUT TRANSACTION FORM
Organization
Tetra Tech, Inc.
Title
6 “Mathematical Models for the Prediction of Temperature Distributions Resulting
From the Discharge of Heated Water into Large Bodies of Water”
1 0 J Author(s)
C. Y. Koh and
Loh-Nien Fan
Project Designation
FWQA C6ntract #14-12-570; DW
Note
22 Citation
FWQA, R & D Report 1613ODW 10/70
23 J Descriptors (Started First)
*Mathematjcal Model, *Outfalls, *Thermal Pollution
25 Identifiers (Starred First)
Thermal outfalls
7 Abst act
L Mathematical models for heated water outfalls were developed for three flow regions.
Near the source, the subsurface discharge into a stratified ambient water issuing from a
row of buoyant jets was solved with the jet interference included in the analysis. The
analysis of the flow zone close to and at intermediate distances from a surface buoyant
jet was developed for the two-dimensional and axisynilietric cases. Far away from the
source, a passive dispersion model was solved for a two-dimensional situation taking into
consideration the effects of shear current and vertical changes in diffusivity..
A significant result from the surface buoyant jet analysis is the ability to predict
the onset and location of an internal hydraulic jump. Prediction can be made simply from
the knowledge of the source Froude number and a dimensionless surface exchange coefficient.
Parametric computer programs of the above models are also developed as a part of this
study.
This report was submitted in fulfillment of Contract No. 14-12-570 under the sponsor-
ship of the Federal Water Quality Administration.
WRtO2 (REV. JULY 1969)
WRS) C
SEND TO WXTER ‘dI’ WTlFlC IN FORMATION CENTER
U.S. DEPARTMENT OF THE INTERIOR
WASHINGTON. D.C. 20240
bs$ractot
hirazi. Mostafa A. I
IflStltUtjotl
EPA/Federal Water Quality Administr i-inn
* GPO: 1969—359339
------- |