EPA-R2-73-133
JANUARY 1973
Environmental Protection Technology Series
A User's Manual for
Three-Dimensional Heated
Surface Discharge Computations
Office of Research and Monitoring
U.S. Environmental Protection Agency
Washington, D.C. 20460
-------
RESEARCH REPORTING SERIES
Research reports of the Office of Research and
Monitoring, Environmental Protection Agency, have
been grouped into five series. These five broad
categories were established to facilitate further
development and application of environmental
technology. Elimination of traditional grouping
was consciously planned to foster technology
transfer and a maximum interface in related
fields. The five series are:
1. Environmental Health Effects Research
2. Environmental Protection Technology
3. Ecological Research
U. Environmental Monitoring
5. Socioeconomic Environmental Studies
This report has been assigned to the ENVIRONMENTAL
PROTECTION TECHNOLOGY series. This series
describes research performed to develop and
demonstrate instrumentation, equipment and
methodology to repair or prevent environmental
degradation from point and non-point sources of
pollution. This work provides the new or improved
technology required for the control and treatment
of pollution sources to meet environmental quality
standards.
-------
EPA-R2-73-133
January 1973
A USER'S MANUAL FOR THREE-DIMENSIONAL
HEATED SURFACE DISCHARGE COMPUTATIONS
By
Keith D. Stolzenbach
E. Eric Adams
and
Donald R. F. Harleman
Project 16130 DJU
Project Officer
Frank H. Rainwater
National Environmental Research Center
Corvallis, Oregon 97330
Prepared for
OFFICE OF RESEARCH AND MONITORING
U.S. ENVIRONMENTAL PROTECTION AGENCY
WASHINGTON, D.C. 20460
For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402
Price $1.25 Domestic Postpaid or $1 GPO Bookstore
-------
EPA Review Notice
This report has been reviewed by the Environmental Protection
Agency and approved for publication. Approval does not signi-
fy that the contents necessarily reflect the views and policies
ot the Environmental Protection Agency, nor does mention of
trade names or commercial products constitute endorsement or
recommendation for use.
ii
-------
ABSTRACT
In February 1971, a report by K. D. Stolzenbach and Donald R. F. Harleman
entitled "An Analytical and Experimental Investigation of Surface Discharges
of Heated Water" was published. (Technical Report No. 135, Ralph M. Parsons
Laboratory for Water Resources and Hydrodynamics, Department of Civil Engineer-
ing, M.I.T.; also published as Water Pollution Control Research Series Report
No. 16130 DJU 02/71 by the Water Quality Office, Environmental Protection Agency,
Washington, D.C.) The report described above presented a literature review of
previous analytical and experimental research on heated surface discharges,
the development of a predictive theory for calculating three-dimensional temp-
erature distributions in the near-field region and verification of the theory
based on laboratory work at M.I.T. and elsewhere. The above report also con-
tained a listing of the computer program as originally developed.
Subsequent work and experience in using the computer program for the
calculation of temperature distributions have resulted in modifications and
improvements in the program. In addition, a number of runs covering a wide
range of input parameters have been carried out. These runs were done with
two objectives: (1) to explore the limits of applicability of the theory and
(2) to prepare design charts showing the important characteristics of the
near field temperature distribution for heated surface discharges. These
charts will enable the designer to make rapid estimates of surface isotherms
and vertical thickness of surface jets.
This report presents a review of the theoretical background for the
three-dimensional temperature prediction model, a detailed discussion of the
revised computer program and a case study illustrating the procedure for optim-
izing the design of a surface discharge channel. The revised computer program
flow chart, program listing and a sample of the input and output data are given
in the appendices.
One difference between the revised computer program presented in this
report and the original computer program of February 1971 deserves further
comment. The original program contained the possibility of considering a sloping
bottom in the receiving water. The sloping bottom, was assumed to extend downward
iii
-------
in a linear slope from the bottom of the discharge channel. In the revised
computer program it is assumed that the bottom of the receiving water does
not interfere with the development of the surface jet. Thus the revised
program corresponds to S = °° in the original program. The reasons for elim-
inating the bottom slope effect in the revised program are two-fold: (1) the
mathematical model, as originally constituted, did not adequately predict the
point of separation of the heated jet from the bottom or lateral spreading when
the jet is in contact with the bottom; (2) from the standpoint of environmental
impact, it may be desirable to accept the depth of the receiving water as a
constraint and to design the surface discharge channel to minimize interference
of the heated jet with the bottom.
It is recognized that almost none of the operating power plants that
employ surface discharge schemes have been designed to minimize bottom impact.
This fact makes it difficult to compare field data from many existing plants
with temperature predictions based on the mathematical model presented in this
report.
Further analytical and experimental studies of bottom interference with
heated surface jets are presently underway in the M.I.T. Parsons Laboratory and
will be the subject of a future technical report.
Inquiries relating to the availability of the program source deck should
be directed to Professor D.R.F. Harleman, Room 48-335, M.I.T., Cambridge,
Massachusetts 02139.
iv
-------
TABLE OF CONTENTS
Page
ABSTRACT iii
TABLE OF CONTENTS v
LIST OF FIGURES vi
I. INTRODUCTION 1
1.1 Characteristics of Heated Discharges 1
1.2 Surface Discharges 3
II. THE MATHEMATICAL MODEL 3
2.1 Basic Assumptions 7
2.2 Discharge Structure 9
2.3 Integration of the Equations 16
2.4 Solution of the Equations 20
III. THE PROGRAM 39
3.1 Dimensionless Equations 39
3.2 The Computational Scheme 43
3.3 Input Formats 45
3.4 Output Format 46
IV APPLICATION 47
4.1 Schematization 47
4.2 Use of the Program Output 56
4.3 Case Study - Heated Surface Discharge into a 57
Receiving Water Body of Finite Depth
ACKNOWLEDGEMENT 65
REFERENCES 67
LIST OF SYMBOLS 69
APPENDIX I, Flow Chart for Solution of [a^] ^L c± 73
APPENDIX II. Program Listing 75
APPENDIX III. Input and Output 97
-------
LIST OF__FIGURE!3 "u
Figure page
1.1 Schematic of Heated Discharge 4
2d Coordinate Definitions 8
2.2 Discharge Structure H
2.3 Calculated Surface Discharge Parameters: 22.
IF = 4.4, A= 0.35, k/u = 4.2 x 10~ , V/u = 0
o o o
2.4 Calculated Isotherms of AT/AT : IF = 4.4, A = 0.35, k/u = 23
o o o
4.2 x 10~5, V/u = 0
2.5 Centerline Temperature Rise AT /AT and Dilution, 0, in 26
the Stable Region for k/u = V/u = 0
° ° (h+r)
~m £3 "X"
2.6 Maximum Vertical Depth of the Jet, —• — 27
h b
for k/u = V/u =0
o o
2.7a Jet Parameters for IF = 20.0, V/u = 0, k/u =0 29
2.7b Jet Parameters for IF = 10.0. V/u = 0, k/u =0 in
o o ' o J
2.7c Jet Parameters for IF = 5.0, V/u = 0, k/u =0 31
2.7d Jet Parameters for IF = 2.0, V/u = 0, k/u =0 •}?
o o o
2.7e Jet Parameters for IF = 1.0, V/u = 0, k/u =0 11
O O o
2.7f Jet Parameters for IF = 20.0, V/u = 0.025, k/u =0 34
2.7g Jet Parameters for IF = 10.0, V/u = 0.05, k/u =0 35
2.7h Jet Parameters for IF = 5.0, V/u = 0.1, k/u =0 36
2.7i Jet Parameters for F = 2.0, V/u = 0.1, k/u =0 77
AQ
4.1 Coefficient of Thermal Expansion for Fresh Water
-| r^ "I
3 = - ^ (°F ) as a Function of Temperature T(°F)
4.2 Discharge Channel Schematization i--,
vi
-------
Figure Page
4.3 Limitations on Maximum Jet Thickness by Bottom Topography 52
4.4 Two Layer Flow in the Discharge Channel 53
4.5 Temperature Distribution Plotting Aid 58
4.6 Schematic of Case Study Problem 59
4.7 Calculated Contours for the Case Study Example 64
vii
-------
I. Introduction
1.1 Characteristics of Heated Discharges
The discharge of heated condenser water into natural bodies of water
can be broadly classified into two groups: surface discharges and submerged
discharges. The latter class includes single port submerged discharges and
multi-port diffusers. From the standpoint of minimizing the environmental
impact of thermal effluents, the surface discharge has a number of advantages
when compared with submerged discharges: (1) By careful design of a surface
discharge it is possible to avoid temperature rises and high velocities along
the bottom of the receiving water. (2) The time of travel of organisms entrained
at the condenser water intake is short. (3) Heat dissipation from the surface
of the receiving water is high because of the tendency of the discharge to
form a stratified surface layer.
A typical open cycle system withdraws water from a natural water body
through an intake structure and passes the flow through the turbine condensers
where it undergoes a temperature rise before being returned to the receiving
water through a discharge structure. In-plant technological considerations dic-
tate condenser temperature rises of from 10-30°F with correspondingly large
cooling water flov/s depending upon the amount of rejected heat.
The use of natural water bodies and coastal waters for disposal of
waste heat must take into account the effect upon the environment of the flows
and temperature rises induced in the receiving water. There is general agree-
ment that water temperature increases which approach the sub-lethal range of
impaired biological activity should be avoided. Practically all regulatory
agencies provide this type of protection through controls on both maximum
temperatures and allowable temperature rises and by designating the size of
mixing zones for various types of receiving waters. In addition to environmental
considerations the intake structure must be located and designed to prevent sig-
nificant recirculation of heated water. Separation of the intake and discharge
or the use of selective withdrawal structures are the most common techniques.
The temperature distribution induced in the receiving water by a heated
discharge is determined by the characteristics of the discharge structure and
by the local ambient heat transfer processes. Close to the point of discharge
the momentum of the discharged water creates jet-like mixing of the heated and
ambient water. Within this "near field" region the temperature and velocity
-------
of the discharge decrease because of dilution by entrained water. The magni
tude and extent of the dilution is determined primarily by the nature of the
initial discharge flow, its submergence, velocity, dimensions and temperature
rise above ambient. Mixing increases with increasing discharge momentum and
decreases with increasing temperature rise. The greater the submergence of
the discharge below the water surface, the lower the temperature rise at the
surface will be after mixing. Mixing may also be affected by the presence of
physical obstructions which tend to block the supply of dilution water. Sur-
face discharges entrain a flow at least equal to the discharge flow in the
near field, often up to twenty times as much.
Beyond the near field mixing region the discharge velocity and turbulence
level are of the order of ambient values. In this "far field" region further
entrainment does not occur and the temperature distribution is determined by
natural turbulent convection and diffusion. Ultimately all of the rejected
heat contained in the discharge passes to the atmosphere through the water sur-
face, a process driven by the elevated surface temperatures. These far field
heat transfers are highly variable, being determined by local water currents,
wind, and meteorological conditions.
Consideration of either environmental impact or recirculation on discharge
design must start with the near field temperature distribution. Far field pro-
cesses are generally an order of magnitude less efficient in reducing temperatures
and the far field temperature distribution tends to be more dependent upon the
total amount of waste heat rather than the discharge design. In contrast, a
wide range of dilutions is achievable in the near field through diverse types
of discharge structures. It has been common practice to analyze near field temp-
eratures by constructing scale models of the discharge structure. There is a
pressing need for analytical models which relate the discharge characteristics
to the flow and temperature distribution in the receiving water. Furthermore,
since the analysis must almost always be performed in advance of actual plant
construction, the analytical model must be totally predictive, containing no
undetermined phenomenalogical coeffients.
This report describes the basis, structure, and use of a predictive model
of the three-dimensional behavior of surface discharges of heated water. Emphasis
is placed upon the assumptions underlying the theory, the scope of its validity,
-------
the nature of its limitations, and the proper application to actual discharges,
For similar treatments of submerged discharges, the reader is referred to
references 1 and 2.
1.2 Surface Discharges
The theory presented herein considers a discharge of heated water from. a.
rectangular open channel at the surface of an ambient body of water of infinite
extent in which a current may be flowing (see Figure 1.1). The three-dimensional
temperature distribution depends upon the mixing between the discharged and
ambient water and upon the rate of heat transfer to the atmosphere at the water
surface.
Experimental investigations of three-dimensional buoyant surface jets
have been reported by Tamai (3), Wiegel (4), Jen (5), Stefan (6) and Hayashi (7).
Buoyant surface discharges are distinguished from non—buoyant turbulent jets
by lateral gravitational spreading and by reduction of vertical entrainment
such as described by Ellison and Turner (8) for the two-dimensional case. The
net result of these two processes is that the velocity and temperature distribu-
tions are much wider than deep with the increased surface area raising the
possibility (as suggested by Hayashi) for significant surface heat loss. Pre-
vious analytical treatments (Hoopes (9), Motz (10) ) of heated surface discharges
have failed to take into account, in a single three-dimensional theory; the
roles of buoyancy, initial channel shape, turbulent entrainment, and surface
heat loss upon the temperature distribution.
In this treatment the discharge is assumed to be a free turbulent jet
with a well defined turbulent region in which velocity and temperature are
related to centerline values by similarity functions. An unsheared core region
is accounted for. Turbulent entrainment is represented by entrainment coefficiencs
as first introduced by Morton, Taylor, and Turner (11) and applied to other
buoyant jet problems by Morton (12) and Fan (1) among others. A major contribu-
tion of this work is the treatment of lateral buoyant spreading by incorporating
an assumed distribution for the lateral velocity into the set of integrated
governing equations. Surface heat loss is assumed to be determined by a single
heat loss coefficient as defined by Edinger (13). In the presence of a cross
current in the receiving water the jet is deflected by entrainment of ambient
lateral momentum.
The basic philosophy behind the formulation of the theory is that the
-------
Surface Heat Loss
m
Heated Discharge
Side
' r i SUr
Discharge
'.'. hanne I
Ambient Cross Flow
BMH
Turbulent
Ent ra inment
Turbulent
Entrai nment
Turbulent
Entroinment
Fig. l-l Schematic of Heat
-------
solution should match known non-buoyant jet behavior as the buoyancy terrnp
go to zero. The effects of buoyancy are obtained by use of the basic eq.Jatio^r.
and some judicious assumptions about the jet structure. In this way, the theor
contains no undetermined coefficients and comparison of the theory with observa
tions involves no curve fitting.
Because of the dependency of the treatment upon jet theory, the predic-
tions of the model are valid only in the region where turbulence resulting from
the discharge dominates ambient turbulent processes, i.e. in the near field.
-------
II. The Mathematical Model
2.1 Basic Assumptions
The model considers a discharge Q of heated water at temperature T
and density p from a rectangular open channel of depth h , width 2b , and
initial angle 6 at the surface of a receiving body of water at temperature
T , density p and of large extent laterally and longitudinally. It is
assumed that the bottom of the receiving water does not interfere with the
vertical development of the surface jet. A non-uniform current V may be present
in the receiving water. It is assumed that this current is parallel to the
shoreline; however, its magnitude may vary in the offshore direction as shown
in Figure 2.1. Far from the jet the water surface, n> is uniformly at z = 0.
The flow in the receiving water is characterized by its velocity, density,
pressure and temperature. These four variables are related by equations expressing
conservation of mass, momentum and heat and an equation of state. The solution
of the equations must be developed by setting certain terms in the basic equations
equal to zero on the basis of the following assumptions:
a) Steady flow: -~- = 0
dt
b) Large Reynolds number: viscous terms negligible
c) Boussinesq approximation: density gradients only important in
pressure terms
d) Hydrostatic pressure: p = -| pgdz
) No jet induced motion at large depths: — ^- = —*- = 0
as z -»- -o>
33 9
f) Boundary layer flow: -— « — - and — -
oX dy oZ
Pa~po
g) Small density differences: - «1
pa
h) Mild jet curvature: V/u <1
With the above assumptions the basic equations of mass, momentum and heat con
servation may be simplified to the following form:
Mass Conservation
in + iz jw (2>1)
3x 9y 3z
-------
Cross Flow = V = V(x)
oo
ATr
Top
View
Fig. 2.1 Coordinate Definitions
-------
x-Momentum Conservation
= _
P
J
3x 3y 3z p& J 3x 3y gz
z
y-Momentum Conservation
wv
3x 3y 3z J 3y
Z
Heat Conservation
f JT 2 ^9 (2>3)
J 3y gx
3uT ivT wT 3v'Tl 3w'T
'l
3x 3y 3z 3y 3z
where :
x,y,z = coordinate directions relative to the jet centerline
(see Figure 2.1)
u,v,w = mean velocity components
turbulent fluctuating velocity components
T = mean temperature
1' = turbulent fluctuating temperature
6 = centerline deflection angle (see Figure 2.1)
g = gravitational acceleration
'''
g = — -^- where p = density
p 31
2.2 Discharge Structure
The governing equations (2.1-2.4) assumed for heated surface discharges
may not be solved without further manipulation since the turbulent transfer
terms are not determined. The technique used to develop the solution is to
assume a structure for the velocity and temperature within the discharge, and
boundary conditions at the outer edges, leaving as unknowns only certain values
such as the centerline velocity and temperature. The governing equations may
then be integrated over a cross-section perpendicular to the discharge centerline
This procedure eliminates the unknown turbulent terms and yields a set of first
order differential equations which may be solved for the variables describing
the discharge behavior.
9
-------
The assumed structure of the discharge is shown in Figure 2.2. The
longitudinal velocity and temperature distributions are taken to be as
follows where n is the water surface elevation and u and AT = T - T are
c c c a
the centerline velocity and temperature rise above ambient at z = r\, y = 0:
u = u + VcosO
c
AT = AT
u=uf(?)+ Vcosf
c z
AT = AT t(c )
u = u f(£ ) + VcosG
AT = AT
0< y
-------
Velocity Temperature
Distribution Distribution
LL
Top
V iew
Region
4 \i Region 2
vb\!
Cross Section
Fig. 2.2 Discharge Structure
-------
The above functions were proposed by Abramovich (14). They have the desirable
properties that a distinct jet boundary is defined at t. = 1.0 and that the
velocity and temperature distributions are not identical but are related by
1/2
t = f as implied by Taylor's vorticity transfer theory.
The y momentum equation (2.3) expresses a balance between the lateral
density gradients and lateral convective motions. It may be shown that for
buoyant jets these lateral movements and thus the pressure gradient may be of
the same order as convective transport in the x direction. The balance between
lateral turbulent fluctuating pressures and Reynolds stresses, which is found
in non-buoyant jets, is of second order and is neglected here. Since there are
finite (but second order) lateral velocities in a non-buoyant jet, the lateral
velocity v in equation (2.3) must be intrepeted as the buoyant spreading velocity
in excess of the non-buoyant value. To enable an integration of equation (2.3)
over the jet cross section, the distribution of the lateral spreading velocity,
v, must be specified. The lateral velocity is geometrically related to the longi-
tudinal velocity by:
— = tan cf> (2.7)
where is the lateral jet streamline angle from the centerline in excess of
the non-buoyant value. A distribution for u has already been given; it thus
remains to choose a reasonable distribution for $ such that the following condi-
tions are satisfied: (1) tan <|> = 0 at y = 0 and (2) tan (j> = ( 4^- - e) at y = b+s
where e is the lateral spreading rate of a non-buoyant jet under the same conditions
(cross flow, channel size, etc.). Since the gravitational spread is induced by
r, rn
the lateral temperature gradient, the y dependence of — is used to distribute
tan between y = 0 and |y| = b + s. The result is:
W - 4- fdb N 1/2 | ,
v ~ I ^x ~ e) u?y s<|y|
-------
The above distribution for v insures that for a non-buoyant jet in which
-j^- = e the gravitational spreading velocity is identically zero everywhere.
Depending upon whether r or s are non-zero at a given centerline dis-
tance, the jet cross section will have 1, 2, or 4 regions (see Figure 2.2)
on each side of the centerline (y = 0). To permit integration of the governing
equations over each region separately, the velocities and the turbulent transfe
of heat and momentum are specified at the boundaries of the regions. At the
centerline of the jet, symmetry implies that there is no net transfer of mass,
momentum or heat:
v = u'v' = v'T' = 0
y = 0
(2.9)
At the boundaries between the regions the velocities are assumed to be:
u'w' = 0
w = w
w = w f (£ )
w = 0
u'v' = 0
V = +V
- s
v = +vbf(
v = 0
0<|y|
-------
u .Hi + v .£11 = w
Sx 3y
u'w' = 0
z =
(2.11)
The transfer of heat through the water surface is assumed to be proportional to
the surface temperature rise above ambient:
w'T' = k(T - T )
3
z = n
(2.12)
The coefficient of heat loss, k, has units of velocity and is thus a kinematic
quantity. It is related to the surface heat exchange coefficient, K, defined
by Edinger (13) by:
k =
pc
(2.13)
where p = density and c = specific heat of water. The determination of k(or K)
for a particular case is discussed in a later section.
The outer boundaries of the jet is where entrainment of ambient water
occurs but across which no heat is transferred. The boundary conditions are:
u'w' = w'T' = 0
W = W - VCOSS -r—
e ax
w = wef (cy) - Vcose
u'v' = v'T' = 0
0<|y
-------
The velocities, w£ and v , are manifestations of the entrainment of ambient
fluid into the turbulent region. In plane and axisymmetric non-buoyant jets,
it is known that the entrainment velocity is proportional to the local center-
line velocity. Note that Equation 2.8 gives v=0at |y| =s+b where the
boundary condition 2.14 specifies |v = v . This is because the entrainment
1 e
velocity v is an order of magnitude less than the spreading velocity v and its
contribution to the integration of equation (2.3) is assumed to be balanced
by the turbulent terms, which were neglected in the y equation. Ellison and
Turner (8) have demonstrated that the vertical entrainment is a function of
the gross Richardson number, _ c , in two-dimensional jets. In this study
2
u
c
the entrainment velocities are assumed to be given by
v
e
= a
f BgA
rc
L u
BgAT h
^= a, exp <-C ^f (2.15)
c
Ellison and Turner's data indicate a value of C = 5.0 as appropriate. The
entrainment coefficients, a and a > are to be determined such that the solution
for the non-buoyant case (T = T ) agrees with the experimental observations
O 3.
that the growth of a non-buoyant turbulent region is symmetrical:
db dh
dx dx
ds _ dr
dx dx
(2.16)
For non-buoyant jets discharging into a quiescent receiving water the spreading
rate e , is constant. In these cases, Abramovich gives e = .22 for the
» o o
similarity functions, f and t, used here.
15
-------
The boundary condition at x = 0 is related to the discharge channel
geometry, the flow rate Q , and the initial discharge temperature TQ.
r = h
s = b
h = b = 0
u = u
AT = AT
c o
x = y
2h b
o o
+Vcos6
0
x = 0
(2.17)
2.3 Integration of the Equations
With the velocity and temperature distributions and boundary conditions
stated in the previous section, the equations of motion may be integrated over
the y-z cross section of the jet. The x momentum and mass conservation equations
are integrated over each of the four possible jet regions on each side of the
centerline plane. (Integration over both sides of the jet results in redundant
equations because of the assumed symmetry.) This yields eight equations. The
y-momentum and heat equations are integrated over the entire half-jet cross-
section, yielding two more equations. With this choice of integrating limits,
the terms in the integrated y-equation represent a balance between the lateral
gravitational force and the lateral spreading of the discharge. Another equa-
tion is generated from the y-equation by integrating it over the entire cross
section of the jet on both sides of the centerline. In this case the lateral
spreading terms all drop out, being anti-symmetrical, and the remaining terms
give the rate of deflection of the jet in the presence of a cross current. The
equation set, completed by a simple geometrical relationship between the coor-
dinates (x,y) referred to the jet centerline and the fixed centerline coordinates
(x,y) is given in Table 2.1. Further details in the derivation of these equa-
tions are given in Stolzenbach and Harleman (15).
16
-------
Region 1; continuity rs - — [u + VcosOl + rv - sw =0
dx c s r
Region
2: continuity s -—- [h(u 3^ + Vcos6)J + (11 + VcosO) ^ + w - a ul + v,hl = 0
3: continuity r \-~- [b(u I, + VcosO)] + (u + Vcos0) -^--v +aul-wbI1=0
|_dx xcl J x c dx s ycj r 1
Region 4:. continuity — hb(u I. + VcosO) + (u I, + VcosO) b 4^ + h 4^ + (w, - a u )I,b
dx|_cl J c 1 L dx d:!y n sz c 1
Region
- (v,- - a u ) I,h = 0
b y c 1
Region 1: x momentum (u + Vcos6)[2rs — [u + VcosO] + rv - sw ] + 3g s {— (AT ^-)+I r -r— (AT h) ] = 0
C CLX C S JL Q.X C £, 3 G.X. C
ro o O A -y-
~ [h (u^^"I0 + 2Vcos6u I + V cos 6)] + [u + VcosB] -7-
"^ C J- C Q.X
[s [h (%:
QAI n , -j
+ w [u + VcosG] - a u VcosG + Bg [I. —•—— + I,AT h ^L ] + v h (u I +Vcos6I1) = 0
re sz c 4 dx 3 c dx J b c 2 1
dAT h2
Table 2.1 Integrated Equations for Deflected Buoyant Jets
-------
2 22 2 ds
Region 3: x momentum r 1-^- [b(u !„ + 2Vcos9u I, + V cos 6)J + [u + VcosO] -3— -
'J--: c 2 c 1 c dx
v
[u + Vcos6] + a u Vcos6 - w,b [u !„ + Vcos9I,]+ gg •
c ycjncz 1 ^
I, dAT br
3 c
o J
dAT bh
2 I
+ AT ( ~ + I,hr) ^ =
dx c 2 3 dxj
j 99 999 9 99
Region 4: x momentum 3— [hb(u I. + 2Vcos6u I, + V cos 6)] + [u I. + 2Vcos6u I, + V cos Q]
dx c 2 cl c2 cl
00
b ~ + h ~\ + Iw,b - v,h] [u I0 + VcosOlJ - u I [a b - a h]Vcos6
|_dx dxj h b c 2 1 clsz y
+ eg 1,1
t
^ 2 d 7
•3-4 dx— + Z4ATch dl + X3 ATchb d^ ! - °
Jet y momentum
~
QX
[u 2bl,(r + hi ) + 2Vcos6u bl (r + hi ) + V2cos20b(r+h) ]
Ct) Z C J X
= 0
Table 2.1 (cont'd) Integrated Equations for Deflected Buoyant Jets
-------
Jet Heat
~- u AT (s + bl_) (r + hi.,) + VcosSAT (s + bl.) (r + hl_)
dx |_ c c 7 / c j j J
Jet Bending
+ kAT (s + bI0) = 0
c j
u (s + bI0) (r + hi.) +-2Vcoseu (s + bL.) (r + hi,) +
c / 2 ell
Jet x position
V2cos26(s+b) (r + h)
~ - sine = 0
dx
,] |f - ujsine [-
u Vsine -a (s + bln) + a (r + hi,) = 0
c I sz 1 y .
Jet y position
dx
- cose = o
•""I
\
*3
1
= 5 «•
0
1 _
- J '*«
o
1
= J t(c:
1 2
;) d^ = J (1 - c3/2) d? = .4500
0
^ 3/2 4
;) d? = V (1 - t. ' ) d^ = .3160
o
1
)d£ = J (1 - ?3'2)dc = .6000
1 1
^-J I "
o 5
1
o
1 _
J6 ' I £ (5:
1 1
tc) d^d^ = r r
j -j
o ^
c1'2 dC - 1 a -
o
1
) C1/2d? = J (1
) d^d? = .2143
- .2222
d? = .1333
f(?) t(c)
(1 - ?) d? = .3680
Table 2.1 (cont'd) Integrated Equations for Deflected Buoyant Jets
-------
The values of the non-buoyant (AT = 0) entrainment coefficients which
satisfy equations 2.16 in addition to Table 2.1 are:
ay = - 0
i
a = . 4-2, for s = 0
y 2
a = (I -I_)e for r> 0
z l/o
(2.18)
for r = 0
a =
z 2
where I, and I? are integration constants as given in Table 2.1.
2.4 golution of the Equations
The thirteen equations in Table 2.1 are a first order system of dimensional
differential equations in x for the variables: u , AT , h, b, r, s, Q, x, J, v »
C (— o
v, , w , and w, . The actual solution proceeds first by writing the equations in a
dimensionless form by normalizing each variable by the characteristic values: u ,
AT = T - T , and Ai b . The solution is then determined by the following dimen-
o o a o o
sionless parameters:
u
•^ = — = initial densimetric Froude number
A = h /b = aspect ratio
oo r
k/u = heat loss parameter
V/u = cross flow parameter
The computer program which solves the equations is described in Section IV along
with instructions for its use. The output consists of values of u/u , AT /AT ,
function of
In addition to these variables which describe the jet structure (u , AT , . . ., etc.)
other interesting dependent quantities which are functions of x may be defined. In
the buoyant jet the vertical entrainment is a function of the local densimetric
Froude number (an inverse Richardson number), ]F :
20
-------
3F = (2.19)
/ggAT h
where u , AT , and h are local values at a given distance from the origin.
The total flow in the jet may be determined by integrating the x velocity, u,
over the jet cross section. The ratio of the flow at a given x to the initial
flow is the jet dilution, D, or Q as labeled by the output,
u (r + I,h)(s + I b) + Vcos6(r + h) (s + b)
D - — (u + vcose )h b (2'20)
O O O 0
Similarly the effect of surface heat loss may be evaluated by calculating the
ratio of convected excess heat flow in the jet to the initial excess heat flow,
HT. 4
u AT (r + I,h)(s + I^b) + Vcos6AT (r * h)(s + b)
HT = -£—S Z 1 9 (2.21)
AT (u + Vcose ) h b
O O O O O
Finally a dimensionless time of travel along the jet centerline from th^ enc7
of the discharge channel to x is computed.
(2.22)
The structure of a heated surface discharge is shown for a particular theoretical
calculation in Figures 2.3 and 2.4. The main features are:
1) A core region in which the centerline velocity is constant
and the centerline temperature rise decreases very slightly.
The dilution, D, and the local densimetric Froude number, 1^,
do not vary greatly in this region. The magnitude of FL is
much larger than IF because the initial, depth of the turbulent
region, h, is zero. There is no significant surface heat
loss in the core region.
21
-------
N>
NJ
Core Region
Stable Heat Loss
Entrainment Region Region Region
ATC
ATn
HT
50
100
Fig. 2.3 Calculated Surface Discharge Parameters- Fo =4.4, A = 0.35, k/uo = 4.2 x I 0~5, V/uo =0
-------
0.9 OB 0.6
0.7
OS 0.6 07 05
Vertical Section
(y 0)
30 40 50 60 70 80 90 100
0 10 20
Fig. 2.4 Calculated Isotherms of AT/AT0 • F0-4.4, A=0.35,k/u0
4.2 x IO'5, V/u0 0
23
-------
2) An entrainment region in which the centerline velocity and
temperature drop sharply, approximately as 1/x, as in a non-
buoyant jet. The jet spreads vertically by turbulent processes.
The lateral growth is dominated by gravitational spreading at a
much greater -rate than the vertical turbulent spread. Because
of this large ratio of lateral to vertical spread, the jet reaches
a maximum depth beyond which the bottom boundary rises to maintain
mass conservation (see Figure 2.4). Local densimetric Froude
numbers in this region decrease rapidly and the dilution
rises sharply as a result of entrainment (see Figure 2.3).
Surface heat loss remains negligible in this region, i.e.,
HT = 1.0.
3) A stable region in which vertical entrainment is inhibited
by vertical stability as indicated by the local densimetric
Froude number which is of order one or less. The jet depth
continues to decrease because of lateral spreading. The
small jet depths reduce the lateral entrainment, and the
dilution and centerline temperature remain relatively constant
in this region. The centerline velocity, however, drops sharply
as a consequence of the large lateral spread. The surface
temperature pattern is dominated by the wide, constant tempera-
ture stable region.
4) A heat loss region marking the end of the stable region. The
lateral spread is sufficiently large to allow significant surface
heat transfer and the temperature begins to fall again. Once
surface heat loss becomes significant, the rate of temperature
decrease is very rapid. However, at this point the centerline
velocity is so low that the discharge may no longer be considered
i
as a jet. In general it may be stated that surface heat loss, as
determined by the value of k/u , is not an important factor in
reducing the discharge temperature within the region of jet-like
entrainment. Beyond the stable region, temperature distribution
is controlled by passive diffusion processes acting upon the buoyant
plume.
24
-------
5) The primary effect of-a current in the receiving water is to deflect
the discharge without significantly affecting the dilution processes
up to the point where the jet centerline velocity excess, u , is
reduced to the magnitude of the cross current component in the
centerline direction. For u < Vcos6 the ambient water may have a
turbulence level comparable to the discharge turbulence and the basic
assumptions of the theory will not be valid. This region is properly
considered part of the far field.
The dilution and corresponding centerline temperature achieved in the
stable region are of interest since the stable region constitutes a significant
portion of the surface temperature distribution. Figure 2.5 is a plot of the
dilution and surface temperature in the stable region as a function of IF with
values of A indicated. It is clear that the ultimate stable dilution in a buoyant
jet depends primarily upon IF and to a lesser extent upon A.
The maximum depth of the discharge, h //h b , is also of interest and
6 max oo* '
maybe determined from the theoretical calculations. Figure 2.6 indicates that
the maximum vertical penetration of the jet is a function of IF with a very small
dependence upon A.
For quick calculation involving the maximum jet penetration, centerline
dilution, or centerline temperature rise, the results from Figure 2.5 and 2.6
can be condensed into three simplified formulas involving a new parameter, 3Fj ,
defined by
IF' = IF A1/4 = ——^—— ^2.23)
° ° / Ap1/2
./rjL, £,\ (h b )
Pa ° °
I" is merely a "Froude number" whose characteristic length is the scaling length,
(h°b )1/2, rather than a depth. The numerical results of Figures 2.5 and 2.6
o o
are then approximated by
AT _ ' * _- - (2.24)
AT"'
c
/(iF» ? + 1 ~ IF' (for IF1 >3)
V o o o
25
-------
—135
25
(F0- Densimetric Froude Number
Fig. 2.5 Centerline Temperature Rise ATC/ATO and Dilution, D,
in the Stable Region for k/u0 =V/u0 = 0.
26
-------
0 I I I I i i I i i I i i i I i i i i i I i i i i i i i i |
fFo - Densimetric Froude Number
Fig, 2.6 Maximum Vertical Depth of the Jet,
for k/u0 = V/u0 =0
27
-------
D = 1.4 /(IF1 )2 + 1 = 1.4 I" (for IF' > 3) (2.25)
stable V o o o
= o.42 IF' (2.26)
Ch b )1/2
o o
For IF' >3, h = (h + r) and (2.26) reduces to
o ' max max
hmax. ., = 0.42 IF' (2-27)
(h b )1/2
o o
Equations 2.24 and 2.25 relate to properties of the jet in the stable
region (stable centerline temperature rise and stable dilution respectively)
and (2.26) relates to a maximum property of the jet. For the spatial history
of the discharge or for information relating to lateral spreading or centerline
deflection under a current, use of the generalized plots of the theoretical
solution is suggested. Figures 2.7a to 2.7i give the basic calculated discharge
properties for a range of values of IF , A, and V/u with k/u set equal to zero.
As discussed previously, a non-zero value of k/u will have only secondary effect
upon the temperature distribution and setting k/u = 0 will always y±eld a
slightly conservative result (i.e. higher temperatures). The next two chapte-s
discuss the computer program and its application to actual discharge design.
28
-------
1000
T—I I I I-
Jet half width
b + s
100
1000
Fig. 2.7a Jet Parameters for F0=20.0, V/u0 =0 , k/u0 =
29
-------
1000
Jet half width
b+s
1000
Fig. 2.7b Jet Parameters for TF0= 10.0, V/u0 =0 , k/un=0
30
-------
ATC
ATn
O.I
0.05
0
h + r
hobo
10
1000
100
b-hs
'hobo
10
100 1000
1—I Mill 1 1 I I I I I I-
1.0
Jet half width
b + s
i i i i i i i i I I I I i i i I i I i I I I I I
10 100 1000
X
Fig. 2.7c Jet Parameters for F0= 5.0, V/u0 = 0 , k/u0=0
3l
-------
AT
.0
5
0
h+r
I 0
1000
00
b+s
0
100
1000
A= I.O/
Jet half width
b + s
2.0
10
100
1000
Fig. 2.7d Jet Parameters for FQ=2.0, V/u=Q , k/u =0
32
-------
ATC
ATn
1.0
.5
1000
h+r
hnbo 5
10
000
100
b+s
'o
10
Jet half width
b +s
TvTb
ouo
10
100
1000
Fig. 2.7e Jet Parameters for F0 = 1.0, V/u0 = 0 , k/u0=0
33
-------
10
1000
100
0
100 1000
T I I I I I I
Jet half width
b + s
Center line
deflection
1000
Fig. 2.7 f Jet Parameters for F0=20.0, V/u0-0.025, k/u0=0
-------
I.OF
0.5
ATC
AT0
O.I
0.05
0
1 1 1—I I I I I I 1 1—I—I MIL1
h + r
h i.o2.0'o.fc
10
1000
100
hnb
ouo
10
Jet half width
b +s
Center line
deflection
100
1000
Fig. 2.7g Jet Parameters for IF = IO.O, V/u = 0.05, k/u=0
0
35
-------
100 1000
- 2.0 Q5
1.0
I 0
IOOOF
00
Jet half width
b +s
hnb
ouo
A=2.0
Centerline
deflection
ho^o
10
100
1000
Fig. 2.7h Jet Parameters for F0=5.0, V/u0 = 0.1, k/uo=0
36
-------
ATC
AT.
1.0
.5
000
0
h + r
10
1000
100
hob
ouo
I 0
Jet half width
b + s
Centerline
deflection
«%/
y
10
ooo
Fig. 2.7 i Jet Parameters for F0 = 2.0, V/u0 =0.1 , k/u0= 0
37
-------
III. The Program
This chapter is a detailed description of the program for computing
the solution to the heated surface discharge equations. The solution method
is described including the basic equations, and the computational scheme.
Lastly, the computer program input and output format are discussed.
3.1 Dimensionless Equations
The computations are performed with the following dimensionless variables:
x = x/,/h b
o o
x = x/A b
n o
y = y/A b
J 00
u = u /u (3.1)
c o
"AT = ATC/ATQ
h = h/A b
o o
b = b/A b
o o
r = r/A b
o o
s = s/A b
o o
V = V/u
o
The equation set in Table 2.1 is reduced by eliminating the internal velocities
v , u, , w and w, as follows:
sb r h
a) Sum all the mass conservation equations to form a
total mass conservative equation.:
Vcos0 (s+b)(r+h)J - ul « rs+bl.) - a (r+hl,) I - 0 (3.2)
39
-------
b) Sum all the momentum conservation equations to form a
total -momentum.conservation equation:
— u2(s+bI9)(H-hI9) + 2uVcos8 (s+bl )(r+hl ) + V2 cos2 6 (s+b) (r+h)
dx L Z i
A~1/2 IF
d ^ + rhl, + h2I.)1 - uVcosefa (s+bl )- a (r+hl
/ j 4 L sz -1- y -1-
= 0
(3.3)
c) The total heat equation;
^ Tu
dx L
If" (s + bI7)(r + hi ) + AT" Vcos6 (s + bl ) (r + hi )
' ' ~> ~> J
AT (s + bi3) = o
(3.4)
d) The y-momentum spreading equation;
_1 ("ik _ e]|u
L
2bi (r+hl ) + 2uVcos6 bl
dx dx
V2 coS20b(r+h)
l
jj
AT (
= 0
(3.5)
e) The y-momentum deflection equation:
dx u (s+bI2)Cr+hI2)+2uVcoseCs+bI1)CJ+hI1)+V cos e Cs+b)(i+h)
0 (3.6)
40
-------
f) Combine the region 1 mass and momentum conservation
equations to eliminate v and w :
s r
(u+Vcose) — (u+Vcose) + IF 2 A"1/2
dx °
dx
:±^ + I ""•*• ° = 0
dx" 3 dx J
(3.7)
g) Combine the momentum and mass conservation equations from
regions 1, 2, and 3 to eliminate v , v, , w and w, :
S D u El.
(rb + ih) (u-I +Vcosel, ) ~
L dx
- ~ ) + Vcose) dVc°s6
I — — I
+ uVcose (I, - ~ ) (s ~ + r — ) + u Cl - Y- ) — l^s (u+Vcose)]
1 dx dx 1 dx
T h2 , T — r dr N , 1T dAT br2 2- dAT bh
H I^AT h —- ) + —I. — 1- I,, r +•
dx dx ^ J dx dx
"AT4 r2 + rhlj ^
Z J dx
ayr) — = 0
(3.8)
h) The geometrical relationships:
dx
dx
- sine = 0
(3.9)
- cos 0=0
dx
-------
Using the similarity functions defined in the theory, the integration con-
stants have the values:
I = .4500
I2 = .3160
I3 = .6000
(3.10)
I, = .2143
4
I5 = .2222
Ifi - .1333
I = .3680
The entrainment coefficients are as follows assuming e - .22:
a = .0295 r > 0
A/AC ~ A ct = a exp l-i
a = .0495 r= 0 sz z r L
z
(3.11)
a = -.0495 s =0
y
It should be noted that the computer program may be used for other choices of
similarity functions and spreading rate than those assumed here by simply
changing the program statements which specify the values of I -,-!-, and e . The
entrainment coefficients are automatically computed from these variables using
equations (2.18)
The above nine equations may be solved for the x derivatives of x, y, u,
AT, h, b, r, s, and 6. (Note that V is given as an input and thus is known.
The "initial" conditions at x = 0 are:
u = 1 - VcoseQ(V taken at origin) r = A1//2
AT = i ; = A-i/2
e = e
42 o
(3.12)
-------
are
3.2 The Computational Scheme
The computer program is conceptually simple. After all inputs
read the variables are initialized to their values at the discharge origin.
The solution proceeds by advancing along the discharge centerline (increasing
x) using the calculated derivatives of each variable to calculate its behavior.
The differential equation technique is a fourth order Runge-Kutta scheme which
has been taken, in modified form, from the IBM Scientific Subroutine Package.
The program consists of a main program and five subroutines:
MAIN; This program reads the input data, sets up the initial conditions
for each calculation, and calls the subroutine SRKGS which per-
forms the calculations.
SRKGS ; This subroutine is a modified version of the IBM Scientific Subroutine
Package DRKGS for solving a system of differential equations. The
form of the equation system is:
dy.
[a ] -J- = c (3.13)
^-J dx -1-
where a., is a coefficient matrix, y. is a vector of the variables
(u, AT, etc.) and c . is a vector of the constants in the i equations.
Subroutine SRKGS advances the solution by successively calling sub-
routine FCT which solves the system of equations for -- . The
results of the calculations are periodically printed by calling sub-
routine OUTP.
FCT_: This subroutine uses the current values of the variables, computes
the coefficient matrix, a.., the vector c., and solves the resulting
ij 'd . !
system of linear equations for -^p by calling routine SGELG. The
following details concerning this routine may be of interest: „
d fdbl d^b
a) To reduce the equation set to first order the equation — — - 2
2 dx IdxJ dx
is added where ^-| is computed from the y-spreading Equation (3.5)
db dX
and — — becomes a variable.
dx
b), Because of their relatively simple form, the y-def lection Equation (3.6)
and the geometrical relationships (3.9) are not included in the matrix
43
-------
solved by SGELG but are solved separately within FCT.
c) When an ambient current is present, (V^O) the value of c
is computed by setting 7? = 0 and solving a reduced set
of equations from which the heat equation (3.4) is omitted
, , db _ dh
and the y-spreading equation (3.5) is replaced by — - dx -
The full buoyant equation set is then solved using the cal-
culated value of e. When there is no ambient current, (V=0)
the spreading rate is e = .22, and the full equation set
is solved directly.
SGELG; This subroutine is a modified form of the IBM Scientific
Subroutine Package routine DGELG. J-t solves the system of
linear equations by Gauss-Jordan reduction.
QUIP; The values of the calculated variables (y.) are printed in
formatted form. Certain other quantities of interest are also
printed (see Section 3.4).
CROSS: This subroutine is called by the other routines whenever
the velocity of the cross flow is required. The routine uses
the input values V, to V0 to compute the cross flow as a function
— .
of it. This subroutine may be modified by any user to any par.ticu-
lar functional form with the only requirement being that the
routine place the value of the cross flow at the current value
~ dV
of x in V and the value of- — in DV. The cross flow function
dx
used in the present version is described in the next section.
The computations for a particular case will be terminated for one of the
following reasons:
1) the limit on x specified by the input has been reached (see
next section). This is the normal termination.
2) u < V cos 9: This termination occurs when the centerline velocity
excess is reduced to the same magnitude as the ambient current in
which case the basic assumptions of the theory are not valid.
3) u < .02: The limiting value .02 for u is an arbitrarily chosen small
number below which the assumptions of jet-like behavior are not valid.
44
-------
4) The total dimensionless momentum, M = u (r + hi )(s + bl ) +
2uVcos6(r + hi ) (s + bl.) + V2 cos20(r + h)(s + b) + F~2A~1/2
— — — 1 —2 —9 °
AT (s + bl ) (-• r + rhl, + hi.), which should be constant for
V V
•^- = 0 and nearly constant for — near zero has deviated from its
o
u
o
initial value by more than 25%. This termination indicates that
the computation is accumulating a large numerical error.
3.3 Input Formats
Input Data Format
Number of cases to be calculated 13
One Set IF, A, k/u , 6 , JL , ERR, STEP 2F10.5, F10.7, 4F10.5
.. , O O O Li
for each
Calculation V , V?5 V , V. , V , V , V , VQ 8F10.5
J_ ^ -j ^4 .J O / o
u u
where IF = initial densimetric Froude number = =
o
SAT gh
o o
P °
a
A = aspect ratio = h /b
oo
k/u = surface heat loss parameter
6 = initial angle (in degrees) between the discharge centerline
and the boundary of the ambient region
x = the value of x/i/h b at which the program should terminate
L o o
for this case
ERR = maximum allowable average error in the variables at each time step
STEP = interval of x at which variable values should be printed
V -V0 = constants describing the cross flow. In the present program
1 ,8
version only V, - V,. are used as follows
T r ^ 12!
V = V1 + V2 exp I- V3|V4X - Vj J
r » ^2^
- v (v. x - v )
|_ 3 4 5 J
A sample input listing is given in Appendix III.
45
dx 4 5234
-------
3.4 Output Format
The output is paged for each calculation with the basic input variables
printed on the heading of each page. The values of the variables are printed
in column form under the following headings
X = x
H = h
•B = b
R = r
S = i
- - 1/2
L o -2
Q = u(r + hl^Cs + bl±) + VcosG (r + h) (s + b)
M = u (r + hI2)(s + bI2) + 2uVcos0(r + hi ) (s + bl )
—2 2 — — — — 9 —1 /9 1—9
+ V cos (r + h)(s -F b) + IF ZA a/^ AT(s+bIQ) (^r + rhl. +
o j 2 3
U u
T = AT
HT = AT [u(r + hly) (s + bl?) + Vcos8(r + h) (s + b)]
V = V
XP = x
YP = y
THD = 6
TM =
= f -
u
An example of program output is shown in Appendix III.
46
-------
IV. Application^
The theoretical model presented in the preceding sections permits
determination of the behavior of a heated surface discharge as a function of
relatively few controlling parameters: F, \/bQ, k/u , V/u . The geometry
of actual field configurations are rarely as simple as those assumed in .the
theory and judgement must be applied in the schematization of the discharge
to a form for which the dimensionless parameters may be given. Similarly,
the output of the program must be interpreted in the light of the basic assump-
tions of the model. The following sections discuss in detail the generation
of input data for the computations and the construction of the temperature
distribution from the program output. Finally, a case study is presented as
an example of the use of the theory and computer program.
4.1 Schematization
This section is a discussion of the data requirements and schematization
techniques for preparing input-to the program. The following are the physical
data needed as input to the theoretical calculations.
The ambient temperature, T , is assumed in the theory to be constant
a
in space and time. Actual ambient receiving water temperatures are often
stratified vertically or horizontally and may be unsteady due to wind, tidal
action or diurnal variations in solar heating. The natural stratification of
the ambient water may be increased by accumulation of heat at the x^ater surface
if the discharge is in a semi-enclosed region. The value of the ambient tempera-
ture, for a given initial temperature rise, AT , determines the initial density
difference between the discharged and ambient water. Also the effectiveness
of the entrainment of ambient water in reducing the discharge temperatures is
a function of the ambient stratification.
I
If the temperature differences resulting from ambient stratification in
the vicinity of the discharge are the same magnitude as the initial discharge
temperature rise, the theoretical model of this study should not be applied
without further development to account for the ambient stratification. The
theory is valid if an ambient temperature, T , may be chosen which is represen-
cl
tative of the receiving water temperatures, that is, if the temporal or spatial
variations in ambient temperature do not differ from T& by more than a few
degrees Fahrenheit.
47
-------
The initial discharge temperature rise, ATQ, is determined from the
discharge temperature, T , and the chosen ambient temperature, T&. The
value of AT will be equal to the temperature rise through the condensers only
if the intake temperature is equal to T . Actual discharge temperatures should
3.
be relatively steady unless the power plant has several different condenser de-
signs in which case the discharge temperature will vary with the plant load.
Use of the theory requires that a constant value of AT be specified, the choice
being based on the likely steady value of the discharge temperature,
The initial relative density difference, Ap /p may be related to T
—- O 3. Si
and AT by Ap /p = gAT where 3 is a function of T as given in Figure 4.1 for
o o a o
fresh water. The curve of 3 vs. T is not linear but only a small error will be
introduced if T is set equal to T + AT /2. A more accurate value of Ap /o
M a o o Ka
may be obtained by using tabulated values of water density vs. temperature, or
in the case of salt water, water density vs. salinity and temperature.
The initial condenser water discharge velocity, u , is a function of the
power plant condenser water pumping rate and the discharge channel area.
The discharge channel geometry is important since the theory uses the
square root of one half of the discharge channel flow area as a scaling length.
Calculation of the discharge densimetric Froude number, IF , requires specifi-
cation of the initial depth, h ; and the aspect ratio, A, requires both h
and the initial width, b . The following procedure is suggested for any channel
shape:
a) Let hQ be the actual maximum discharge channel depth
so that calculation of IF is not affected by the schema-
tization.
b) Let bQ be such that the correct discharge channel area
is preserved:
channel area
__
o
Then the aspect ratio is given by:
h 2h 2
A - —2- - °
b channel area (4.2)
As an example, the aspect ratio of a circle is 8/ir = 2.55.
48
-------
.0002
.000 I 5
.0001 -
.00005 —
Note : Function is not
quite ! inear
i i i M i i 11 i i i 11 i 11 11 i i 11 M i 11 j i i 11 11 i i 11 i i 11 i.i i 11 I i i 11 i 11 11
40 50 60 70 80 90 100
T(°F)
Fig. 4-1 Coefficient of Thermal Expansion for Fresh Water
Q - JL 5L£— (op"1) as a Function of Temperature T(QF)
P 3T
49
-------
The discharge channel geometry may vary with time if the elevation
of the receiving water changes due to tidal motion or other causes. In
this case separate calculations for each elevation of the receiving water
must be made, assuming that the steady state theory predicts the instan-
taneous temperature distribution.
Very complicated arrangements of the discharge channel relative to
the boundaries of the ambient receiving waters are beyond the capabilities
of the theory of this study. If the discharge channel is nearly at right
angles to the solid boundaries, the theory may be used as developed, or if
the discharge is directed parallel to a straight boundary, the discharge may
be schematized by assuming that the solid boundary *Ls the centerline of a jet
whose discharge channel is twice the width of the actual discharge (see Figure
4.2). The width b is then twice the width calculated by the procedure discussed
o
previously. However, if it is not clear whether the jet will entrain water from
one or both sides or if irregularly shaped solid boundaries will deflect the
jet or distort it from the form assumed in the theory, a meaningful schematiza-
tion is not possible.
The densimetric Froude number IF and the aspect ratio A of the s'urface
o
discharge channel should be chosen to be consistent with the bottom topography
of the receiving water body. If optimum dilution is to be obtained, the ver-
tical development of the surface jet should not be limited by the bottom of
the receiving water. This is generally considered to be a. desirable objective
by the aquatic or marine biologist inasmuch as it avoids exposure of benthic
organisms to high velocities and temperature rises. Designs which contain minimum
interaction between the surface jet and the bottom topography are shown on
the left side of Figure 4.3. Examples of surface jets in which there are substantial
interaction? in relation to the bottom topography are shown on the right side of
Figure 4,3. ,In the latter case, vertical entrainment and dilution are reduced
by the interference of the jet and the bottom.
If the discharge channel densimetric Froude number, IF is less than
o
unity, a wedge of ambient water will intrude into the discharge channel and
the heated flow will be forced to obtain a densimetric Froude number of unity
at the discharge point (see Figure 4.4). The depth of the heated flow in the
presence of a wedge, h *, is given by:
50
-------
bn L d
I////
7 7 7 7 T
00 = 9 0
GQ near 90°
Discharges with Entrainment from Both Sides
i i i / r
Discharges with Entrainment from One Side
—" f
Schematization Possible
Schematization not Possible
Because of Irregular Geometry
Discharges with Obstructions
Fig. 4.2 Discharge Channel Schematization
51
-------
Minimum Interaction
Substantial Interaction
Rg. 4.3 Limitations on Maximum Jet Thickness by Bottom
Topography
52
-------
Ln
= i.o
Fig. 4.4 Two Layer Flow in the Discharge Channel
-------
2/3 (4.3)
h o
o
where h and IF are based on the channel dimensions. A flow area may be
o o
calculated based on the depth h * and the aspect ratio calculated as des-
cribed previously in this section.
The surface heat loss coefficient, k, may be estimated from local meteor-
ological variables, principally wind speed, and the ambient water temperature.
Care should be taken to use values for k which are appropriate for local condi-
tions. Edinger (13) and Ryan (16) discuss methods for determining values of
K = pck.
The cross flow velocity, V, in the receiving water may be measured
directly or estimated from flow measurements in the case of a channel flow.
The theory accepts as input values of V/u as a function of x/A b (see
o o o
Chapter 3) .
Once the schematization of the discharge configuration is achieved, the
theoretical calculation is performed by the computer program described in
Chapter 3. The inputs to the program for each calculation are IF A = h /b
^ O 00
k/u and V/u (as a function of x//h b ) .
o o o o
It should be noted that the theoretical model may not necessarily be
applied to any arbitrary set of input parameters. In general the value of
k/uQ has little effect; however, depending upon the value of the other para-
meters, IFQ , A, and V/UQ, and upon the specified error size, ERR, the program
may fail to generate a solution. This problem may take one of the following
forms:
1) Little or no advancement of solution, i.e. the step size
remains very small.
2) Abrupt discontinuities noticeable either in the values of the
variables or in the total mass, momentum, or heat.
3) Unstable behavior in one or more variables, culminating usually
in an overflow error.
The chance of the program encountering one of the above problems increases
with decreasing 1FQ and A and increasing V/UQ. The value of ERR produces no
consistent result except that too large a value of ERR most often yields
54
-------
unstable behavior or discontinuities, and too small a value may prevent
advancement.
The causes of this anomolous behavior in certain cases are not totally
understood at the present time, but are thought to be related to the following:
1) accumulation of round-off error where the solution should be
stable, i.e. a totally numerical problem,
2) genuinely unstable solutions, probably caused by the deter-
minant of the differential equation coefficient matrix being
zero or near zero.
The first of thes'e could probably be solved by paying greater attention to round
off error generation than is now done. To the extent that the second cause has
no physical significance it may also be desirable to eliminate this problem by
numerical manipulation. However, the governing equations for the heated dis-
charge are similar to open channel flow equations which possess critical flow
points at which the equations strictly have no solution and the coefficient
matrix is zero. It is not clear to what extent the (degeneracy of th& three-
dimensional heated discharge equations corresponds to critical flow behavior;
the answer will only come by observing actual discharges in the laboratory
and field. Thus for the time being the limitation on the computational range
of this theory must be accepted.
Of course the values of IF , A, and V/u are pre-determined by the case
being considered, leaving ERR as the only free parameter. Figures 2.7a to 2.7i
indicate the range of each variable and the various combinations of values for
which calculations may be successfully performed. Table 4.1 gives a rough
guideline as to the best value of ERR to be used for different cases. In
general if the solution is not advancing, a larger ERR might be tried and if
the solution appears unstable a smaller value should be tried. Experience has
shown that the behavior of the solution may also depend upon the type of com-
puter used. The results quoted herein are based on calculations done on an
IBM 360 either -65 or -75. Users employing different systems than this are
encouraged first of all to try several values of ERR in an attempt to make an
unsuccessful calculation work and secondly to take up the challenge of illumin-
ating more clearly than is done here the properties of the governing equations,
the physical relevance of these properties, and their efficient numerical treat-
ment .
55
-------
Table 4.1
Guideline to Choice of Maximum Error Value. ERR
IF
o
1-2
0.1 - 0.5
.005 with no
crossflow
0.5 - 1.0
.005
l.U - J..O
.005
z. u — °°
.005
.05 with cross-
flow
.005 with no
2 - 5 crossflow
.05 with cross
flow
5-10 .01 .01 .01 .01
10 - °° .01 .01 .01 .01
4.2 Use of the Program Output
An example of the program output is shown in Appendix III. All of the
outputs are in dimensionless form and must be transformed by the following steps
1) Multiply x, y, x. r, s, b. and h by yhT b
_ ° °
2) Multiply u by u
3.) Multiply AT by AT.
4) Multiply TM by /h b /u
The calculated isotherms may be constructed as follows:
1) Locate the centerline using x and y.
2) Assign values of T = T + AT along the centerline.
C 3. C
3) Plot the boundaries of the core and turbulent regions using
s, b, r and h.
4) Assign temperatures within the discharge using the assumed
temperature distribution (Equation 2.5).
56
-------
Figure 4.5 may be used as an aid in the construction of isotherm contours;
lines of constant r = id~s 2-r
T h ^z " ~TT~ are Pl°tted and then the temperature
contours are drawn by referring to the figure. The values of temperature rise
above ambient are expressed as fractions of ATQ, the initial temperature rise.
If the surface area within a given temperature rise AT* is of interest,
the following formula may be used:
where
-
*
(1 -
•y-
A* = 2hobo J
s = s//h b
o o
b = b//h b
o o
AT = AT /Ai b
C 00
AT,= AT./AT
x..= the value of x//h b where AT = AT,
* o o *
(4.4)
As given in the program output as a function of x.
A simple computer program may be written which performs the above integration
numerically.
Finally, the travel time along the centerline to a given temperature',
t • •••
Ta +> AT*j is §iven by ™ (5r,^ - •
It is important to note that although the theory calculates the tempera-
ture distribution out to very small values of temperature rise above ambinet,
these extreme regions of the discharge will be subject to distortion by random
ambient processes especially wind stresses.
The calculated temperature distribution must be interpreted with this
in mind and no critical significance should be given to the exact calculated
position of isotherms of small temperature rise.
4.3 Case Study_-_H_eated Surface Discharge into a Receiving Water Body
of Finite Depth
Consider a proposed discharge into a large body of water at ambient
temperature T with a uniform depth H as shown on Figure 4.6. A constraint is
3.
57
-------
ATC
AT
l-n
OO
Fig. 4.5 Temperature Distribution Plotting Aid
-------
hn
-------
placed upon the discharge velocity such that UQ < UQ. Because the theory
applies only to discharges which are not influenced significantly by solid
boundaries in the receiving water, it is desirable that h^^ < H where h^
is the maximum vertical penetration of the discharge (see Figure 4.6). The
choice of h as the "critical" boundary of the jet is an arbitrary one but
max
is probably conservative considering the assumed structure of the jet which
sets u = 0 at that point. Entrainment should not be significantly affected
by contact with the bottom at the point of maximum depth only.
A third constraint may be imposed by local topography. For instance
h must be less than or equal to the water depth, and 2bQ may be constrained by
some critical width. Thus in general h <; R and bQ< S-
The discharge is characterized by an initial temperature rise, ATQ,
and a condenser water flow, Q . The problem treated in this case study is to
design a discharge channel for maximum ultimate (stable) dilution (see Figure 2.5)
while meeting all the imposed constraints. This case study will illustrate the
generation of temperature contours and isotherm areas from the computer output
in addition to treating the more specific, but often relevant, problem of optim-
izing a channel design with one or more constraints imposed.
The definition of IF' and the following approximate formulas may be fe-
o
called from Chapter 2 for IF1 greater than 3.
, (h b
o o o
IF' = ff A= - r - -_' (4.5a)
(4.5b)
2(g6AT)L/2(hb )5/4
00
1/4u 5/4
Up _ , .
1/2 — I7T~ (4.5c)
AT
^p) ~ I?; (4.6a)
c
s
60
-------
Ds = I-* n (4.6b)
h
max
= 0.42 IF' (4.7)
1/0
(hb J1/2
o o
Equation (4.5) defines IF^ , (4.6a and b) relates centerline dilution and overall
dilution to l^and (4.7) relates the maximum jet penetration to IF' . Note that
for a given Apo/p and QQ both the maximum stable dilution and the maximum pene-
tration are dependent only on discharge velocity u and not on the aspect ratio.
This provides freedom in selecting channel geometry.
The constraint that h < H may be restated using (4-7) as
6.8
_
O
and from (4.5c) the constraint that u <; U may be restated as
o - o
1.19 U 5/4
Because ultimate dilution increases monotonicly with IF' (see Eqn. 4.6),, the
channel should be designed witn IF' equal to the smaller of the two expressi
in (4.8) and (4.9). From (4.5b) it may be shown that
Q4/5
h b = -r-7= - 7^ - -_ and (4.10)
u = ^-^- (4.1D
o
where IF'is picked from (4.8 and 4.9). If Chere is a constraint on either the
width or depth of the discharge structure, this constraint may be included with
(4.10) to determine h and b . If there is a constraint on both width and depth
o o
61
-------
such that h < R and b < S and the product RS < h b as determined by (4.10),
o o o o
then a discharge channel cannot be designed with the desired constraints.
Example: Design a discharge channel no more than 12 feet deep to
achieve maximum stable dilution o± a neated discharge of
2,000 cfs and 15°F temperature rise into a water body 31 ft.
deep at an ambient temperature of 70°F. Maximum discharge
velocity is U =6 ft/sec. The area within the 4°F surface
isotherm should be calculated.
1. The value of g is taken from Figure 4-1 using T = 77.5°F, g = .000144°F
and (ggAT ) = .069 ft/sec .
2. Using (4.8) the constraint on depth is
IF' < 6.8(31)5/3(.069)1/3 , ,
o - = 5'3
2000
Using (4.9) the constraint on velocity is
, ,
b.4
r— r~ - , ••/,
C.069) ' (2000) '
Hence IF1 should equal 5.3 and the jet should just touch the bottom.
3. From (4.10) and (4.11)
h b = 2000^ = 2
24/5(5.3)4/5(.069)2/5
1/2
(h b ) = scaling length = 13.9 ft.
O O
and U0 = = 5-2 ft/sec
oo
4. No constraint has been placed on b so any number of combinations of h b
0
will be satisfactory. As an example, h = 11 ft. and b =17.5 feet
o o
satisfies the requirement on channel depth and results in a 3F =60-
o '
aspect ratio, A S 0.6; and IF' = 5.3.
62
-------
Theoretical calculation (compute! output) for these values are shown
in Appendix III.
Figure 4.5 is used to plot the isotherms in two planes (horizontal at
the surface and vertical at the j = .27 and XA = 4'5_. The computations are organized as follows:
1.131 1.52 .585
2.62 1.37 1.16
5.24 1.44 2.39
10.5 1.48 5.27
« • (
41.9 0 41.5
AT
AT^/AT
s+b(l-ATA/AT)2/3
.971 .278 .8
.940 .237 .7
.884 .305 .7
.692 .390 ,7
. •
9 *
.282 .957 .1
05 1.99 2.25
98 2.30 5.67
84 3.31 14.3
19 5.27 39.3
o •
'
22 5.0! 267.
The first and last column may be multiplied by the scaling length, /hQbQ = 13.9 ft.
and plotted against each other to produce the 74° isotherm shown inFig. 4.7. The inte-
gration in Equation 4.4 results in an AA = 2(192) (267) = 103,000 sq. ft. or
2.4 acres.
63
-------
-500k
500
1000
x(FT)
b0=!75ft AT0 = 15° F
u0= 5.6 ft/sec
Fig. 4.7 Calculated Contours for the Case Study Example
64
-------
ACKNOWLEDGEMENT
This study Was partially supported by the Environmental Protection
Agency under the National Thermal Pollution Research Program. The coopera-
tion of Mr. Frank H. Rainwater, Chief of the National Thermal Pollution
Research Program, is gratefully acknowledged.
Additional support for the study was received from the Public Service
Electric and Gas Company of Newark, New Jersey, as part of a study of thermal
discharges for the proposed Atlantic Generating Station, an offshore floating
nuclear power facility.
Mr. Patrick Ryan, Research Assistant in the R.M, Parsons Laboratory
for Water Resources and Hydrodynamics made substantial contributions in the
development of the relationships used in Chapter IV. Numerical computations
were carried out at the M.I.T. Information Processing Service Center. Our
thanks to Miss Kathleen Emperor who typed the report.
65
-------
References
1. Fan, L. N., "Turbulent Buoyant Jets into Stratified or Flowing Ambient
Fluids", W. M. Keck Laboratory of Hydraulics and Water Resources,
California Institute of Technology, Report No. KH-R-15, June 1967-
2. Fan, L. N. and N.H. Brooks,"Numerical Solutions of Turbulent Buoyant Jet Problems",
W. M. Keck Laboratory of Hydraulics and Water Resources, California
Institute of Technology, Report No. KH-R-18, 1969.
3. Tamai, N. and others, "Horizontal Surface Discharge of Warm Water Jets",
ASCE Journal of the Power Division, No. 6847, PO 2, October 1969.
4. Wiegel, R. L. and others, "Discharge of Warm Water Jet Over Sloping Bottom",
Technical Report HEL-3-4, University of California, November 1964.
5. Jen, Y. and others, "Surface Discharge of Horizontal Warm Water Jet",
Technical Report HEL-3-3, University of California, December 1964.
6. Stefan, H. and F. R. Schiebe, "Experimental Study of Warm Water Flow
Part I", Project Report No. 101, St. Anthony Falls Hydraulic Laboratory,
University of Minnesota, December 1968.
7. Hayashi, T. and N. Shuto, "Diffusion of Warm Water Jets Discharged
Horizontally at the Water Surface", IAHR Proceedings, Ft. Collins,
Colorado, September 1967.
8. Ellison, T. H. and J. S. Turner, "Turbulent Entrainment in Stratified
Flows", Journal of Fluid Mechanics, Vol. 6, Part 3, October 1959.
9. Hoopes, J. A. and others, "Heat Dissipation and Induced Circulations
from Condenser Cooling Water Discharges into Lake Monona", Report No. 35,
Engineering Experiment Station, University of Wisconsin, February 1968.
10. Motz, L. H. and B. A, Benedict, "Heated Surface Jet Discharged into a
Flowing Ambient Stream", Department of Environmental and Water Resources
Engineering, Vanderbilt University, Nashville, Tennessee, Report No. 4,
August 1970.
11. Morton, B. R., Taylor, G. I. and J. S. Turner, "Turbulent Gravitational
Convection from Maintained and Instantaneous Sources", Proc. Roy. Society,
London, A 234:1-23, 1956.
12. Morton, B. R., J. Fluid Mech. 5:151-63, 1959.
13. Edinger, J. E., Duttweiler, D. W. and J. C. Geyer, "The Response of Water
Temperature to Meteorological Conditions", Water Resources Research, 4:5,
1137-43, 1968.
67
-------
14. Abramovich, G. N., The Theory of Turbulent Jets, TheM.I.T. Press,
M.I.T., Cambridge, Massachusetts, 1963.
*
15. Stolzenbach, K. D. and D.R.F. Harleman, "An Analytical and Experimental
Investigation of Surface Discharges of Heated Water", R. M. Parsons Labora-
tory for Water Resources and Hydrodynamics, Technical Report No. 135,
Department of Civil Engineering, M.I.T., February 1971.
16. Ryan, P. J. and K. D. Stolzenbach, Chapter 1: "Environmental Heat Transfer",
Engineering Aspects of Heat Disposal from Power Generation, (D.R.F. Harleman,
Ed,.), R. M. Parsons Laboratory for Water Resources and Hydrodynamics, Depart-
ment of Civil Engineering, M.I.T., Cambridge, Massachusetts, June 1972.
*
Also published as Water Pollution Control Research Series Report No. 16130 DJU 02/7.
by the Water Quality Office, Environmental Protection Agency, Washington, D. C.
68
-------
List of Symbols
A - discharge channel aspect ratio, h /b
o o
a.. - coefficient matrix in the governing differential equation set
b - horizontal surface distance from core boundary to jet boundary
b - one half the width of rectangular discharge channel
C - coefficient in the exponent of Ellison and Turner's vertical
entrainment velocity function
c . - vector of the constants in the governing differential equation jet
D - ratio of flow in the jet to the initial flow = dilution
D - dilution in the stable region of the heated discharge
s
QL _ lateral spread of the turbulent region of a jet
CLX
ERR - maximum allowable average roundoff error for each step in
the numerical computation
IF - densimetric Froude number of the discharge channel = —
° uc /^Po
IF - local densimetric Froude number in the iet = /—gh
L /- 7- v Ap o
Ap_ Sh
1/4
IF' a characteristic Froude number = IF A V p
o o
f - similarity function for velocity =(!-<; )
g - acceleration of gravity
HT - ratio of heat flow in the jet to the initial discharge heat flow
H - maximum allowable vertical penetration of the jet
h - vertical centerline distance from core boundary to jet boundary
h - depth of the discharge channel
o
h* - depth of the heated flow at the point of discharge if a cold
water wedge is present
h - maximum value of h obtained in a heated discharge
max
1
i - f f(e)de = .4500
Jo
-i
-i
= .3600
= .6000
69
-------
1 1
ft t(t)d?dC = .2143
J
f f(?)C ' d? = .2222
Jo
r 2 1/2
I,- - J f (£)e de = .1333
6
1
J f(e)t(?)d? = .3680
K - surface heat loss coefficient
k - kinematic surface heat loss coefficient
P - pressure
Q - discharge channel flow
R - maximum allowable depth of discharge channel
r - vertical distance from the jet centerline to the boundary
of the core region
S - maximum allowable half-width of discharge channel
s - horizontal distance from the jet centerline to the boundary
of the core region
STEP - increment in x//h b at which numerical output is printed
T - temperature
T - ambient temperature
3.
T - jet centerline surface temperature
c
T - temperature of the heated flow in the discharge channel
TA - a particular temperature of interest
AT - temperature rise above ambient in the jet, T - T
cL
AT - surface temperature rise above ambient at the jet centerline, T -T
C C cl
AT - temperature difference between the discharge and the ambient water,
T -T
o a
* * a
TM - dimensionless time of travel along the jet centerline,
U. X
o r dx
/h b
o o
3/2
similarity function for temperature = (l-£ )
70
-------
u,v,w - velocity components in the coordinate system relative to the
centerline of a deflected jet
u,v,w - velocity components in the fixed coordinate system
UQ - maximum allowable velocity in the discharge channel
u - velocity in the discharge channel
u* - velocity of the heated flow at the point of discharge if
a cold water wedge is present
u - surface centerline jet velocity
V - ambient crossflow velocity
V- -V - ambient crossflow velocity components as input to the numerical program
v - lateral velocity of the entrained flow at the jet boundary
v - lateral velocity in the jet at y = s+b and -(h+r)
-------
C - either of £ or t,
y z
6 - angle between the jet centerline (x axis) and the y axis
6 - angle between the discharge channel centerline. and the y axis
p - density of water
p - density of the ambient water
3.
p - density of the heated discharge
Ap - difference between the ambient water density and the water density, p& - p
Ap - difference between the ambient water density and the density of
the heated flow in the discharge channel, p - p
3. O
superscript ' - indicates a turbulent fluctuating quantity
superscript - indicates a dimensionless quantity. Note that variables
printed by the computer are dimensionless but are
capitalized.
subscript s - indicates a jet property in the stable region
72
-------
Appendix I Flow Chart for Solution of [aij] ——
INPUT F0 , A, K/u0 , B0 , XF, ERR, STEP, V(I)
X Xf
A
X = X + AX
*
X =
*
dyi
Compute aij , Cj , — -*- , V
*
Compute yj*
X AX
i
AX- -AK-
ave. truncation error
-------
DOUBLE PRECISION OERY,Y,PRMT,X,TH,XP,YP,U,6,T,h,BX,S,R,XLI^,AUX,XX
DOUBLE PRECISION A ,C , CCGS,DSIN,UELX
REAL II,12, 13,14,15, 16,17, 18,M
INTEGER P,RR,REG,VFG,FLAG
DIMENSION PRMT(5) ,Y( 1C) ,AUX( 8, 10),OtRY( 10) , A(1C,10),C( 10),VtL(8)
EQUIVALENCE (Y{8),TH),
-------
XP = 0. PGM10037
YP=0. PG-M10J38
CALL CROSS PGM10039
U=1.0-V*COS(THU PGM 10040
IER=0 PGM10041
TM = 0 PGM10042
XI=0 PGM100t3
NPAGE=0 PGM10044
NLINE=50 PGM1Q045
REG=1 PGM0046
FI=(1./(FR**2) )*S PGM10047
RMJ=1.0*R*FI*.5 PGM10048
PRMT{1) = 0. PGM10049
PRMT(2) = XLIM PGM10050
PRMT (3 J=. 00001 PGM10051
IF (VEL(IJ) 2,1,2 PGM1GUJ3
1 NOIM=7 PGM10054
VFG=2 PGM10055
GO TO 3 PGMlOOiiC
2 NDIM=10 PGM10057
VFG=1 PGMlOOiJiJ
3 DO 4 N=1,NOIM PGM10059
DERY(N)=1./NDIM PGM10060
4 CONTINUE PGM10061
CALL SRKGS PGM10062
IFiK-KK) 7,8,8 PGM10063
7 K=K+1 PGML0064
GO TO 10 PGM10G65
8 CALL EXIT PGM100t>6
9 STOP PGM10067
END PGM10068
-------
c
c
c
c
SUBROUTINE SRKGS
DGU4LE PREC I SIOlM PRHT ,Y,OLRY , AUX ,A,fi ,C , X ,X E,\l C , h , A J t B J , C J , RL, R2,
lCELT,DABS,XLIMtXX, AUN ,CUN, DC OS , OS IN, DEL A
RFAL I 1, 12, 13, 14,15, 16,17, Id , ','•
INTEGER P,KR , REG, VFG , FLAG
DIMENSION A (4) , B(4), C (4 j , AUN ( 10, 10) ,CUN{ 10)
DIMENSION PRMH5),Y(1C) , AUX{ d,10),!)EKY{103,VEL(i)
COMMON DERYT Y, PRMT , X , XL I M, AUX, XX f AUN, GUN ,OELX
COMMON THIO,THI ,ERi; , V t'L , 1 1 , I 2 , 13 , 1 4, I 5 , I t , I / , I 8 , E P S, V , 3X S, UT , F I
COMMON PtKf<,NPAGF:,NLINE,KEG, NU I'w, ,VFG , FL AG , I HLF, :> TE Pf TM , X I
DC 1 1=1, ND I N
AUX( d, I ) = . 0666666 6 tact 66666 7 DO*DtRY ( 13
X=PRMT (1 )
H=PRNT (3 }
PRMT (5)=0.00
CALL FCT
ERROR TEST
IF
-------
C PREPARATIONS OF FIRST RUNGE-KUTTA STEP SRKGG037
CU 3 I = 1,NDIM SRKGQ03t!
JHJX{1,I)=Y(I) SRKG0039
AUX( 2,IJ=DERY( I ) SRKG0040
AUXOt I) = O.DC SRKG0041
3 AUX(o,I)=O.DO SRKG0042
IREC=0 SRKG0043
h=H+H SRKG0044
IHLF=-1 SRKG0045
ISTEP^O SRKG0046
IEND=0 SRKG0047
C SRKG0048
C SRKGOQ49
4 IF ( ( >+H-XENDJ*H) 7,6,5 SRKG0050
C START OF A RUNGE-KUTTA STEP SRKG0051
5 H-XEND-X SRKG0052
6 IEND=1 SRKG0053
C SRKG0054
C RECORDING OF INITIAL VALUES OF THIS STEP SRKG0055
7 CALL GUTP SRKGOOSu
36 IFiPRMT(5))40f8,40 SRKG0057
8 ITEST=0 SRKG0058
9 ISTEP^ISTEP+1 SRKG0059
C SRKG0060
C SRKG0061
C START OF INNERMOST RUNGE-KUTTA LOOP SRKG0062
j=l SRKG0063
CO 91 I=1,NDIM SRKG0064
AUX(3,I)=0.00 SRKG0065
91 AUX(6,I)=O.DC SRKG0066
10 AJ=A(J) SRKG0067
BJ = B(J) SRKG0068
CJ=C(J) SRKG0069
DO 11 I=1,NDIM SRKG0070
R1=H*OERY(I) SRKG0071
R2=AJ*(Rl-BJ*AUX(o, IJ ) SRKG0072
-------
c
c
c
c
c
c
c
11
12
13
14
15
16
17
18
=Y(I )+R2
R2=R2+R2+R2
AUX(6, I)=AUX (6,1 )+R2-CJ*Ri
IF(J-4)12,15,15
IF( J-3) 13,14,13
X=X+.5DO*H
CALL FCT
GO TO 10
END OF INNERMOST
SRKC0073
19
20
21
22
RUNC-E-KUTTA LUGP
ACCURACY IS IMPOSSIBLE
TEST OF ACCURACY
IF( ITEST ) 16, 16,20
IN CASE ITEST=0 TESTING OF
DO 17 I=1,NDIM
AUX( 4, I)=Y( I )
ITEST=1
ISTEP=I5T£P+ISTEP-2
IHLF=IHLF + 1
X=X-H
H=.500*H
DO 19 I=1,NOIM
Y( I )=AUX(lf I >
DERYd )=AUX( 2, I )
£UX<6, I)=AUX(3t I)
GO TO 9
IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE
IMOO=ISTEP/2
IF( ISTEP-IMOC-IMOD)21,23,21
CALL FCT
DO 22 I=1,NDIM
4UX(5,I)=Y( I)
AUX( 7,I)=DtRV( I )
GO TO 9
SR KG 00 75
SRi KG 00 76
SRKGJJ77
SRKG007o
SRKG0079
SRKG0080
SRKG0081
SRKG0082
SRKG0063
SRKGODSn
SRKQ 0036
SRKGOJcJ7
SRKG0038
SRKGOOdV
SRKG0090
SRKG0091
SRKG0092
SRKGC093
SRKG0094
SR KG 00 96
SRKG0097
SRKG0093
SRKG0099
SRKGU100
SRKG0101
SRKG0102
SRKG0103
SRKG0104
SRKG0105
SRKG0106
SRKG0107
SRKG0103
-------
C SRKGOlJSi
C COMPUTATION OF TEST VALUE UELT SRKGuilO
23 CELT=O.DO SRKGOlli
DO 24 I=1,NDIM SRKGOll^
24 CELT=DELT+AUX(S, n*D0ES(AUX(4,I)-Y( I )) SRKG0113
IF(DELT-PRMT(4))23,28,25 SRKG 01 14
C SRKG0115
C ERROP IS TOO GREAT 5RKG0116
25 GO TO 26 SRKGOL17
26 CO 27 I=1,NOIM SRKG3113
2.7 AUX(4,I)=AUX{5,1) SRKG0119
ISTEP=ISTEP + ISTEP-4 SK
-------
ISTEP=lSTEP/2 SRKG0145
H=H+H SRKG01*6
33 IMCD=ISTEP/2 SRKGQ147
IF( ISTEP-lMOO-IMOL)}4f 24,4 SRKG0148
34 GO TG 35 SRKG0149
35 IHLF=IHLF-1 SRKG0150
IS!EP=ISTEP/2 SRKG0151
h^H+H SRKG0152
GO TO 4 SRKGOlbJ
C SRKGOlS't
C RETURNS TO CALLING PRCGRAM SRKG0155
37 IHLF=12 SRKGOlSo
GO TO 3y SRKG01J7
38 IHLF=13 5RKG0153
39 CALL OUTP SRKG01b9
40 RETURN SRKG0160
END SRKG0161
00
-------
SUBROUTINE OUTP
DOUBLE PRECISI ON OERV , Y,PRMT,X,TH,XP,YP,U,8,T,h,BX,S,R,X LIM,AUX
COUBLE PRECISION A,C , XX,DCOS,OS IN,OELX
REAL I 1,12,13,14,15, 16,17,13,M
INTEGER P,RR,REG,VFG,FLAG
CIMENSION PR NT (5),Y( 1C) , AUX ( 8,10),DERY{lC),A«10,10)fC{103,VEL(8)
EQUIVALENCE (Y(8),TH),(Y(9),XP),{Y(103,YF),(Y{3),B),(Y{1),U)
EQUIVALENCE {Y<2),T),{Y{4),H),
-------
2021 FOPM/iT {/2X, 13HFRCUOE MJMBEk , <3X , 1 2HASPFCT RATIC,7X,
118HANGLE OF 01 SCHARG E ,6X , 14HRCU.MD JF F E KR CR , 7 >• , ^HX L I'« ,
214HHEAT LOSS CUEF)
URITF (P,203) FRf AS , TH D , ERR ,XLI I",E
203 FORMAT(1X,4( F9.5,12X) ,F10,5, 11X.FC.7//)
VvRITE (P,300) (I,VELiI)fI = lf3)
300 FORMAT(3(3X,3HVFL,I 1 , 1H= ,Fo. 3 ) / / )
V«RITE (P,199)
199 FaRMAT(4X,lHX,BXflHHi€X,lHB»8X,lhR,8X, 1H S , 6X , 2HFL t 7X ,
11HQ, 7X, 1HM, 5Xf 1HU,5X,1HT,4X,2HHT , 5X , 1H V , t X,
27X,2HYP,5X,3hTHD,6X,2hTM/)
150 NLINE=NLINE+1
THD=180.0*TH/3.14159
CO
FL=FL/(AS**,5)
FL=ABS(FL)
FL = FL**,5
FL=FR/FL
)*{R+H*I1 J^V*DCOS
I 2)*(R + H* I2)+Q*V*DCGS(TH)
y=M+U*V*DCOS
95
96
97
200
130
159
I7)*( R+H*I7) +T*V*UCOS( TH ) *
CELX^X-XI
DTM=2*OFLX/( L+UI )
TH=TM+DTM
IF INDIM-7H6 ,96,97
XP=X
CONTINUE
URITE (P.200) XtH»B,RtS,FL,Q,H,U,T,HT,V,XP,Y
FORMAT ( /( IX, 03,3), 51 l>,F'j,3)»2{lX,iJ3.3) ,lX,F
XI=X
LI=U
IF(ABSIM-RM I H.25) 1£C,18U,153
IF (U-V*DCOS (THS ) 16 0,160, 159
IFJU-0.02) 1£0,160,17C
PtTHD,
i».l,2X
T.M
,Od
CUT^O J37
UUTPOJ3o
OUT P.JJ39
JUTPJ041
•JUTP0042
OUTP J043
UUTPO-J45
OJTP Jo-^o
•JIJTP D-J47
GUTP0049
GDTPOO'32
.'JUTP0053
OUT POO i^
OUTPOJiid
uUTPOJbY
,1-UTPGJbo
J'JTP j J3V
OJTPJJ60
GUTPC061
G JTPOJ&2
uUT P j
•5UTPOJ&7
•JUTPJOo:1.
OUTPOJby
UUTP0070
OUTO j071
OUTPOJ7Z
-------
158 URITE (P,400) GUTPOQ73
400 FQRMAT{/10Xf43HMOMEMLM HAS EXCEEDED BOUNDS RUN TERMINATES) OUTP0074
GO TO 190 OUTP0075
160 kRITE IP ,401) OUTP0076
401 FORMAT*/10X.43HJET VELOCITY HAS BECOME TCC SfMLL RUN TERKI/gATES) OUTP0077
190 PRMT(5)=1.0 OUTPOJ7d
170 RETURN OUTP0079
END OOTPOJ80
oo
-------
SUBROUTINE CROSS CRQSOJ01
DOUBLE PRECISION 06 R V , Y , PRHT , X , Th ,/ P , Y P , I, B ,T , h , BX ,5 ,P , X L I M, AUX , XX CRUS0002
COUBLE PRECISION AUN , CUN , OCl) S , OS IN , DEL X CROS0303
RCAL T 1,12, 13,14,15, 16,17, 13 ,M CKOSOOO'+
INTEGER PtRR ,REG, VEG <, FLAG CRGSOOG5
CIMENSIQIN! FRMT(5) ,Y(1C),AUX(6, 10 },ULRY( lC),VEL(bJ CROSJJ06
DIMENSION AUN(10,10) ,CUN(10) CROSOJQ7
EUUI VALENCE (Y(8),THj , (Y(9),XP), (Y( 10) ,YP) ,(Y(3) ,B) ,(Y( 1),U) CROS0008
EQUIVALENCE (Y(2),TJ7(Y(4>,H),(Y{?),BX),{Y(6),S),(V(5),K) CKUS 0009
COHdUH OERY, VTPRMT , X , XL 1 M ,A J X , XX ,AUN ,CUK ,!3EL X CRCSOO 10
COMMON EPS8,G,OV,1ER ,HR,WI,VI , R W I, M , FK , AS , E CROSOOL1
COMMOM THIO, TH I ,ERR, V EL, 11,1 2, 13 ,!'-+,! 5 , 16,17 , lo,EPS, V,i3XS,UT ,FI CHOS0012
COMMON P,Kk, NPAGE,NL INE, REG, ND IM , VF&, F LAG , I H Lf , S TE P, TM , X I CRC1S001J
A=XP*VEL(4) CROS0014
A=A-VEL(5) CROS0015
3 )*VEL(4) CRCS0016
CRUS0017
A=VEL(3)*A CROSOOld
A=VFL(2)*EXP(-A) CRQS0019
DV = V*A CROS0020
V-VEL(1)+A CROS0021
RETURN CROS0022
r CR OS 00 2 3
-------
50
101
10100
10101
10102
10103
10104
1011
1C12
SUBROUTINE FCT
DOUBLE PRECISION OERY ,Y, PRMT , X ,TH, XP , YP , U, B, T , h, B X ,S ,R , XLI M, AUX, XX
DOUBLE PRECISION A ,C , CCOS , OS I N ,D ELX
REAL II, 12, 13, 14,15, It, 17,18, M
INTEGER P,RR, REG, VFG, FLAG
DIMENSION PRMT (5},Y(iC),AUX(8,10),OERYIlC),A(lG,10),CUO),veL(d»
EQUIVALENCE {Y(3J,THJ,(Y(9),XP),lY(10),yp),(Y(3),B),(YIR,WI,VI ,RMI,M,FR,*S,E
THIO , TH I ,ERK , VEL , 1 1 , 1 2, I 3 ,1^ , I 5 , I 6 , 1 7 , I 8, E PS, V, BX S, DT , F 1
P, RR, NPAGE,NL INE, REG, NDIM , VFG , FLAG, 1HLF , STE P, TM , XI
= 1,10
EPSB=EPS
ARG=5*T*H/( U*U*FR*FR J
ARG=ARG/(AS**,5J
IF (ARG) 10100,10 100, 10 101
ARG=C.
IF(ARG-100. U0103, 101C3, 10102
VsIR = C.
GO TO 10104
*IR=EXP(-ARG)
GO TO (1011,1012), FLAG
G=0.
BXS=BX
ex=o.
GO TC 1013
G=FI*T
FCT.NOOJ1
FCTNOJ02
FCTN0003
FCTNOUOt
FCTN0005
FCTNOU06
FCTN0007
FCTNGOOU
FCTNOOOV
FCTKOOIO
FCTU0011
FCTN0012
FCTN0013
FCTN0014
FCTN0015
FCTNJ016
FCTN0017
FCTN0013
FCTN0019
FCTN0020
FCTN0021
FCTN0022
FCTN0023
FCTN0024
FCTN0025
FCTN0026
FCTN0027
FCTN0028
FCTN0029
FCTN0030
FCTN0031
FCTN0032
FCTN0033
FCTN0034
FCTN0035
FCTN0036
-------
HI=WI*WIR FCTN0037
1013 GO TO (102,103,104,1C5), KEG FCTNQ033
102 kl='rtl*{ I 1-12) FCTN0039
VI=V1*U1-I2 ) FCTNU-J40
GO TO 11C FCTN0041
103 WI=WI*(I1-I2) FCTN0042
VI=VI*Il/2 FCTNOQ43
GO TO 110 FCTN0044
104 WI=WI*U/2 FCTH0045
VI=VI*(I1-I2) FCTNG346
GO TO lio FCTN0047
105 U=fcI*Il/2 FCTN0048
VI = VI*Il/2 FCTi'g0049
110 H1=R+H*J1 FCTNJJ30
B1=S+B*I1 FCTNOJ51
/»( 1, 1) = H1*B1 FCTNJ052
A(1,2)=U*I1*E1 FCTNOU53
A(l,3)=U*8i FCT.NI0034
^(l,^)=u*hl FCTNOO-J5
C( 1)=WI*81-VI*H1-U*I 1*H1*BX FC Fi\J056
H3=R+H*I3 FCTN0057
E3=S+B*I3 FCTNQ058
H2=R+H*I 2 FGTr
-------
00
00
1
2
3
4
5
6
120
C(
GO
A(
At
A(
At
GO
A(
GO
At
GO
b (
At
GU
10
At
At
At
C(
CO
•>.
5
5
5
5
5
5
5
4
=
4
4
*
4
)=-
TO
,1)
»2)
,3)
,5)
TC
,4)
TO
,3)
TO
» 4)
,3)
TO
1.-
,1)
I 3)
,43
E*G*83-U*G*I 7*h?*BX
(
1
f
2,3,4), REG
= U
= G*
13
=G
=
5
=
5
=
5
=
=
6
I
=
-
=
F
1
1
1
1
I
«
•
.
•
2/
U
p
o
U
U
)=U*(
TO
<
1
CALL CRO
*
*
*
*(R/2+H*I 3)
0
0
11
{{S*H+R*B^*I.2+R*S*I
u*S*I8+G*J"S*t-*I3 + R*
*U*R:*>I8+G*t-R:SR/2-HH*K
V
i* iO'll'fl' l^^"L'«N'»\'i
I*R-WI*S) *I 2/Il-G*t
20,131), VFG
SS
8)
3 )
DERY{9)=DSIN(TH}
CERY(10)=DCOS(TH)
VCT=V*DERY(1C)
DT=U*U*{ S+B*12 )*(R+H*I2)+2*U*VCT*(S + S
1+VCT*VCT*(R+H)*(S + B )
CT=-VST*(WI*(S+B*I1)-VI*{K+H*I1))/OT
DERY(8)=DT
1 J* ( R -U-*! 1 )
FCTiMOj73
FCTN'J074
FCTNJO 75
FCTN0076
FC TNJ077
FC1 ''-J0u7o
FCT.N0031
FCTN J0d2
FCTN0033
FCI'.'MOObS
FCTI',00 Jo
FCTi\iOJo7
FCTN J038
FCT,MOu-39
FCTiM J.J9J
FCTNOjyi
FCTN J;J92
FCTN0093
FCTN0094
FCTN 00 95
FCTNJJ-S6
EO=S+B
FCTNOJ98
FCTNG01-39
FCTN 01 30
FCTN0101
FCTiM01Q2
FCTN0103
FCTN0104
FCTi\i0105
FCTN0106
FCTN01U7
FCTNJ10B
-------
03
=A (1,2 )+VCT*BO
=A t 1, 3) + VCT*rtO
= A ( 1, 4)+VCT*HO
(1 )-VCT*HO*BX+BC*HO*UT
= A( 2, 1 H2*VCT*H1*81
=A (2,2 )+VCT*(2*L*I 1*81 +VCT *f30 )
=A ( 2, 3)+VCT*(2*L*bl+VCT*BOJ
=A12,4 ) + VCT*(2*L*Hl+VCT*HO 3
(2 )-VCT*(2*U*I 1*HL+VCT*HO) *BX + 2* ( U
=A (3,2 )+G*VCT*I3*B3
^A( 3,
= A (3,
=A(3,5)+Fl*VCT*H3*B3
A( 1,2)
A( 1, 3)
A(1,4}
C(1)=G
A ( 2 , 1)
A(2,2)
A{2,3)
A(2,H)
G ( 2 ) =G
1+VCT*(
A(3,2)
A(3,4)
A(3,5)
GO TC (7,<3,3,9),REG
7 A{5,1)=A(5,1)+VCT
G(5)= (U+VGT)*DT
GO TG 9
8 A(4,1}=A(4,1)+VGT*I1*(R
A (4, 2) = A (4, 2 ) + S*U*VCT*U 1-I 1 / I 2 )
A{4,3)=A14,3 )+S*U*VCT*I8
A14,4)=A (4, 4 J+R*U*VCT*I8
C(4)=C(4)-R*U*VCT*(I 1-I1/I2)*BX
1-K (U*( 2* I 1-1 2/1 D + VC T )*( S*m-R*B )+R*S*U*I 8J*DT
9 GQ TC (130,131 ),FLAG
130 A(l,2)=A(l,2)+(U*I1*H1+VCT*HOJ
A(2,2)=A(2,2 )-KU:*;U*12:«H2+2*U-VGT*[ 1*H 1+ VC T*VC T*H 0 )
GO TG ( 10,10,10,11),REG
10 A (4, 2)= A (4, 2) + R*U*VC1*U 1-12/11)
11 CALL SGELG(5)
EPSO^C (2 )
BX-BX'j
FLAG=2
G=FI*T
^I=W I*W I R
L+VCT *30*MO )*DT
FCTN0109
FCTN0110
FCTN0111
FCTN0112
FCTN0113
FCTN0114
FCTN0115
FCTN0116
FGTN0117
FCTN0118
FCTN0119
FCTN0120
FCTN0121
FCTN0122
FCTN0123
FCTN0124
ECTN0125
FCTN0126
FCTN0127
FCTN0128
FCTN0129
FCTNU130
FCFN0131
FCTN0132
FCTN0133
FCTNU134
FGTNOUo
FCTN0137
FCTNOIBb
FGTN0139
FCTNOl^O
FCTN0141
FCTN0142
FCTNOl-tj
FCTNOL44
-------
GO TO 110 FCTNU145
131 D8=8X-EPSB FCTN0146
A(6,1)=DB*2*U*B*H2*I6 FCTN0147
A( 6,2)=DB*U*U*B*I2*I6 FCTN0148
/(6,3)=D6*U*U*B*I6 FCTN0149
A( 6, 6)=U*U*B*H2*I6 FCTN0150
C(6)=G*HG-U*l*I6*DB*h2*BX FC'i NOlbl
GO TO (132,133),VFG FCTN0152
132 A(6, 1)=A(6, 1) + 2*B*VCT*I5*H2#D8 FCT,N015"3
A(6,2)=A{6,2)+VCT*B*{2*U*I5*I1+VCT)*D6 FCTN0154
A(6,3)=A(6,3)*VCT*B* (2*U*I5+VCT)*OB FCTN0155
A(6,6)=A{6,6)+VCT*8*(2*U*I5*Hl+VCT*HO) FCTNOlSo
C(6)=C(6 )-D8*VCT*{2*U*I5*HH-VCT*hO)*BX*-2:*UB*d*{U*I5:!-Hl+VCT FCINOl j?
1*HOJ*OT FCTN0158
133 CALL SGELG15) FCTN0159
DERY(1)=C(1J FCTN0160
DERY(2J=C(5) FCTNOlol
CERY(3) = BX FCTN0162
OERY(A)=C(2) FCTN0163
CERY(5)=C(3) FCTNOltot
DERY(6) = C{'t} FCTN0165
DERY(7)=C(6) FCTNOlLo
PETURN FCTNOlo?
FCTNOloS
-------
51
52
53
SUBROUTINE SGELG(M)
DOUBLE PRECISION AD » A I , XX,DA tiS , OCOS , OSI N ,SAV E ,S , R , XL IM, AUX
COUBLE PRECISION DERY,Y,PRMT,X,TH,XP,YP»L,B , T,H,BX
DOUBLE PRECISION A,P IV ,TB,TOL,PI VI,DELX
REAL MUiM
INTEGER P,RR,REG,VFG,FLAG
DIMENSION PRHT{5),Y(1C) ,AUX(8,10),UERY{1Oi,AI (100) ,R (10) ,VEL(8J
DIMENSION A(100),AD(1C,10)
DIMENSION SAVE(IO)
EQUIVALENCE (AI(1),AC(1, 1))
COMMON DERY,Y,PRMT,X,>LIM,AUX,XX,A I,R,DELX
COMMON EPS6,G,DV,IER,HR,WI , VI ,R VI, MUN ,F R , AS ,E
COMMON THID,THI ,ERR,V£L, II,I 2,13,14,15,I 6,17,18,EPS,V,
COMMON P,R^,NPAGE,NLINE,REG,NDIM,VFGfFLAG,IHLF,STEP,TM,XI
EPSS=.1E-16
SCALE THE MATRIX
CO 54 1=1 ,H
AMAX=R(I)
DO 52 J=1,M
DA=DA3S{AC( I ,J))
IF(DA-AM AX)52,52,51
AMAX=DA
CONTINUE
DO 53 J=1,M
PQ{ I ,J )=AD( I,J J/AMAX
CONTINUE
Rt I )=R(I )/AMAX
CONTINUE
SAVE THE Y EQUATION -PCv*
999
DO 999 1=1, MM
SAVE ( I)=AD(MM, I)
CONTINUE
CONVERT MATRIX TO
CJLLKN FORM
SGEL0001
SGEL0002
SGEL0003
SGEL0004
SGELOU05
SGELQ006
SGEL0007
SGEL0008
SGEL0009
SGEL0010
SGEL0011
SGELG012
SGEL0013
SGEL0014
SGEL0015
SGELuOlo
SGEL0017
SGELOJ18
SGEL0019
SGELOJ20
SGELOJ2L
SGEL0022
SGEL0023
SGEL0024
SGEL0025
SGELOJ26
SGEL0027
SGELOJ28
SGELU029
SGEL0030
SGEL0031
SGEL0032
SGEL0033
SGELOQ34
SGEL0035
SGELOOSo
-------
1000
1100
c
c
C
C
c
c
c
c
NM=0
DO 1100 K = l,f
CO 1000 L=l,l*
IJ=IJ+1
NM=NM+1
A(IJ)=Al(NM)
NM=NM+10-M
IF(M}23,23, 1
SEARCH FOR GREATEST ELEMENT IN MATRIX
IER=0
PIV=C.DO
DO .3 L=1,«M
T8 = DA8S(A(L) )
IF ITS-PIV)3,3,2
2 PIV=TB
CONTINUE
TOL=EPSS*PIV
A(I) IS PIVOT
ELEMENT. PIV CONTAINS ThE ABSOLUTE VALUE OF All)
START ELIMINATION LOCF
LST=1
CO 17 K=ltM
TEST ON SINGULARITY
IF{PIV)23,23,4
IF( IER)7,5, 7
IF(PIV-TCL)6,6,7
IER = K-1
PIVI = 1.00/A( I)
J={ 1-1 )/M
I=I-J*M-K
SG EL 00 3 7
SUELGG3c>
SGELJ039
SGELOU40
SGEL0341
SGELOOA^I
SGELQ043
SGELOJ4o
SGELJJ47
SGEL0048
SGELOJ49
SGEL005J
SGELOO:>1
SGELOJ32
SGELQ053
SGELOOD4
SGEL J055
SGEL JJ56
SGELOJ57
SGEL0058
SGEL 005V
SGELJJ6J
SGELOJol
SGFLJ062
SGLLOJ63
SGELOQ04
S GEL J Jo 3
SGEL 0066
SGEL006?
SGEL0069
SGELOJ70
SGEL JO 71
SCEL0072
-------
c
c
c
J=J-H-K
I + K IS ROW-INDEX, J+ K CCLUM.M-I N-0 EX Oh PIVbT ELMENT
c
c
c
c
c
c
c
c
c
c
PIVOT ROW RECUCTION
DU 8 L = K,NM,!*
LL=L+I
Tb=PIVI*R(LL)
R(LL)=R(L)
3 P(L)=TB
IS ELIMINATION TERMINATtD
IF(K-M)9,18,18
COLUMN INTERCHANGE IN MATRIX A
9 LEND=LST+M-K
IF{ J)12, II, 10
10 II=J*M
DO 11 L=LST,LEND
TB=A(L)
ROW INTERCHANGE IN RIGHT HAND S I Dt R
A(L)=A(L.LJ
11 MLL) = TB
ROV< INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A
12 DO 13 L=LST,^M,M
LL=L+I
T3 = PIVI*A(LL )
A(LL)=A(L)
13 A(L)=TB
SAVE COLUMN INTERCHANGE INFORMATION
A(LST)=J
ELEMENT REDUCTION AND NEXT PIVOT SEARCH
PIV=O.DO
LST=LST+1
SGEL0073
SGEL0074
SGELJU75
SGEL 0076
SGEL0077
SGEL007d
SGEL0079
SGEL0080
SGELOJ81
SGEL0062
SGEL0083
SGEL0084
SGEL0035
SGEL0086
SGEL0037
SGELOOdS
SGEL0039
SGEL 0090
SGEL0091
SGEL0092
SGEL0093
SGEL0094
SGELOQ95
SGEL0096
SGEL0097
SGEL0098
SGEL0099
SGEL0100
SGEL0101
SGEL0102
SGELCM.03
SGEL0104-
SGEL0105
SGEL0106
SGEL0107
SGEL01C8
-------
J=0 SGELOUv
DO 16 II=LST,LEND SbtLxDilj
PI VI —At I I) SGELJlil
IST = IH-M SGEL3112
J=J+1 SGELOliJ
CO 15 L=IST,MM,M SGELOllt
LL^L-J SGEL0115)
A(L)=A(L) *PIVI*A(LL) SGELOllo
TB=DA8S(A(L)) SGELOL17
1F(TB-PI V)15T15,14 SGELJil,:
14 PIV = TB SGELJil'v
I=L SGEL012J
15 CONTINUE SGEL0121
CO 16 L=K,NH,M SGEL0122
LL=L+J SGEL0123
16 R(LL) = R(LL) +PI VI*RIL ) SGELOl^^t
17 LST=LST + M SGEL012-J
C END Of ELIMINATION LCCP SGEL0126
C SGtLOl^/
C SGEL0123
C BACK SUBSTITUTION ANC BALK I INTERCHANGE SGEL012^
18 IFlM-l)23,22,19 SGELJliJ
19 IST=MM+M SGELOl^l
LST=M+1 SGEL0132
CO 21 1 = 2,M SGEL0133
II=LST-I SGEL0134
IST=IST-LST SGELOlbi?
L=IST-M SGELOlio
L^A(L1+.5DO SGEL0137
CO 21 J=II,NM, M SGEL0133
7B=R(J) SGELJ139
LL=J SGELOltO
DO-20 K=IST,I*'M,M SGtLG141
LL^LL+1 SGEL0142
20 TB = T6-A(K)*R{LL) SGEL01'+3
SGEL01
-------
R(J J-K(K ) SGEL 0145
21 R(K)=TB SGtL0145
22 GO TO (223t221)i t-LAG SGtLJl^?
221 M=iv-H SGELG-1-fci
CG 222 I=1,M SGEL3149
R(MKi)=R(>lM)-SAVc(I } *R ( I ) SGELOlbo
2?.2 CGMINUE SGfcLul-31
R(MM)-'^( MM) /SAV£{ 4M) SGEL 01 32
223 CU 24 1=1, 1GC SGtLOl^S
A(I)=0. SGEL 0154
24 CUMTINUE SGELOIS^
RETURN SGEL01S6
C 3GELJ1S7
C SGEL0158
C EKRGR RETURN SGEL0159
23 IER=-1 SGEL0160
RETURN SGEL0161
END SGELulfc2
-------
Input:
Output:
1
6.0
0.0
0.6
0.0
0.0
0.0
90.0
6.0
500.0
0.0
0.01
0.0
1.0
0.0
0.0
BUCYAM JET CALCULATICNS
PAGE 1
FPCUCE NUMBER
6.CCCOO
ASPECT RATIO
0.60000
ANGLE OF DISCHARGE
90.00000
ROUNDOFF EkRCR
O.C10QC
XLIH
500.COOOO
HEAT LOSi CUEh
.0
VELl =
0.0
VEL2=
0.0
VEL3 =
c.o
VEL4=
0.0
VEL5= 0.0
O.d
VEL7= 0.0
0.0
FL
HT
XP
YP
THD
.0
.1310
.2620
.5240
.7660
.9160
. USD
.1310
.1570
.2100
.2360
.2490
.2620
.2750
.2880
.3150
.3470
.4190
.4120
.5240
.6290
.734D
.8390
. 1C5D
. 1260
. 1470
. iteo
.2100
.2520
Cl
01
Cl
Cl
01
C2
02
C2
C2
02
C2
C2
02
02
02
C2
C2
02
C2
02
02
C2
03
C3
03
C3
03
C3
.1000-03
.1900 00
.3780 00
.7240 00
.1040 01
.1190 01
.1340 01
.1610 01
.1850 01
.2210 01
.2240 01
.2270 01
.2280 01
.2280 01
.2270 01
.2250 01
.2200 01
.2150 01
.2090 01
.2040 01
.2020 01
.2110 01
.1950 01
.1560 01
.1250 01
.1020 01
.3590 00
.6410 00
.5060 00
.1000-03
.5850
.1160
.2390
.3770
.4500
.5270
.6920
.8730
.12SO
.1530
.1660
.1800
.I960
.2110
.2460
.3240
.4150
.5160
.62SD
.8840
.1240
.1740
.3240
.5520
.8640
.1260
.2350
.3820
00
01
01
01
01
01
01
01
C2
02
02
02
02
02
02
02
02
02
02
02
03
03
03
C3
03
04
04
04
.7750 CO
.5210 00
.4500 00
.2180 00
.2930-C1
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.1290
.1520
.1370
.1440
.1530
.1520
.1480
.1300
.1000
.1950
.1030
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
Cl
Cl
Cl
Cl
Cl
Cl
Cl
Cl
Cl
CO
CO
.528E
.124E
.890E
.664E
.574E
.517E
.465E
.391E
.342E
.283E
.267E
.258E
.250E
.244E
.238E
.227E
.211E
.197E
.186E
.177E
.151E
.107E
.868E
.658E
.543E
.470E
.413E
.349E
.304E
03
02
01
01
01
01
01
01
01
01
01
01
01
01
01
01
01
01
01
01
01
01
00
00
00
00
00
00
00
.100E
.109E
.USE
.133E
.161E
. 177E
.L9oE
.232E
.265E
.321E
.347E
.360E
.373E
.387E
.400E
.'*26E
.476E
.523E
.568E
.609E
.683E
.707E
.715E
-720fc
.723E
.723E
.722E
.720E
.716E
01 1.C14 l.COO l.COO l.uOC 0.0
01 1.014 1.C05 O.S71 l.OOG 0.0
01 1.014 1.C04 0.940 1.000 0.0
01 1.014 1.C06 0.684 l.OOC 0.0
01 1.014 1.C07 0.828 I. 000 0.0
01 O.S95 0.931 0.759 C.984 O.O
01 0.997 0.8t7 0.692 0.986 0.0
01 l.CJO 0.726 0.597 0.988 0.0
01 1.002 0.645 0.535 C.99C 0.0
01 1.002 C.539 0*457 0.990 0.0
01 O.S99 0.492 0.424 C.987 0.0
01 0.996 0.471 0.409 0.985 0. C
01 0.996 0.449 0.394 0.985 0.0
01 O.SS7 0.429 0.381 C.935 0.0
01 0.957 0.411 0.368 0.985 0.0
01 O.SS7 0.380 0.346 C.985 0.0
01 0.997 0.329 0.310 C.985 0.0
01 0.997 0.290 0.282 0.985 0.0
01 0.997 0.260 0.260 0.935 0.0
01 0.997 0.235 0.242 0.985 0.0
01 0.994 C.189 0.215 0.984 0.0
01 0.995 0.134 0.208 0.984 0.0
01 0.994 0.104 0.206 0.983 0.0
01 O.SS1 O.C7Q 0.204 C.981 0.0
01 0.938 O.C52 0.202 0.978 0.0
01 0.985 O.C40 0.202 0.974 0.0
01 0.981 0.033 0.201 C.971 0.0
01 0,973 0.024 0.200 0.963 0.0
01 0.965 O.C18 0.200 0.956 0.0
.0
.1310
,2b2D
.5240
.7860
.9130
. 1050
.1310
.1570
.2100
.2360
.2490
.2o2D
.2750
.2880
.3150
. 367U
.4190
.4720
.5240
.6290
.7340
.8390
.1050
. 1260
.1470
. 1630
.2100
.2520
01
01
01
01
01
02
02
02
02
02
02
02
02
02
02
02
02
02
02
02
02
02
03
03
03
03
03
03
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
90.0
90.0
90.0
90.0
90. 0
9o.O
90. 0
90.0
90.0
90.0
90.0
90.0
90.0
90.0
90.0
90.0
90.0
90.0
9u.O
90. 0
90.0
90.0
90.0
90.0
90.0
90.0
90.0
90.0
90.0
.0
.13 IE 01
.26 1C 01
.5226 01
. 73iE 01
.9! 8E 01
.107E 02
.14CE 02
.17£E 02
.267E 02
.317E 02
.345E 02
.373E 02
.t03E 02
.434E 02
.50 IE 02
.648E 02
. 818E 02
.101E 03
.122E 03
.171E 03
.236E 03
. 325E 03
.565E 03
.909E 03
.136E 04
. 194E 04
.342E 04
.542E 04
JET VELOCITY HAS BECOME TOO SMALL RUN TERMINATES
Appendix
Input and Output
-------
SELECTED WATER
RESOURCES ABSTRACTS
INPUT TRANSACTION FORM
1. Report No.
3. Accession No.
4. Title
A User's Manual for Three-Dimensional Heated
Surface Discharge Computations
7. Author(s)
K.D. Stolzenbach, E. Eric Adams, and Donald R.F. Harleman
9. Organization
Department of Civil Engineering, Massachusetts Institute
of Technology, Cambridge, Massachusetts 02139
12. Sponsoring Organization
Environmental Protection Agency
IS. Supplementary Notes
Environmental Protection Agency Report
Number EPA-R2-73-133* January 1973.
5. Report Date
6.
8. Performing Organization
Report No.
10. Project No.
11. Contract/Grant No.
16130 DJU
13. Type of Report and
Period Covered
-I
16. Abstract
The temperature distribution induced in an ambient body of water by a surface dis-
charge of heated condenser cooling water must be determined for evaluation of thermal
effects upon the natural environment, for prevention of recirculation of the heated
discharge into the cooling water intake, for improved design of laboratory scale models
and for insuring that discharge configurations meet legal temperature regulations. This
report presents a review of the theoretical background for a three-dimensional tempera-
ture prediction model, a detailed discussion of the computer program and a case study
illustrating the procedure for optimizing the design of a surface discharge channel.
Flow chart, program listing and a sample of the input and output data are given in the
appendices. The model presented here includes modifications of the report by Keith D.
Stolzenbach and Donald R. F. Harleman, published in February 1971 entitled, "An
Analytical and Experimental Investigation of Surface Discharges of Heated Water."
17a. Descriptors
Waste heat disposal, heated surface discharge, turbulent buoyant jets,
temperature prediction, thermal pollution.
17b. Identifiers
17c. COWRR Field & Group Q5G
18. Availability
19. Security Class.
(Report)
20. Security Class.
(Page)
21. No. of
Pages
22. Price
Send To:
WATER RESOURCES SCIENTIFIC I N FOR M AT ION CENTER
US DEPARTMENT OF THE INTERIOR
WASHINGTON, D C 20240
Abstractor
Massachusetts TnQtitMtft of Technology
WRSICI02(REV JUNEI971)
GP-0 9 13.281
ftU.S. GOVERNMENT PRINTING OFFICE: 1973 514-151/148 1-3
------- |