United States
Environmental Protection
Agency
Corvallis Environmental
Research Laboratory
Corvallis, Oregon 97330
PREDICTION OF INITIAL MIXING FOR
MUNICIPAL OCEAN DISCHARGES
by A. M. Teeter and D. J. Baumgartner
CERL - 043
-------
PREDICTION OF INITIAL MIXING FOR
MUNICIPAL OCEAN DISCHARGES
by A. M. Teeter and D. J. Baumgartner
CERL - 043
Marine Systems Division
Corvallis Environmental Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
200 S.W. 35th Street
Corvallis, Oregon 97330
-------
PREDICTION OF INITIAL MIXING FOR
MUNICIPAL OCEAN DISCHARGES
by
A. M. Teeter and D. J. Baumgartner
Marine Systems Division
Corvallis Environmental Research Laboratory
Corvallis, Oregon 97330
CERL PUBLICATION 043
May, 1979
CORVALLIS ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
CORVALLIS, OREGON 97330
-------
CONTENTS
Page
Abbreviations, Symbols, etc i"11
Acknowledgments . v
1. Introduction
Mixing Zone Concepts . . 1
Ocean Discharges 2
Buoyant Plume Models 4
2. Recommendations
Methods 7
Appropriate Conditions 7
Mixing Zone Specification.
3. Summary and Conclusions
General. . . 13
Simplified Results for Plumes in Stratified Ambients .... 15
Plumes in a Stagnant Environment 16
Plumes in a Flowing Environment . 17
Merging Plumes in a Stagnant Environment 19
Merging Plumes in a Flowing Environment 20
Non-Uniform Stratification 22
Plumes in Uniform Density Ambients .... 22
4. PLUME Model
Theoretical Development 26
Model Description 29
Example Input, Output and Model Listing. . . 30
5. OUTPLM Model
Theoretical Development and Model Description 41
Example Input, Output and Model Listing 42
6. DKHPLM Model
Theoretical Development. ....... 51
Model Description. . 54
Example Input, Output and Model Listing . . 55
References 85
ii
-------
LIST OF ABBREVIATIONS, SYMBOLS, ETC.
B - Radius or characteristic half-width of the plume.
C - Concentration by volume, = 1/Sa.
Ca - Concentration in ambient.
Ce - Concentration in effluent.
Cf - Concentration after initial dilution.
D - Diameter of discharge port.
E - Entrainment.
1/2
Fn - Froude number, = U./(g D A0/po)
J
g - Gravity.
Q - Volume rate of discharge.
q - Volume rate of discharge per unit length of diffuser.
r - Radial distance from the plume's centerline.
s - Distance along plume axis.
Sa - Flux averaged dilution or the total plume efflux divided by the volume
of the discharge.
sg - Length of the zone of flow establishment.
T - Temperature.
T' - Fluctuating component of T.
U - Ambient current speed.
U. - Initial velocity of the jet.
3
u' - Fluctuating component of V.
V - Velocity at some point in the plume in the s direction,
v - Velocity normal to the plume's axis.
i i i
-------
v' - Fluctuating component of v.
x - Horizontal distance in the direction of a jet discharge and ambient
currents or the horizontal distance parallel to a diffuser.
y - Horizontal distance 90° to x.
z - Vertical distance from the discharge depth.
2 - Depth of water.
A - Normalized density disparity, = (p^ - p)/(pQ " Pd).
Ao ~ Initial density disparity of the waste, = p0 - p^.
Ah - Maximum height of rise or height of rise to the trapping level.
p - Density of the plume at some point.
p0 - Ambient density at the level of the discharge.
Pd - Density of the discharge.
Pw - Ambient density at some level.
dp/dz - Absolute value of the ambient density gradient.
1 m3 = 35.31 ft3 = 264 gallons.
iv
-------
ACKNOWLEDGMENTS
The assistance and cooperation of Drs. Lorin Davis, Larry Winiarski, and
Safa Shirazi, and Mr. Walter Frick of the Ecosystem Modeling and Assessment
Branch of the U.S. Environmental Protection Agency's Corvallis Environmental
Research Laboratory are gratefully acknowledged. A similar acknowledgment is
due Bill McDougal of the Marine and Freshwater Branch.
v
-------
SECTION. 1
INTRODUCTION
Initial dilution is the rapid turbulent mixing between wastewater dis-
charged at depth in the ocean and seawater. It takes place close to the point
of release, resulting from physical interactions between these components. The
understanding of this process can be used to predict the initial dilution
attained under given conditions. Such information enables comparison of waste-
water constituent concentrations, subsequent to initial dilution, to water
quality standards and thus can assist in evaluating the discharge's effect on
the marine environment.
The purpose of this repprt is to present procedures for calculating the
initial dilution and describing the zone of initial dilution surrounding a
discharge site. These calculations would be required under regulations pro-
posed by the U.S. Environmental Protection Agency on April 25, 1978 (1) imple-
menting section 301(h) of the Clean Water Act (2), draft regulation implement-
ing section 403(c), and provisions of the State of California's Ocean Plan (3).
MIXING ZONE CONCEPTS
When wastewater is injected into the ocean, it has buoyancy and momentum.
This energy is used in an initial mixing or dilution process. The wastewater
discharge forms a plume which rises and entrains ambient water and is thus
diluted. The mixing process slows markedly when plumes reach a position of
equal buoyancy and momentum with the ambient (no relative energy) or the
surface. Then the dilution process becomes more dependent on ambient oceanic
1
-------
processes. Thus, the mixing zone is a region where rapid mixing takes place
between waste and ambient waters. It can be thought of either as a concentra-
tion isopleth or as the zone where initial mixing takes place. It is not a
"treatment for wastes but a transition region (4). By insuring that water
quality standards are not violated outside this zone, impacts on the marine
environment from the waste discharge should be controlled. While the dilution
of some pollutants may have little effect on their impacts, biological uptake,
concentration, and deposition, and physical processes, like aggregation and
sedimentation, are concentration dependent.
OCEAN DISCHARGES
Ocean discharges of municipal wastes have a wide variety of physical
characteristics which can affect initial dilution. Although single port out-
falls are still used, many outfalls now are multi-port designs to optimize the
dilution and vertical rise attained by their plumes. The numbers of ports and
their orientation can form complex patterns reflecting a great deal of design
effort. The discharged material is almost uniform, with physical characteris-
tics similar to those of fresh water. The marine environment, however, varies
greatly along the miles of our country's shoreline and the physical conditions
at the discharge site have a great effect on initial dilution. The ambient
receiving waters are of a higher density than the waste discharge and are
typically vertically stratified. Currents are usually present but may have
periodicity caused by tidal and other forces.
Discharges initially have momentum and buoyancy as they are injected into
1 /2
the environment. Densimetric Froude numbers (Fn = U./(g D A0/po) ) describe
J
the ratio between these two forces. The higher the initial Froude number, the
more closely the resulting plume resembles a momentum or forced jet. The
2
-------
smaller the number, the more the plume resembles a purely buoyant element.
Because the density difference between the waste and the ambient varies only
slightly and jet velocity is bounded by practical considerations, Froude numb-
ers for ocean discharges generally vary between 0 and about 30. The ends of
this range are characteristic of single-port and multi-port diffusers, respec-
tively. In this range of Froude numbers, buoyancy is likely to dominate the
initial mixing process, making volume rate of discharge, currents, ambient
stratification, and possibly water depth, important. Especially in a current,
buoyancy dominates since the plume's momentum becomes insignificant compared to
the momentum entrained from the ambient.
Existing coastal discharge sites vary in depth up to about 75 meters. The
water columns at these sites are typically stratified vertically by differences
in temperature and/or salinity. This leads to density stratification. The
amount or degree of stratification varies widely in time and space. Proximity
of large sources of fresh water, surface heating, advection, and wind mixing
all can influence density stratification. Salt wedges can form near the mouth
of estuaries with large fresh-water flows. Advection and nonlinear interac-
tions between wind-induced surface turbulence and buoyancy gradients also can
produce pycnoclines, areas of strong density gradients or discontinuities.
Strong pycnoclines can arrest the vertical motion of buoyant plumes, ending the
initial mixing process. Winds of sufficient duration can erode the density
gradient by the turbulence they produce, leaving the ambient well mixed and
allowing buoyant plumes to rise to the surface.
Currents are very complex in coastal areas. These are areas of high
energy received from a number of processes. Currents at discharge sites have
tidal components which can be strong. Oceanic currents impinge on the coast
and intensify on the western edges of the oceans (the Florida Current for
3
-------
instance). Local forces such as wind shear and waves also generate currents.
Instantaneous current values (time scale of minutes), rather than long-term
(time scale of hours, days, etc.) net currents, are important in the initial
dilution process.
BUOYANT PLUME MODELS
Ocean outfalls can be modeled by properly scaled hydraulic models, or by
mathematical models. The methods described here include only the latter,
although in cases of complex geometry or special conditions, physical or hy-
draulic modeling may be appropriate. The characteristics and behavior of
plumes and jets in other aquatic settings, the atmosphere, as well as for
ocean outfalls, have been studied by means of mathematical models. Many
excellent reviews have been published (5, 6, 7).
Early work in this field was concerned with purely buoyant convection (no
initial momentum). Rouse et aL (8) studied the buoyant plume above a contin-
uous heat source in a uniform, stagnant environment. Equations of mass contin-
uity, momentum, and energy were integrated by assuming Gaussian approximations
for the lateral distributions of velocity and temperature across the plume and
with experimentally determined spreading constants. The results were used to
describe the distribution of velocity and temperature as functions of mass
flow and height of rise. Priestly and Ball (9) studied thermal plumes issuing
into a stratified, as well as uniform, environment. They used essentially the
same equations as Rouse et aL for mass, momentum, and energy, but their
lateral Gaussian profiles of temperature and velocity assumed a lateral length
scale.
Morton et a2. (10) studied buoyant plume convection. They introduced an
entrainment function based on local characteristic plume velocity and length
4
-------
scale to replace the spreading ratio. This approach has been adapted by many
later investigators (5). Later, Morton (11) assumed other lateral profiles of
velocity and temperature to account for the difference between eddy transport
of heat and momentum. Another coefficient was introduced and estimated from
experimental data.
Abraham (12) investigated neutrally-buoyant forced plumes, purely-buoyant
plumes, and plumes with combinations of both properties. Equations of mass,
momentum, and energy were integrated using lateral Gaussian profiles recogniz-
ing the difference in transfer rates between momentum and heat (or concentra-
tion). Entrainment was considered by integrating the continuity of mass equa-
tion with appropriate spreading coefficients. Abraham's results included
analytic expressions for plume centerline velocity and concentration as func-
tions of initial plume characteristics, height of rise, and entrainment coef-
ficients for momentum and heat. The limiting cases of neutrally-buoyant
forced jets and purely-buoyant plumes were considered.
Fan (13) reported on plumes issuing at arbitrary angles into a stagnant,
stratified environment. Velocity and buoyancy profiles were used similar to
those used by Morton (11) and included a lateral characteristic length. Solu-
tions were found by simultaneously solving six equations (continuity of mass,
vertical momentum, horizontal momentum, conservation of buoyancy, vertical
geometry, and horizontal geometry).
Fan (13) also investigated buoyant plumes in a flowing, uniform-density
ambient. A theoretical solution was given as a function of plume drag forces
and an entrainment function. The entrainment function was such that drag and
entrainment coefficients varied with plume conditions, a disadvantage in
applications beyond the range of experimental validation (14). Abraham (14)
introduced an entrainment function based on subdomains influenced by initial
5
-------
momentum and buoyancy. With entrainment coefficients defined for these two
subdomains and a single coefficient describing the drag force, the theoretical
solution was in good agreement with the experimental results of Fan.
Plumes from multi-port discharges have been studied less extensively.
Cederwall (15) studied the flow regimes of line sources for discharges in
confined and unconfined environments. Roberts (16) studied the flow regimes
of line sources in unstratified and unconfined environments. Soti1 (17)
proposed a mathematical model for slot jets or continuous line sources in
stagnant, stratified ambients. Koh and Fan (18) proposed a mathematical model
of multi-port diffusers by interfacing single jet and slot jet solutions at a
transition point. A good review of this subject is given by Davis (7).
Three computer models are described in this paper. The PLUME model by
Baumgartner et aK (19) simulates a solitary plume in a stagnant environment
of arbitrary stratification. OUTPLM by Winiarski and Frick (20) models a
single plume in a stagnant or flowing environment of arbitrary stratification.
DKHPLM by Davis (21) describes single plumes in flowing, stratified ambients
which are allowed to merge with identical adjacent plumes to simulate closely
spaced, multiple-port diffusers. The results of these models have been used
to develop analytic solutions to describe the initial dilution of ocean dis-
charges in terms of the principal quantities involved. These three computer
models and the analytic solutions resulting from them will be presented along
with suggestions and example applications.
6
-------
SECTION 2
RECOMMENDATIONS
METHODS
The buoyant-plume phase of waste dispersion can be described by a variety
of published methods. However, to evaluate the water quality impacts of
municipal ocean discharges, the methods presented here are well suited and
recommended for consistency. The analytic methods or computer models can be
used as the situation warrants or requires. When especially high resolution
is necessary, the computer models should be used. For unusual situations or
conditions, other methods will be considered.
Although these methods are considered reliable, they should be continu-
ally evaluated relative to theoretical developments and measurements, espec-
ially quality field measurements.
APPROPRIATE CONDITIONS
The dilution achieved during the initial mixing process is dependent on
ambient and discharge conditions and is, therefore, highly variable. To avoid
the inception of biological effects, the occurrence of pollutant concentra-
tions greater than water quality criteria should be small. In evaluating a
discharge's effect on water quality, therefore, the appropriate conditions to
consider are those which result in the "lowest" dilution and those which occur
at times when the environment is most sensitive. For example, minimum dilu-
tion will be predicted using a combination of maximum vertical density strati-
fication, maximum waste flow rate, and minimum currents for a particular site.
7
-------
Other situations may be more critical to water quality depending on the ambi-
ent water quality, and applicable water quality standards.
Predicting "lowest" dilution relies on the availability of quality meas-
urements or data with which to estimate absolute worst case conditions, and is
therefore not without problems. The statistical uncertainty in estimates of
absolute worst case conditions is generally great. Also there are inherent
biases to some oceanographic measurements. For example, current measuring
instruments have finite thresholds. It therefore becomes difficult to dis-
tinguish low values (which may be as high as 0.05 m/sec) from zeroes with
these data. In estimating environmental conditions, a more reliable estima-
tion can be made at the lowest 10 percentile on a cumulative frequency distri-
bution. Data on ambient density structure are not collected by remote means
so that the frequency of observations is generally low. To increase the
reliability of "worst-case" estimates, data should be evaluated not only from
the discharge site but from nearby coastal areas of similar environmental
setting.
Defining "worst-case" conditions as a combination of affecting conditions
each taken at the.worst 10 percentile on cumulative frequency distributions is
recommended. This approach allows a reliable estimation of these conditions
to be made and prevents the unlikely occurrence of more extreme conditions
from biasing the predictions. The probability of these conditions occurring
simultaneously is much less than 10 percent insuring that the predicted dilu-
tion will be exceeded most of the time. Given the sensitivity of dilution to
external conditions and the ranges of these conditions, the evaluation of data
and selection of model input conditions are critical to the proper evaluation
of water quality impacts from municipal ocean discharges.
The characteristics of the discharge needed for evaluation are:
8
-------
a) port sizes, spacings, depths and orientation,
b) flow rates through ports,
c) wastewater density, temperature and apparent salinity at the point
of discharge.
The hydraulic characteristics of diffusers should be checked and the
diffuser segment with the highest flowrate per unit diffuser length used for
subsequent initial dilution predictions. Apparent salinity refers to the
major inorganic ions of the wastewater which contribute to its density. It
has been found that a salinity, based on conductivity ratios for seawater,
gives an acceptable value. Alternately, an apparent salinity which corres-
ponds to a measured density can be used. Apparent salinities of municipal
effluent generally range from 1 to 5 °/oo.
As mentioned, the principal quantities to consider in dilution prediction
are volume rate of discharge (Q) or volume rate of discharge per unit diffuser
length (q), ambient stratification (dp/dz), initial density difference between
waste and ambient (A0), and ambient currents (U). These quantities should be
considered at:
a) times of seasonal maximum waste outflow,
b) times of seasonal maximum and minimum stratification,
c) times of low ambient water quality,
d) times of low net circulation or flushing, and
e) times of exceptional biological activity.
These periods have durations generally of 1 to 3 months. The quantities se-
lected to represent these periods should reflect conditions of the worst 10
percent of occurrences during these periods; highest volume rate of dis-
charge, extreme ambient stratification, lowest initial density difference, and
9
-------
lowest current speed. If a reliable estimate can't be made from the available
data, a conservative estimate should be substituted, e.g. zero currents.
To select conditions at the worst 10 percentile, they must be ranked in
some way. Current speed data usually consist of discrete values and can be
ranked into cumulative frequency distributions by computer or by hand. It is
more difficult to rank existing data on ambient stratification and discharge
rate. A suitable value to use for the volume rate of discharge is that value
exceeded 2 to 3 hours during the average day for the period. The average over
the period of daily values exceeded 2 to 3 hours per day or the value exceeded
2 to 3 hours on a day which has average total flow over the period could be
used, depending on the form of the flow data. Data from the last 2 years or
longer should be used to insure that the flows are representative.
The worst stratification is that which is greatest over the height of
rise of the plume. This condition will result in the lowest dilution if other
conditions are constant. If the density gradients are uniform with depth they
can be ranked numerically. If not, measured profiles can be input to the
computer models PLUME or OUTPLM and the results used to rank the profiles.
Alternatively, the average or overall density gradient over the height of rise
of the plume can be estimated by the method given in the next Section (Non-
Uniform Stratification) and these average gradients ranked. The density dif-
ference between the ambient and the effluent will probably vary only a small
amount and will probably be due to changes in the ambient density. The dens-
ity difference can, therefore, be selected by ranking the ambient density at
the level of the discharge.
In unusual cases, the discharge ports may resemble thin (< 1/10 D),
sharp-edged orifices. In these cases, the port diameter may be taken as 0.62
times the actual diameter of the orifice and the discharge velocity adjusted
10
-------
accordingly. This will increase the computed Froude number of the discharge.
Where the orifice is rounded or thick, as is generally the case, no correction
should be made.
MIXING ZONE SPECIFICATION
After initial dilution, the concentration of waste constituents (Cf) is a
function of the average dilution achieved (Sa) and their concentration in the
ambient (Ca) and the effluent (Ce):
Cf = Ca + (Ce-Ca)/Sa
If the effluent has been adequately treated and disposed of in an environment-
ally appropriate area, Cf's for its various constituents should be at or below
water quality standards.
The zone surrounding the discharge site geometrically bounding the criti-
cal initial dilutions is termed the zone of initial dilution (ZID) to disting-
uish it from other mixing zone definitions. It is a concentration isopleth.
The ZID describes an area in which inhabitants, including the benthos, may be
chronically exposed to concentrations of pollutants in excess of water quality
standards or at least to concentrations greater than those predicted for the
critical conditions described above. The ZID does not attempt to describe the
area bounding the entire initial mixing process for all conditions or the area
impacted by the sedimentation of organic material.
The purpose for defining a zone around the discharge point is to disting-
uish biological impacts caused by exposure of organisms to concentrations
higher than Cf from those impacts caused by the re-concentration of pollutants
by physical, chemical, and biological means. The former is unavoidable with-
out pre-dilution or treatment of the waste to meet water quality standards.
The latter indicates that dilution alone is not sufficient to ameliorate the
11
-------
environmental and ecological impacts of the discharge and additional treat-
ment, relocation of the discharge, or influent source control may be neces-
sary.
The ZID can be specified by analyzing prediction results for a range of
conditions. However, it can be more simply specified using the maximum height
of rise predicted for the critical conditions as a radial distance measured
horizontally from the outfall diffuser or port. This distance will often
equal the depth of water at the discharge site. During periods of higher
currents, the plume trajectory will be more horizontal and initial dilutions
higher than predicted for the critical conditions. The dilution achieved over
that portion of the trajectory within the ZID, however, will be approximately
equal to the initial dilution predicted for the critical conditions.
12
-------
SECTION 3
SUMMARY AND CONCLUSIONS
GENERAL
The plume attributes of interest here are the dilution attained at the
end of convective ascent through ambient waters, and the plume's position.
For our purposes, the poorest or lowest dilutions are important. Dilution is
defined as the total volume of a parcel or sample divided by the volume of
waste contained within it. This is equivalent to the inverse of the volume
fraction (or concentration) of waste that a parcel or sample contains.
The behavior of buoyant plumes can be mathematically modeled by properly
considering mass, momentum, and energy. Some form of entrainment function
must be assumed and fitted to experimental data. Other assumptions generally
made are that all flows are steady and incompressible, pressure is hydrostatic
throughout, the plume iis fully turbulent and axisymmetric, and that turbulent
diffusion dominates and is significant only in the radial direction. Distri-
butions of velocity and concentration may also be assumed. The physics of
this problem can be solved in various ways. Systems of differential equations
can be integrated across the plume, reducing the equations to one independent
variable, length along the axis of the plume. Difference equations can be
used with Lagrangian coordinants. Simultaneous equations or three-dimensional
numerical schemes also have been used.
Mathematical models of buoyant plumes are systems with internal variables
(mass, momentum, and energy), external variables (discharge characteristics,
ambient vertical density, and currents) and internal mechanisms (entrainment).
13
-------
These models have behavior which is affected by principal quantities. Compu-
ter model studies can be used to describe plume behavior of interest in terms
of effecting principal quantities. For stagnant conditions, the principal
quantities are initial density difference between the waste and ambient (A0),
ambient stratification (dp/dz), and flow rate or flow rate per unit length for
multi-port diffusers (Q or q). Current speed (U), if considered, is also a
principal quantity. Froude number is not important to the dilution attained
in the neighborhood considered (0 < Fn < 30), especially in currents. The
Froude number and orientation of the discharge affect the plume's trajectory
in stagnant conditions but have only a small effect in a flowing ambient. The
spacing of multiple ports, for a given discharge rate per diffuser length,
does not significantly effect the dilution attained by merging plumes in a
flowing ambient. In stagnant conditions, port spacing has a moderate (perhaps
20 percent) effect on initial dilution.
Closed-form solutions for the behaviors or attributes of interest (dilu-
tion at maximum height of rise and height of maximum rise) have been found in
semi-theoretical, state-transition modeling. These solutions have been mapped
or tuned to the behavior of the computer models. Plume behavior can, there-
fore, be solved analytically in terms of the principal quantities. These
solutions also allow the importance of the principal quantities to be seen
directly. The computer models are still useful for detailed analysis or spec-
ial conditions. In many cases, however, the analytic models are quite ade-
quate to calculate dilution and maximum height of rise.
The discharge configuations which will be considered include single
ports, and multiple ports spaced such that they interact at some point during
convective ascent. These configuations are termed, respectively, plumes, and
merging plumes. Some other restrictions and assumptions must be made. The
14
-------
discharge site is assumed to be at least ten port diameters deep and uncon-
fined. Flow in the ambient is assumed to be uniform. Each of the three
models used in this study is suited to certain combinations of conditions.
Recommendations are given in Table 1. Application criteria for the analytic
or simplified models are summarized on Figures la and lb at the end of this
section.
TABLE 1. RECOMMENDED MODEL APPLICATION.
Discharge Ambient Ambient Surfacing Analytic Computer Models
Configuration Density Currents Plumes? Models PLUME OUTPLM DKHPLM
Single Plumes Stratified Stagnant Yes X X
" "No XX
11 Flowing Yes X X
No X X
Uniform Stagnant Yes X X
11 Flowing Yes X X
Merging Plumes Stratified Stagnant Yes X X
No X X
" Flowing Yes X X
11 No X X
SIMPLIFIED RESULTS FOR PLUMES IN A STRATIFIED AMBIENT
Analytic solutions are given in this section for dilution and height of
rise of plumes and merging plumes. For stagnant conditions, the height of
rise is taken to the trapping level. The trapping level is where the plume
first loses its positive buoyancy. Above this level the plume is expected to
"pool" and spread horizontally. For flowing conditions, height of rise is
taken to the centerline of the plume where it first becomes horizontal.
Dilutions predictions are made for these locations.
The results described cover ambient conditions including vertical density
stratification, both stagnant and flowing, and the previously described dis-
15
-------
charge configurations. Vertical density stratification is normally found in
the marine environment and, since it adversely affects dilution, should be
considered. A summary of the analytic solutions covering the conditions and
configurations listed in Table 1 will be given. Following these, sections on
non-uniform stratification and single plumes in uniform-density ambients are
presented. SI units (kilograms, meters, seconds) are used.
Plumes in a Stagnant Environment
The dilution for a single port discharge into a stably-stratified and
quiescent ambient is given by::
Sa = 0.46 Q-1/4 A03/4 (dp/dz)"5/8 (1)
where Sa is the flux averaged dilution at the trapping level. The height of
the trapping level can be determined by:
, Ah = 4.6 (A0 Q)1/4 (dp/dz)"3/8 (2)
The dilution of the plume is most affected by initial density difference
(A0) and ambient density gradient (dp/dz), and secondarily by the volume rate
of discharge (Q). Flux of buoyancy and ambient density gradient both affect
height of rise. The height of rise (Ah) given by equation (2) is a maximum
value for the trapping level. This expression does not consider the angle of
discharge from the vertical or the Froude number which can reduce the actual
trapping level height. If, in application, it is found that the height of
rise given by equation (2) exceeds the water depth (Z) at the discharge site,
dilution can be calculated by:
Sa = 0.10 Z (A0/Q)1/2 (dp/dz)"1/4 (3)
This relationship assumes linear dilution with height of rise which is not
strictly observed in stagnant conditions. As the water depth becomes smaller
than Ah, equation (3) predicts dilution which becomes less than those pre-
16.
-------
dieted by the PLUME program. If greater accuracy is needed the PLUME program
should be used.
Equation (1) is in good agreement with model results for jets with Froude
numbers of up to about 50. For jets with Froude numbers above this value the
PLUME program should be used.
Example—
A discharge has a maximum flow rate of 0.10 m3/sec and a density of
1000 kg/m3. The ambient at the discharge site has a density at the surface of
1020 kg/m3 and, at 30 m depth, a density of 1025 kg/m3 during periods of maxi-
mum stratification. From equations (1) and (2):
Sa = 0.46 (0.10)"1/4 (25)3/4 (0.166)~5/8
= 28
Ah = 4.6 ((0.10)(25))1/4 (0.166)~3/8
= 11.3 m
If the depth over the discharge were only 10 m and other conditions the same,
then from equation (3):
Sa = 0.10 (10.) (25/0.10)1/2 (25)1/2 (0.166)"1/4
- 25
Plumes in a Flowing Environment
The dilution and height of rise of a single plume discharged into a
stably-stratified, flowing ambient are:
Sa = 3.0 (U/Q)1/3 (A0/(dp/dz))2/3 (4)
Ah = 2.3 (A0 Q/U(dp/dz))1/3 (5)
17
-------
The dilution of the plume is most strongly dependent on the ratio of initial
density difference to ambient density gradient. Equation (4) will predict
increased dilution (greater than equation (1)) when:
U > 0.0036 (AoQ)1/4 (dp/dz)1/8
Neither the dilution nor the height of rise are sensitive to discharge angle
or typical Froude numbers. The maximum height of rise is for the centerline
so that the plume actually extends one radius up from this level. The plume
1/2
radius can be estimated by: B = (SaQ/nU) If Ah + B exceeds the depth of
1/3
water at the discharge site (when Z < 3.3 (A0Q/U(dp/dz)) ), the computed
dilution can be corrected by:
Sa (corrected) = Sa • Z/(Ah + B) (6)
where B is the plume's radius or half width. By substitution into equation
(6), the dilution of a plume which reaches the surface is:
Sa = 0.92 Z (U/Q)2/3 (Ac/(dp/dz))1/3 (7)
Equation (7) is applicable when U > 0.036 (A0Q)^ (dp/dz)^®, otherwise
use equation (3).
Example—
In the previous example, A0 = 25 kg/m3, Q = 0.10 m3/sec, dp/dz =
0.166 kg/m3/m. If the lowest 10 percentile of current speed observations is
0.10 m/sec, what is the dilution and height of rise?
Sa = 3(0.10/0.10)1/3 (25/0.166)2/3
= 85
Ah = 2.3((25)(0.10)/(0.10)(0.166))1/3
= 12.2 m.
1 /?
The radius can be estimated as: B = (56(0.10)/n(0.10))
= 5.2 m
18
-------
If the water at the discharge site is only 10 m deep, dilution as cor-
rected by equation (6) is:
Sa (corrected) = 85(10)/(12.4 + 5.2)
= 49
Alternately, dilution can be re-computed using equation (7):
Sa = 0.92 (10) (0.10/0.10)2/3 (25/0.166)1/3
= 49
Merging Plumes in a Stagnant Environment
Brooks (22) gives expressions for the dilution and height of rise of slot
plumes in stagnant, stratified ambients. Slightly re-arranged, re-dimensioned,
and tuned to DKHPLM for consistency with the other analyses presented here,
dilution can be predicted by:
Sa = 0.32 q~1/3A02/3(dp/dzf1/2 (8)
The height of rise of the merging plume is:
Ah = 5.9 (A0q)1/3 (dp/dz)"1/2 (9)
The height of rise is predicted from equation (9) is the maximum trapping level
height.
Widely spaced plumes which merge after a significant rise will attain a
higher dilution (perhaps 20 percent) than predicted by equation (8).
If the height of rise of the plumes exceeds the depth of water, dilution
can be calculated by:
Sa = 0.054 Z q"2/3 A01/3 (10)
This dilution will be lower than predicted by simulation model DKHPLM because
of the assumed linearity of dilution with height of rise and because the trap-
ping level is assumed to be the maximum. If more accuracy is needed, DKHPLM
should be used.
19
-------
Example—
A diffuser discharges 0.02 m3/sec/m in 50 m of water. If the efflu-
ent has a density of 1000 kg/m3, A0 = 25 kg/m3, and the ambient density gradi-
ent is 0.05 kg/m3/m, the dilution by equation (8) will be:
Sa = 0.32(0.02)"1/3 (25)2/3 (0.05)~1/2
= 45
The maximum height of rise to the trapping level, from equation (9), will be:
Ah = 5.9 ((0.02)(25))1/3 (0.05)"1/2
= 20.9 m
If the water depth were 10 m, the dilution, from equation (10), would be:
Sa = 0.054 (10) (0.02)~2/3 (25)1/3
= 21
Merging Plumes in Flowing Environment
For multi-port diffusers where plumes interact or for slot plumes, the
principal quantities which affect dilution and height of rise of the plume are
initial density difference (A0), discharge rate per unit length of diffuser
(q), current speed (U) and ambient density gradient (dp/dz). The dilution from
these discharges can be found by:
Sa = 2.7 (A0U/q(dp/dz))1/2 (11)
and the height of rise by:
Ah = 1.3 (A0q/U(dp/dz))1/2 (12)
Like single plumes, the dilution and height of rise of multiple plumes in a
current are not sensitive to Froude number and discharge orientation. Appar-
ently the port spacing or level of interaction also have little effect on the
values ultimately attained by the plume.
Currents will increase the dilution attained by merging plumes when:
20
-------
U > 0.014 (A0q)1/3
If currents are not perpendicular to the diffuser, the discharge rate (q) can
be adjusted to consider the length of the diffuser projected by the current.
Below a relative angle of about 45 degrees, however, prediction should be
based on stagnant conditions.
The height of rise of merging plumes is about equal to the half-width.
1 /2
If Z < 2.6 (A0q/U(dp/dz)) the plume will surface and its dilution, cor-
rected by equation (6), calculated by:
Sa = 1.0 Z U/q (13)
The dilution of surfacing, merging plumes is affected by currents when:
U > 0.054 (A0q)1/3
These analyses, as well as DKHPLM, assume that the diffuser length is very
long compared to the height of rise. Otherwise edge effects complicate the
problem.
Example—
A discharge of 12.5 m3/sec 300 MGD) is made in 50 m of water.
The ambient density during maximum stratification varies from 1024 kg/m3 at
the surface to 1025 kg/m3 at the bottom. The diffuser is 1250 m long so that
the discharge rate per unit length of diffuser (q) is 0.01 m3/sec/m. Density
of the effluent is 1000 kg/m3. If the lowest 10 percentile of current speed
occurrences is 0.05 m/sec, the dilution will be about:
Sa = 2.7((0.05)(25)/(0.01)(0.02))1/2
= 213
And the expected height of rise will be:
Ah = 1.3((25)(0.01)/(0.05)(0.02))1/2
= 21 m
21
-------
If the discharge site were only 20 m deep, equation (13) predicts the dilution
to be:
Sa = 1.0 (20)(0.05)/0.01
= 100
Non-Uniform Stratification
The preceding analyses are for uniform, or monotonical, density gradi-
ents. Often density gradients deduced from field measurements are not
strictly monotonical but can be approximated as such. For cases where the
vertical density distribution is interrupted by a sharp discontinuity or
pycnocline, other approximations are necessary. Brooks (22) gives one method.
The height of rise can first be estimated by using the density gradient of the
lower layer. If this height of rise extends to the upper layer, the density
gradient is taken as the average gradient to that point. It may be necessary
to re-estimate maximum height of rise if the first adjustment resulted in a
large change in the average density gradient. The density gradient can then
be re-estimated. By this method, an average density gradient over the ex-
pected height of rise of the plume can be estimated and dilution calculated.
This procedure is useful for evaluating stratification data and identifying
the most adverse stratification. If the density structure of the ambient is
very complex, or if more accuracy is needed, then PLUME or OUTPLM models
should be used.
PLUMES IN UNIFORM DENSITY AMBIENTS
Conditions may exist at a discharge site which mix the ambient.to a uni-
form state. Normally these are not the conditions of interest since ambient
stratification is more critical to the initial dilution process. The analytic
-------
solutions for single plumes in stratified ambients, because of their form,
approach infinity as the ambient density gradient approaches zero. Thus, they
may over-estimate dilution when stratification is slight. Different analytic
solutions are therefore necessary for single plumes in uniform ambients. The
analytic solutions for merging plumes in stagnant and flowing conditions,
equations (10) and (13), are applicable to uniform density ambients.
Analytic solutions will be given here for cases of single plumes in stag-
nant and flowing ambients with little or no stratification. For a single
plume in a stagnant, uniform ambient:
Sa = 0.074 Z5/3 Q~2/3 (14)
1 /4
For a single plume in a uniform ambient flowing at U > 0.1 (Q/Z) , dilution
can be calculated by:
Sa = 1,6 Z2 U4/3 Q_1 (15)
Equation (15) should be applied when:
dp/dz < 0.18 Z"3 if2 Q A0 (16)
The criteria presented for the application of equation (14) in Figure la is of
limited value because of the uncertainty in predicting dilutions in stagnant,
stratified ambients when the height of rise is limited and much less than the
trapping level height. In very siight stratification, this is the case.
However, in either the flowing or the stagnant case, results from equation
(14) or (15) should be compared with equation (3) or (6) when stratification
is slight. The lesser result is correct since stratification should reduce
the dilution of single plumes.
23
-------
PLUMES
Stagnant
F1owi ng
U < 0.0036 (A0Q)1/4 (dp/dz)1/8 < U
Stratified Sa =
0.46 Q"1/4 A03/4 (dp/dz)"5/8
Sa =
3.0 (U/Q)1/3(A0/dp/dz)2/3
Z < 4.6 (A0Q)1/4 (dp/dz)"3/8 Z < 3.3 (A0Q/U(dp/dz))1/3
U < 0.036 (A0Q)1/4 (dp/dz)1/8 < U
Surfacing
Sa =
Sa =
0.10 Z (A0/$)1/2 (dp/dz)"1/4
0.92 Z (U/Q)2/3 (A0/dp/dz)1/3
dp/dz < 3.
3 z"8/3Z2/3 Ao2
dp/dz < 0.18 Z"3 I
l"2 A0Q
U < 0.1 (Q/Z)1/4 < U
1
Uniform
Sa = 0.074
z5/3 q-2/3
Sa = 1.6 Z2 Q"1 U4/3
Figure la. Application criteria summary for the simplified analysis of plumes
in various ambients.
24
-------
MERGING PLUMES
Stagnant
Flowing
1
U < 0.014 (A0q) < U
Stratified
Surfacing/
Uniform
1/3 (dp/dz)"1/2 Z < 2.6 (A0q/U(dp/dz))1/2
1/3
U < 0.054 (A0q)
Sa =
2.7 (A0U/q(dp/dz))1/2
Sa = 1.0 Z (U/q)
-2/3 A 1/3
Sa =
-1/3 A 2/3
(dp/dz)
-1/2
Figure lb. Application criteria summary for the simplified analysis of merging
plumes in various ambients.
25
-------
SECTION 4
"PLUME" MODEL
THEORETICAL DEVELOPMENT
The computer model PLUME (5, 19) considers a buoyant plume issuing at an
arbitrary angle into a stagnant, stratified environment. Two zones of plume
behavior are considered. Profiles of velocity and concentration are considered
to be uniform or "top hat" at the port orifice. The region required to develop
fully established profiles is the zone of flow establishment. Beyond this
point, similarity profiles are assumed for velocity and concentration. This
allows the equations for mass continuity, momentum, and density disparity to be
integrated across the plume, reducing them to one independent variable, the
length along the axis of the plume.
Zone of Flow Establishment
The length of flow establishment (sg) has been found to depend on the
initial Froude number (12). As initial Froude numbers approach infinity,
s = 5.6 D. At low Froude numbers, this distance can be much shorter. The
e
centerline value of concentration is assumed equal to the initial discharge
conditions. The inclination of the plume's axis at sg is determined by linear
interpolation using 0, sg, and Fn. The centerline velocity at sg depends on
the Froude number and, at low Froude numbers, may be greater than the initial
velocity. To avoid estimation, the velocity at s (V ) is determined by s^
e me e
and empirical spreading coefficients.
-------
Zone of Established Flow
The following assumptions are made:
a) flow is steady and incompressible,
b) the plume is fully turbulent,
c) turbulent diffusion is significant only in the radial direction,
d) pressure is hydrostatic throughout,
e) plume flow is axisymmetric.
It is also assumed that the distributions of velocity and concentration are
Gaussian across the plume, so that:
V/Vm = exp (-k (r/s)2) (17)
C/C^ = exp (-k|j(r/s)2) (18)
where V and C are the velocity and concentration, respectively, at some radial
distance, r, and the subscript m refers to the centerline. The distance along
the plume's trajectory is s, and k and n are the empirical coefficients.
The integrated governing equations are:
conservation of mass,
d/ds /°° V r dr = -r. u. (19)
o b b
conservation of momentum,
d/ds f° V2 r dr = sine /° (p0o_p)g r/pd dr (20)
conservation of density disparity,
d/ds f V A r dr = /" V r dr • dp^/ds (21)
O O
(Po-Pd)
conservation of pollutants,
d/ds f V C r dr = 0 (22)
O
Substituting in the velocity profiles (17) and (18), and simplifying
yields:
27
-------
momentum,
d/ds (Vm3 s3/k3/2) = 3 g((p0-pd) Vm Am s pk3/2 sine (23)
density disparity,
d/ds (Vm Am s2/k(p + 1)) = (l/(Po-Pd)) Vm s2/k-dPo6/ds (24)
conservation of pollutant,
Vm Cm s2/k(|J + ]) = °2 Vo Co/4 (25)
The equations can be cast dimensionlessly as:
momentum,
dE*/ds* = 3/Fn2-((1 + M)/pk1/2)s* R* sine (26)
density disparity,
dR*/ds* = sVk1/2 • E*1/3 (dp£/ds*) (27)
conservation of pollutant,
Cm = (k1/2 (M + 1 ))/(4E"tl/3 s*) (28)
where s* = s/D, V* = Vm/V0, E* = (VJ s*/k1/2)3, R* = V* Am s*2/k(|j + 1) =
e*1/3 sVkl/2 (m + 1)( and = Poo/(Po-pd)
The initial conditions for the zone of established flow are obtained by evalu-
ating E* and R* at s* = sg*. It is assumed that ambient stratification is
negligible in the zone of flow established so A„ = C and R * = 1/4. Since
33 m m e
1/2 3
Cme = C0, Eg* = (k Cm + l)/4 se*) The values of p and k are mean values
of the limiting cases of purely-buoyant and purely-momentuni jets found by
Abraham (12).
The angle of the plume (6) can be evaluated by considering conservation
of horizontal momentum,
d/ds f° V2 cos0 r dr = 0 (29)
O
which combined with equation (17) yields
6 = cos'1 ((Ee*/E*)2/3 cos0e) (30)
28
-------
MODEL DESCRIPTION
The zone of flow established is first solved depending on the Froude
number:
s* = 2.8 Fn2/3 Fn < 2
e
s* = 0.113 Fn2 +4. 2 < Fn < 3.1
e
s* = (5.6 Fn2)/(Fn4 + 18)1/2 Fn > 3.1
Then other initial conditions, 9e, Eg*, and Re*, are solved for. Equations
(26) and (27) are then solved at steps along the plume's axis by Runge-Kutta
approximation. E* is used to calculate the inclination 6 and concentration C.
Then the position of the plume's axis can be incremented by dx/ds = cos0 and
dz/ds = sin0, where x and z are horizontal and vertical coordinates, respec-
tively.
The centerline concentration is inverted to describe the centerline dilu-
tion. To convert the centerline dilution to a flux averaged dilution, the
distribution of concentration must be weighted by distribution of velocity.
With the assumed distributions, the flux average dilution can be found by:
Sa = 1.77 (1/CJ or = 1.77 Sm (31)
m m
where S is the centerline dilution from PLUME. Calculation of dilution is
m
terminated when the density of the plume equals the density of the ambient.
This is normally somewhat below the maximum height of rise but is where simi-
larity ends and above which the plume will spread and become passive, possibly
interfering with further dilution.
Computation of time and plume diameter have been added to the original
program. Those changes are set off by asterisks in the program listing given
here. The centerline velocity, averaged over each step length, is divided
into the step length to obtain time. The plume diameter is found by D = 0.308
s.
29
-------
EXAMPLE INPUT, OUTPUT, AND MODEL LISTING
Input
The example problem is that of a single port with a diameter of 0.25 m
discharging 0.1 m3/sec of effluent with a density of 1000 kg/m3. The dis-
charge is horizontal. The discharge site varies in density from 1020 kg/m3 at
the surface to 1025 kg/m3 at the discharge depth of 30 m.
A schematic of the computer-punch input cards for this problem follows.
30
-------
I 1
t
1
5
Header Information
)
8
1
Plume Test. Case
r
s
tr
1
L-.,
~ 1 1
A
; T i i ii i ; 1 r i ( i i i i
-n !—rr-
TTTT"f
r riTTTi—t— • 11 i —
i
ii "¦ f i i i i 1 i : ¦ ! i i
' I
tf
I 1
^rii?rL
ir
_L_!
JSj
£• 4
I
P o P G £
I 8 ! f g
q CJ pn 3
•n
f S® g s
« 5 c
3* I ?
« S.2
© 3,8
« 0
tfl 3
5- 5 3
(9
X
CU
3 O
¦O O) T3
-«-s r
Q Q.C
3
T3
-------
CO
ro
DEPTH
1 13 F5.0
F10.0
F10.0
Fin.n
i 1°'
METERS
ANGLE
DEPTH
Units, 1=MKS, blank^FPS
Number of ambient density points to be specified
Port angle measured from the horizontal
Port diameter, m.
Depth of port
Density of effluent, gm/cm
Blank
Total volume rate of discharge
Number of ports
Print out interval
filMIT WTTIT1
hh ro m
-------
H..HL —
DP(r)
o.o
30.0
RHO(I) |
1.0200 |
1.0250
F10.0
F10.0
Im
lC>!
SAL
Card 3
F10.0
Mi
Card 4
rm i
icj
in
nnnmi
1 !
DP (I) - Depth measured from the surface
RH0(I)- Density (gm/cm3), or temperature (°C) if Sal specified
SAL - If density not specified in RH0(I) then salinity, o/oo.
-------
Output
The first part of the model output repeats the input and provides other
computed initial conditions. The Froude number given is the square of the
definition used here. The discharge velocity and characteristics of the zone
of flow establishment are also given in this section.
The next section of the output traces the plume and its characteristics
along its trajectory. T is a time based on the centerline velocity. S is the
distance along the centerline. X and Z are distances measured horizontally in
the direction of the discharge and vertically from the water's surface to the
center of the plume, respectively. Elevation is the vertical distance to the
centerline of the plume from the level of the discharge. D is the diameter of
the plume and DILN is the centerline dilution. The dilution and height of
rise to the trapping level are also given. As mentioned, equation (31) is
used to convert this to a flux averaged dilution. In this case, Sa = 1.77 Sm
= 27.3. The example output and program listing follow.
34
-------
S3 PLUME VERSION 2.3 9/12/77
GJ
cn
BUOYANT plume in a density stratified mediae
PLUME TEST CASE
CASE NO. 1 INITIAL CONDITIONS ....
UNITS: mks
PORT ANGLF ...... . . 0
FROUOE NUMRER 67.7
LENGTH FOR FLOW ESTABLISHMENf . . . 1.40
INTEGRATION STEP LENGTH .076
PRINTOUT INTERVAL . 3.00
XO 1.39
ZO .......... 29.84
DISCHARGE DENSITY . ... 1.00000
PORT DEPTH 30.00
FLOWRATE 1.OOOOOE-O1
NUMBER OF PORTS 1
DISCHARGE VELOCITY 2.04
PORT DIAMETER 2.50000E-01
DENSITY STRATIFICATION DEPTH RHO
0 1.02000
30.00 1.02500
T
S
X
Z
D
ELEV
THETA
2.35
4.29
4.06
28.85
1.32
1.15
36.2
6.43
7.25
5.99
26.62
2.23
3.38
59.1
11.51
10.22
7.27
23.95
3. 15
6.05
60.1
17.64
13.18
8.30
21.17
4.06
8.63
70.3
25.61
16. 14
9.41
18.42
4.97
11.58
63.1
33.91
18.12
10.74
17.07
5.58
12.93
0
LAST LINE ABOVE IS FOR MAXIMUM HEIGHT OF RISE.
TRAPPING LEVEL IS 21.5 WITH DILUTION OF 15.4321
HEIGHT OF RISE= 28.3 PERCENT OF DEPTH
OILN
3.3354
7.0703
11.6926
-------
0001 PROGRAM PLUME
0002 COMMON G»FK.FM,COSTH.SINTH.COSTHE.DS.Cl.C2»E13.FLAG.GRAV
0003 DIMENSION ZD(SO).00(50).RHO(50).DP(50).CMT(5)
0004 INTEGER NDC,TRAPPD,FLAG.CHGOEN.METERS
0005 WRITEI6.110)
0006 110 FORMAT("0S3 PLUME VERSION 2.3 9/12/77")
0007 NOCASE=0
0008 CUBRT = 1,/3.
0009 TW0THD=2./3.
0010 20 READC5.222)(CMT(I),1=1.5)
0011 222 FORMAT(5A8)
0012 IF (EOF(5)) CALL EXIT
0013 READ(5.73) NDC.METERS,NPTS.ANGLEtDIA.DEPTH,RHOJ.RFD.Q.FN.PS
0014 IF(EOF(5)) GOTO 7*
0015 73 FORMAT(2I1»I3»F5.0»7F10.0)
0016 FD=A9S(RFD)
0017 DO 102 1=1.NPTS
0018 RE AD(5.75) OP(I).RHO(I).SAL
0019 IF(EOF(5)) GOTO 999
0020 IFISAL.EQ.O.) GOTO 102
0021 RHO(I)=1.*.0 01»SIGMAT(SAL.RHO(I)>
¦0022 102 CONTINUE
¦0023 CALL ASORT(OP,RHO,.NPTS)
10024 999 IF(NPTS.NE.1) GOTO 76
10025 NPTS=2
10026 DP(2)=DEPTH
0027 RH0(2>=RH0<1)
0028 76 N0CASE=N0CASE*1
0029 75 FORMAT(8F10.0)
0030 6RAV=32.172
0031 IF(METERS) GRAV=9.80665
0032 DO 55 1 = 1.NPTS
0033 IF(DP(I).GE.DEPTH) GOTO 56
0034 55 CONTINUE
0035 WRITE(6.59) NOCASE
0036 59 FORMATC'-NO DENSITY INFORMATION FOR JET LEVEL. EXECUTION FOR".
0037 « » CASE NO.".12." DELETED.")
0038 T0 20
0039 56 NP=I
•0040 NM=I-1
0041 RHOB =(OEPTH-DP(NM))»(RHO(NP)-RhO(NM))/(DP(NP)-DP(NM))*RHO(NM)
0042 niSP=RHOJ-RHOH
0043 00 54 I=1.NM
•0044 J=NP-I
0045 ZD(I)=(OEPTH-OP(J))/DIA
• 0046 54 r>G(I) = (RHO(J.l)-RHO(J> ) »D I A/(DI SP» (OP (J* 1) -DP (J)) )
:0047 IF(NDC) GOTO 58
'0048 U0=Q/(FN«.7853962»0IA*DIA)
>0049 RFD=UO»UO»RHOJ/<-DISP»DIA»GHAV)
¦0050 FO=A0S(RFD)
10051 58 .IF (FD.LE.4.01.0R.FD.GE.9.99) GOTO 61
'0052 S=.113»FU*4.
>0053 GOTO 62 j *-o
10054 61 IF (FD.LE.4.01) S = 2.8«FD»«CUBRT /
10055 lF if/fO.Q^ <1.3 ^ S = 5". 0/SV(?.f IF0¥
10056 62 TEMP=ATANU.416667«S/FD)
) 00 5 7 THETA0=.Oil1«ANGLE«(1.5708-TEMP).TEMP
10058 C0STHE=C0S(THETA0)
10059 SINTHE=SIN(THETAO)
10060 DSI=DEPTH/(177.»DIA)
10061 IF(OGU) .EO.O.) GOTO 77
10062 DGTEMP=.01/DG(1)
J0063 IF(DGTEMP.LE.O.) DGTEMP=-DGTEMP
36
-------
00064' "hS I =. 12» 1.*» f#2^® ( ACtT5lrb'*l DGTEMPj )•
00065 77 NPO=IFIX(PS/(DSI»DIA)~.5)
00066 N=0
00067 Z=S»SINTHE
00068 X=S»COSTHE
00069 F=<4./S)«<»3
00070 R=.25
00071 C1=E**TW0THD
00072 C2=•75/RFD
00073 IPTS=1
00074 6=PG1IPTS)
00075 7LIM = ZDUPTS>
00076 CHGDEN=0
00077 FLAG=0
00078 TRAPPD=0
00079 XP-X*DI A
00080 7P=DEPTH-Z«DIA
00031 nSIP=DSI«DIA
00062 SP=S»DIA
00083 C
00084 T=2.«SP/(U0«tl.»336.67«(S«S)))
00085 C o»oo«oooo(H>oooo»o»ooo«oooooo»ooo4ooeooooooo»ooeoi>
00086 C
00037 C WHITE INITIAL CONDITIONS
00088 C
00089 WRITE(6.23) (CMT(I>* I = 1,5).NOCASE
00090 23 FORMAT C •• 0»A BUOYANT PLUME IN A DENSITY STRATIFIED "
00091 «»MEDi A»«"«<»°»/"0"»5A8/
00092 »"-CASE NO. »,I2»" INITIAL CONDITIONS
00093 IF (NDC) WHITEI6.104)
00094 104 FORMAT.("0 UNITS: NON DIMENSIONAL")
00095 IF(.NOT.NDC.AND.METERS) WRITEI6.105)
00096 105 FORMATC'O UNITS: MKS")
00097 TF(.NOT.NDC.AND..NOT.METERS) WRITE(6»106)
00098 106 FORMATC'O UNITS:FPS")
00099 WRITE(6.108)ANGLE.RFD.SP.DSIP«PS.XP,ZP.RHOJ,DEPTH
00100 10d FORMAT(
00101 »40H- PORT ANGLE
00102 «40H0 FHOUDE NUMBER F7.1/
00103 *4OH LENGTH FOR FLOW ESTABLISHMENT . . ..F0.2/
90104 «40H INTEGRATION STEP LENGTH F9.3/
00105 »40H PRINTOUT INTERVAL .F8.2/
00106 «40H XO »F8.2/
00107 °40H 20
00108 «40H DISCHARGE DENSITY .F11.5/
00109 «40H PORT DEPTH F8.2)
00110 IF(.NOT.NDC) WR1TE(6,60 IQ.FN.UO,01A
00111 60 FORMAT(
00112 «40H FLOWRATE 15.5/
00113 »40H NUMBER OF PORTS »F5.0/
00114 «40H DISCHARGE VELOCITY «F8.2/
30115 °40H PORT DIAMETER »E15.5)
00116 WRITE (6»24) ((DP(I)tRHO(I)>.I = 1.NPTS)
00117 24 FORMAT(DENSITY STRATIFICATION DEPTH RHO"//
00118 «10<24X.F6.2,F11.5/>)
00119 C
00120 C SET INITIAL CONDITIONS
00121 C
00122 DS=DSI
00123 WRITE(6.109)
00124 C
00125 109 FORMAT < lH-,6X»"T".9X,»S".9X.,,X'»,9X«»Z"t9X.»D"»7X»"ELEV'S
00126 •6X."THETA"f4X."DILN")
00127 C
00128 GOTO 16
00129 11 DS=DSI
37
-------
urso
45
DElX^COS fH»DS
0131
PELZ=SINTh«OS
'0132
DELE=FK
'0133
OELT=FM
0134
SP0S2=S*DS*.5
0135
no 10 1=1,2
0136
CALL SDEWIV(SPD52,E».5«FK,R«.5*FM)
0137
0ELX = DELX*2.»C0,iTH«DS
>0138
0ELZsUELZ*2.#SINTH«DS
'01 39
0ELE=DELE*2.#FK
10140
0ELT=DELT*2.»FM
>0141
10
CONTINUE
'0142
CALL SOERIV(S*OS»E*Fk»R*FM)
'0143
ZLAST=Z
>0144
7INCR=(DELZ*SINTH«DS)/6.
>0145
Z=Z*ZINCW
>0146
IF(CHGOEN) GOTO 41
0147
IF(Z.GT.ZLIM) GOTO 40
<0148
4 3
X = X*(DELX*COSTH»OS)/6.
;0149
C
« « o o » « « o « o o e » « o « « o « « o « » « « « ft o tt
>0150
E0LD=E
>0151
C
ft»»e»»»«»»«»»«»»ft»«ft»»»«fto««ft
>01 52
E=E*(DELE*FK)/6.
J0153
R=R*(0ELT»FM)/6.
10154
C
10155
0T=.218«DIA»0S/(U0«(E0LD,l«(l./3.)/S*E®«0159
IFIE.LE.O.) FL AG=1
>0160
c
10161
c
THIS STOPPING CRITERIA IS BASED ON VELOCITY GOING TO ZERO
>0162
c
>0163
IF(FLAG) GOTO 13
>0164
16
CALL SDERIV(S,E»R)
>0165
IF(TRAPPD) GOTO 14
>0166
IFIR.GT.O.) GOTO 15
10167
RAT =R/(RO-R)
>0168
F13TRP=E13*(E13-E130)
>0169
7TRAP=DEPTH-(Z*ZINCR»RAT)»OIA
>0170
SMTRAP=.245»(S*DS«RAT)«E13TRP
'0171
TRAPPD=1
>0172
GOTO 14
>0173
15
RO=R
>0174
F130=E 1 3
>0175
14
N=N* 1
>0176
IF(N^(N/NPO)*NPO.NE.O) GOTO 11
>0177
13
XP=X«D1A
>0178
ELEV=Z«OIA .
>0179
?P=OEPTH-ELEV
>0160
SP-S°01 A
>0181
C
>0162
R=.308«SP
>0183
C
o u « o • ft « ft # ft ft o « a « » » « o « ft « a » a o « ft ft ft » ft o a ft ft ft ft ft ft ft ft « o ft « ft ft ft
10184
nlLN=.245»S<»E13
>0185
THETA=ACOSF(COSTH)«57.2958
>0186
IF(.NOT,TRAPP01 GOTO 7?
10187
C
fta«ftao»«ft»«««a0«aa«a««afte««»
JO 188
WRITE(6,101) T«SP»XP»ZP»B,ELEV,THETA
10189
C
a«««eft««»«ftftft««««attftft««efte«a«a«»»«ft«e«««a«»fteo««»
10190
GOTO 71
10191
C
JO 192
72
WRITE(6,101) T,SP,XP,ZP,8,ELEV,THETA,0ILN
>0193
101
FORMAT(6F10.2,F10.1,F10.4)
10194
C
-» e « a a e » « » a « • a « « « » « « « « « » « « « « « « « a « » « » « o « « « a a « « « « » « «
50195
71
IF I .NOT.FLAG) GOTO 11
38
-------
J0196
30197
00198
30199
30200
30201
30202
J 020 3
30204
30205
30206
30207
30208
30209
30210
30211
90212
90213
90214
0O21S
00216
00217
00218
00219
00220
00221
00222
00223.
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00236
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
PlJfcf'l M=100(2TRaP/0ePTm® 100.)
WRITE(6.65)
65 FORMAT("-LAST LINE ABOVE IS FOR MAXIMUM HEIGHT OF RISE.")
67 FORMAT{'STRAPPING LEVEL IS'SFlO.lt
• •' WITH DILUTION OF".F10.4)
TF(TRAPPO) WRITEI6.67) ZTRAP.SMTRAP
IF(.NOT.TRAPPO) WRITE(6.6B)
IF(TRAPPD)WRITE(6»69) PI1EPTH
69 FOHMAT (" HEIGHT OF RISE = ".F5. 1." PERCENT OF DEPTH")
68 FORMAT("ATRAPPING LEVEL NOT REACHED")
GOTO 20
C
C FIND NEXT STRATIFICATION AND RECOMPUTE LAST STEP IF
C NECESSARr
C
40 OS=DS»(ZLlrt-ZLAST)/(Z-ZLAST)
CALL SDERmS.E.R)
CHG0EN=1
Z=ZLAST
GOTO 45
41 CHGDEN=0
IPTS»IPTS*1
IF(IPTS.GT.NM) GOTO 42
G=DG(IPTS)
ZLIM=ZD(IPTS)
GOTO 43
42 WRITE(6»44)
44 FORMAT("OPLUME HITS SURFACE"/
«» HEI6HT OF RISE3100.0 PERCENT OF DEPTH")
flag=i
ROTO 43
74 CALL EXIT
END
SUBROUTINE SDERIV(StE.R>
COMMON G.FK »FM» COSTH.SINTHtCOSTHE•DS»C1»C2»E13»FLAG»GRAV
INTEGER FLAG
F13=0.
TF(E.LE.O.) GOTO 3
E13=E" < 1./3.)
COSTH=COSTHE*C1/(E13»E13)
IF(C05TH.LT.1•) GOTO 1
C
C THIS STOPPING CRITERIA IS BASED ON PLUME BECOMING HORIZONTAL
C AGAIN
C
FLAG=1
SINTH=0.
COSTH* 1 .
GOTO 3
1 SINTH=5QRT(1.-C0STH»C0STH)
3 FK=C2«S»R°SINTH«DS
FM=.109»S®E13*G«DS
RETURN
END
FUNCTION S1GMAT(SAL.T)
SIGO=<((6.8E-6#SAL)-4.B2E-4)®SAL*.8149)«SAL".093
R=1.E-6«T»(<.01667»T-.fll64)»T*18.03)
A*.001°T*((.0010843*T-.09818)»T*4.7867)
SUMT=(T-3.9fl)«(T-3.9fl)o(T*283.)/(503.57*(T*67.26))
SIGMAT=(SIGO*.l324)•{l.-A+B'MSIGb-.1324))-SUMT
RETURN
END
SUBROUTINE ASORT (DP»RMO» 'Kj PTS)
DIMENSION DP(1).RHO(l)
NESTED=NPTS
L=NESTED-1
39
-------
00262
00263
00264
00265
00266
00267
00268
00269
00270
-0021L
0027 S
00276
00277
no,10 M=1«L
nested=nested-i
no 10 I=l.NESTE0
IF(OP 11».LE.DP
DP 11 * 11=DUMHY
nUHMY=KHO(I)
RHO( Ii =RHO(I ~ 1)
JPHO(1*1)=DUMMY
GOTO 10
""CONTINUE"
RETURN
END
40
-------
SECTION 5
"OUTPLM" MODEL
THEORETICAL DEVELOPMENT AND MODEL DESCRIPTION
The computer model OUTPLM (18, 23) considers a plume element or puff. By
following this element, the characteristics of a continuous single plume in a
flowing ambient are described. The original model has been adapted for ocean
discharges. Density, velocity, temperature, and salinity are assumed to be
average properties of the element. The quantities of mass, momentum, and
energy are conserved. Two explicit functions for entrainment are used. An
equation of state is used to calculate the density of the ambient and the
plume element at each time step.
Initially, the plume segment has a horizontal velocity component (compu-
ter variable U) and a vertical component (V). It has temperature (T) and
salinity (S) and is injected into a uniformly flowing (UW) ambient. The
segment has a radius (B) and thickness (H) and is tracked in two dimensions,
depth (Z) and horizontal distance (X).
Entrainment brings ambient mass (plus momentum, temperature, and salin-
ity) into the plume element. Entrainment is assumed to consist of two mechan-
isms. One mechanism (EINS) is due to the impingement of current on the plume.
That mass which crosses the exposed area of the plume element projected by
currents is entrained. Since the plume element is assumed to be cylindrical,
two corrections to the projected area are included in EINS. They are a radius
correction and a correction for the curvature of the plume trajectory (see
(23)). The other mechanism (ZWEI) is aspiration entrainment which captures
-------
0.1 times the product of the external area of the plume element and its shear
velocity. The assumption that the total entrainment is the larger of these
two mechanisms gives good agreement with data.
In the computer program, the entrained mass (DM) is added to the ele-
ment's mass (PM) to become the new mass (SUM). The new temperature and salin-
ity of the element are the averages of the old values and the entrained ambi-
ent values (TA and SA) weighted by their relative masses. The new density
(DEN), and thus buoyancy, creates a vertical acceleration on the plume seg-
ment. Since the element is considered to be one of a train, each following
the preceding element, drag is assumed to be negligible and the vertical
velocity is incremented (by DV). The horizontal velocity is found by conser-
vation of horizontal momentum with the entrainment (DM) having the velocity of
the current. The segment length (H) is changed in proportion to the total
velocity (VEL) to conserve mass and pollutant. The radius (B) is changed to
correspond to the new mass and density. Dilution (DIL) is calculated by
comparing the initial volume to that of the element. The program terminates
execution when the vertical velocity reaches zero, the surface is reached, or
length scales or execution step limits are reached.
EXAMPLE INPUT, OUTPUT, AND MODEL LISTING
Input
Consider, for example, a horizontal, single port discharge into a flow-
ing, stratified ambient. The maximum discharge velocity is 2.0 m/sec. The
discharge port radius is 0.125 m. The temperature and salinity of the waste
are 15.0°C and 1.09 °/oo, respectively. The ambient at the surface is 15.0°C
and 27.2 °/oo. At 30 m, the level of the discharge, the ambient is 15.0°C and
42
-------
33.71 °/oo. Currents are assumed to be 0.10 m/sec. The input punch cards for
this case follow.
43
-------
¦¦jul
uw
T
15.00;
F10.3
F10.3 i
Si
1 .09!
F10.3
e :
125
F10.3 i
1.0
F10.3
A ;
o.io
F10.3
125
F10.3
minmmi
riTTTTTTTTT
V - Vertical velocity of the plume, m/sec.
UW- Current speed, m/sec.
T - Temperature of the plume, °C.
S - Salinity of the plume, o/oo.
H - Thickness of the plume element C^B), m.
E - Impingement entrai.nment coefficient, (recommended =
A - Aspiration entrainment coefficient (recommended = 0
B - Plume radius, m.
1.0)
.1)
5 moo
p x co cr
I n> -s —I
* 3 CL -a
i "a r—
¦1
13
"a
*
beamu
J.
-------
u
KSCALE
YSCALE
NPTS
HERB
] |
tR| |
2.0
500
50C
2
j
1000 i 25
1
! F10.3 1 F5.0
F5.C
15
15 i 15
\
\
!
{
i
*
1
?
vi
i! j
TTCT
1
i _
I
i
1
!
.-ULillL
| |
Mill
1 | i 1 ! 1 M ; i! ;
1 1 : i 1 i ! I i i ! f
Hi
-i I ;
°j ! !
Ml
O
. i ¦ rj;=
1 !oi
¦H
J . i
i is
1 j i -• ^
j... ®!.. i j ;c. 4
cn
U - Horizontal Velocity of the plume.
XSCALE - Number of diameters in the X direction at which run will be cut-off.
YSCALE - Number of diameters in the Z direction at which run will be cut-off.
NPTS - Number of points at which ambient is specified.
ITERB - Number of steps allowed.
IR - Interval of printout.
mOO
x o> cr
o> —I
3 Q- -o
r~
-T03
CD
3
-o
C
c+
-------
DP (I)
SACI) :
TA r :
27.20 <
15.00
Card 3
33.71 i
5.00
Card 4
F10.3
F10.3
F10.3 i
DP(I) - Depth where ambient is specified, m.
SA(I) - Salinity, o/od.
TA(I) - Temperature, °C.
m
C
O
X
a>
c=
0J
-s
-H
3
a.
"U
TD
i—
OJ
3
a>
t—i
a>
3
r+
XJ
O
c
•
c+
-------
Output
Model output and program listing follow. The line above the header
echoes the first two input cards. Ambient stratification is then given with
input depth, salinity, temperature, and computed sigma-t densities (density,
kg/m3, —1000). Other discharge conditions are given. K is a velocity
ratio, = U/UW. Q is the volume rate of discharge. The first data line in the
listing are the initial conditions. The horizontal distance in the direction
of the jet discharge is the X dimension. Z is depth measured from the water's
surface. B is the radius, THICK is the plume element thickness along the
centerline. MASS is the mass of the plume element being considered. EINS and
ZWEI are the previously mentioned entrainments. DENDIF is the plume's density
difference with the ambient. DILUTION is an average volumetric dilution, the
inverse of CONC. TEMPDIF is the temperature difference between the plume
element and the ambient at that level.
47
-------
0 .100 15.000 1.090 .125 1.000 .100 .125 2.000 500 500 2 1000 25
• ••»••••• •••••••••••••••• • • » «
OUTFALL BUOYANT (JET) PLUME IN FLOWING. STRATIFIED AMBIENT
0S3 VERSION 2.3 5-16-79
AMBIENT STRATIFICATION- DEPTH,M SALIN TEMP.C SIGMaT
0 27.20 15.00 20.01
30.00 33.71 15.00 25.00
K FROUOE (i CURRENT
2.00E 01 8.O0E 00 9.82E-02 1.00E-01
MODEL INPUT (LINE 1) AND MODEL OUTPUT
x Z B THICK MASS EINS ZWEI DILUTION DENDIFF HOR VEL VER VEL TOT VEL TEMPO IF
OE 00 3.00E 01 1.25E-01 1.25E-01 6.14E 00 1.12E-02 1.91E Ul 1.00E 00 2.50E 01 2.00E 00 OE 00 2.00E 00 OE 00
4.44E-03 3.00E 01 1.25E-01 1.24E-01 6.18E 00 1.12E-02 1.91E 01 1.01E 00 2.48E 01 1.99E 00 5.43E-04 1.99E 00-2.33E-10
1.21E-01 3.00E 01 1.47E-01 1.06E-01 7.30E 00 2.97E-03 4.97E-02 1.18E 00 2.10E 01 1.70E 00 1.36E-02 1.70E 00 2.33E-10
2.67E-01 3.00E 01 1.74E-01 9.02E-02 8.68E 00 4.27E-03 5.91E-02 1.40E 00 1.77E 01 1.44E 00 2.76E-02 1.44E 00 OE 00
4.41E-01 3.00E 01 2.05E-01 7.69E-02 1.03E 01 6.21E-03 7.03E-02 1.66E 00 1.49E 01 1.23E 00 4.22E-02 1.23E 00 1.40E-09
CO 6.50E-01 3.00E 01 2.42E-01 6.57E-02 1.23E 01 9.13E-03 8.36E-02 1.98E 00 1.25E 01 1.05E 00 5.7SE-02 1.05E 00 9.31E-10
8.99E-01 3.00E 01 2.84E-01 5.64E-02 1.46E 01 1.36E-02 9.95E-02 2.34E 00 1.05E 01 8.99E-01 7.49E-02 9.02E-01 2.33E-10
1.20E 00 2.99E 01 3.34E-01 4.86E-02 1.74E 01 2.05E-02 1.18E-01 2.78E 00 8.83E 00 7.72E-01 9.38E-02 7.77E-01 2.33E-10
1.55E 00 2.99E 01 3.90E-01 4.22E-02 2.06E 01 3.13E-02 1.41E-01 3.31E 00 7.42E 00 6.65E-01 1.15E-01 6.75E-01 6.98E-10
1.97E 00 2.98E 01 4.54E-01 3.70E-02 2.45E 01 4.78E-02 1.68E-01 3.93E 00 6.23E 00 5.75E-01 1.37E-01 5.91E-01 9.31E-10
2.47E 00 2.97E 01 5.26E-01 3.28E-02 2.92E 01 7.25E-02 1.99E-01 4.66E 00 5.22E 00 4.99E-01 1.62E-01 5.25E-01 6.98E-10
3.03E 00 2.94E 01 6.03E-01 2.96E-02 3.47E 01 1.08E-01 2.37E-01 5.54E 00 4.36E 00 4.36E-01 1.87E-01 4.74E-01 9.31E-10
3.65E 00 2.91E 01 6.86E-01 2.73E-02 4.13E 01 1.56E-01 2.83E-01 6.59E 00 3.62E 00 3.82E-01 2.11E-01 4.37E-01 1.16E-09
4.32E 00 2.87E 01 7.72E-01 2.56E-02 4.91E 01 2.18E-01 3.37E-01 7.83E 00 2.98E 00 3.38E-01 2.31E-01 4.09E-01 1.40E-09
5.02E 00 2.82E 01 8.64E-01 2.43E-02 5.84E 01 2.96E-01 4.01E-01 9.31E 00 2.43E 00 3.00E-01 2.47E-01 3.88E-01 1.40E-09
5.76E 00 2.75E 01 9.64E-01 2.32E-02 6.94E 01 3.92E-01 4.77E-01 1.11E 01 1.94E 00 2.68E-01 2.57E-01 3.71E-01 4.66E-10
6.52E 00 2.68E 01 1.07E 00 2.22E-02 8.26E 01 5.11E-01 5.67E-01 1.32E 01 1.51E 00 2.41E-01 2.62E-01 3.56E-0X 6.98E-10
7.32E 00 2.58E 01 1.20E 00 2.12E-02 9.82E 01 6.62E-01 6.74E-01 1.56E 01 1.13E 00 2.19E-01 2.60E-01 3.40E-01 6.98E-10
8.15E 00 2.48E 01 1.34E 00 2.01E-02 1.17E 02 8.06E-01 7.42E-01 1 .86E 01 7.98E-01 2.00E-01 2.52E-01 3.22E-01 6.98E-10
8.97E 00 2.38E 01 1.52E 00 1.87E-02 1.39E 02 9.59E-01 7.75E-01 2.2IE 01 5.11E-01 1.84E-01 2.36E-01 2.99E-01 1.16E-09
9.79E 00 2.27E 01 1.73E 00 1.71E-02 1.65E 02 1.14E 00 7.84E-0) 2.63E 01 2.69E-01 1.71E-01 2.14E-01 2.74E-01 9.31E-10
1.06E 01 2.17E 01 1.99E 00 1.54E-02 1.96E 02 1.36E 00 7.57E-01 3.13E 01 7.07E-02 1.59E-01 1.87E-01 2.46E-01 2.33E-10
1.15E 01 2.08E 01 2.31E 00 1.36E-02 2.33E 02 1.61E 00 6.82E-01 3.72E 01-8.51E-02 1.50E-01 1.57E-01 2.17E-01 2.33E-10
1.23E 01 1.99E 01 2.70E 00 1.18E-02 2.78E 02 1.92E 00 5.61E-01 4.42E 01-1.97E-01 1.42E-01 1.24E-01 1.88E-01-b.98E-10
I.32E 01 1.93E 01 3.I6E 00 1 .02E-02 3.30E 02 2.28E 00 4.24E-01 5.26E 01-2.66E-01 1.35E-01 9.14E-02 1.63E-01-1.40E-09
1.41E 01 1.88E 01 3.69E 00 8.93E-03 3.93E 02 2.71E 00 3.19E-01 6.25E 01-3.00E-01 1.30E-01 6.01E-02 1.43E-01 4.66E-10
1.52E 01 1.84E 01 4.26E 00 7.98E-03 4.67E 02 3.23E 00 3.07E-01 7.44E 01-3.09E-U1 1.25E-01 2.64E-02 1.28E-01 9.31E-10
1.61E 01 1.83E 01 4.49E 00 7.71E-03 5.04E 02 3.48E 00 6.63E-01 7.97E 01-3.04E-01 1.23E-01-2.07E-03 1.23E-01 2.33E-10
NUMBER OF STEPS= 636
-------
00001
00002
00003
00004
00005
0000b
U0007
000UB
oooov
00010
00011
00012
00013
00014
00015
0001b
0001 7
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
0003S
00036
00037
00038
00034
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
0005H
00059
00060
00061
00 062
00063
PROGRAM OUTPLM
DI MENS I ON DP(30) ,SA(30) ,TA(30).DENP(30)
DATA (G=9.«) , < I W0=2. I ,
6 FORMAT(BF10.3*/,IF]0,3*2F5.0,3151
9 FORMAT(3F10,3)
14 FORMAT 126X,4F6.2/)
READ (10,61 V,UW,T,S»H,E,A,B,U,XSCALE»ZSCALE,NPTS,ITERB, IR
WRITE(11*51 V*UW)
8 CONTINUE
WRITE (11,14) ((OP(I)iSAtl)»TA(I),D£NP(1)> ,1=1,NPTS)
DT =OD=DB=1.
X=BSAVE=ZWEI = ZERO
Z=DPINPTSI
OEL = T-T A(NPTS)
AK = 99<59.
DIL = 1 •
BA=B
BO=B«TWO
VI =VEL=SUR7(V»V*U»U)
OENA=1000.«SIGMAT(SA(NPTS),TA(NPTS))
DEN=1000.*SIGMAT(S*T)
RHO=OEN
0EN01FF=DENA-DEN
PM=PMO=PI»H«R»H»OEN
(J=HM/DEN»YEL/H
IF(UW.NE.ZERO) AK=VEL/UW
FR=\/EL/SURT (DENl)lFF/OEN*TWO»B»G)
7 FORMAT!/"' "/
1" K FROUDE 0 CURRENT"/1X4E9.2//
2" MODEL INPUT (LINF II ANO MODEL OUTPUT X
3 Z B THICK MASS EINS ZWEI DILUTION DENDIFF
4 HOrt VEL VEH VEL TOT VEL TEMPOIF")
WRITE(11,7) AK
ZWEI=TW0«P1»DENA«H«B»ABSIVEL-UW*U/VEL)«A«0T
DM=E1NS
IF(ZWEI.GT.EINS) DM=ZWEI
GAMMA=OM/.0069555/PM
DT=DT/GAMMA
EINS=EINS/GAMMA
49
-------
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
0007b
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
oooaa
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00103
00109
00110
00111
00112
00113
OUl 14
00115
dm=dm/gamha
13 FORMAT (1X.14E9.2)
IFtJ.EO.1)WHITE(11t13)X.Z.Pt.H,PM,EIN5tZWEI
1•OIL tDENOIFF.U.V»VEL«DEL
SUH=PM*DM
UO=U
U= /IH.NE.1)G0 TO 99
IF (RATIOZ.GT.ZSCALE.OR.RATIOX.GT.XSCALE)GO TO 999
98 RATIOX=X/BO
WATI0Z=Z/dO
DIL = (RHO»PM)/
-------
SECTION 6
"DKHPLM" MODEL
THEORETICAL DEVELOPMENT
The computer model DKHPLM (7, 21, 24) is an approach to the problem of
merging plumes which considers the detailed dynamics of the process instead of
simplifying the problem to a idealized slot plume or to a combination of plume
and slot plume. DKHPLM considers three zones of plume behavior; zones of flow
establishment, established flow, and merging. The first two zones are based
on the analyses of Hirst (25, 26) for a plume in a stratified, flowing en-
vironment. In the zone of merging, neighboring plumes are superimposed. This
allows a smooth transition as single plumes begin to compete for dilution
water then gradually merge with their neighbors. Equations for the conserva-
tion of mass, pollutant, energy, and momentum are developed for all three
zones. Entrainment is an explicit function dependent on the local Froude
number, plume spacing, excess velocity, and ambient velocity. Similar lateral
profiles, a 3/2 power approximation of a Gaussian, are assumed for velocity,
concentration, and temperature. These profiles are superimposed in the merg-
ing zone. A complete theoretical development of this model is beyond the
scope of this paper but can be found in (21, 24, 25, and 26). The following
is a brief summary.
Zone of Flow Establishment
All quantities are assumed uniformly distributed in the plume at the
point of discharge. In the zone of flow establishment, these uniform profiles
-------
change to similar profiles as the boundary layer diffuses inward to the cen-
terline of the jet. The rate at which the profiles of velocity, concentration
and temperature develop may vary. The integrated forms of the governing
equations are:
conservation of mass,
d/ds /°° V r dr = -Tim r (r v) = E (32)
o
conservation of energy,
d/ds J°° V (T-T ) r dr = - d T /ds J°° V r dr - lim r -» » (r v'T') (33)
O . 00 O
conservation of pollutant,
d/ds f° V (C-C ) r dr = -d C /ds J°° V r dr - lim r -> °° (r v' C') (34)
O 00 00 O
conservation of momentum in the s equation,
d/ds f° r dr = U E sinBj cos02 + J°° g (p - p)/pH) r dr sin02
o o 00 Q
- lim r a (r u' v') (35)
where 0X is the horizontal angle between the centerline and the x axis and 02
is the angle between the centerline and the horizontal. Two.additional inte-
gral equations have been developed from equation (35) to describe momentum in
two additional plume coordinants. These "natural" coordinants of the plume,
described in (25), are converted to conventional three-dimensional Cartesian
coordinates for model output. Implicit in the derivation of.these equations
are the assumptions that:
a) flow is steady in the mean,
i ¦ >
b) flow is fully turbulent,
c) fluid is incompressible and density variations are included only in
the buoyancy terms,
d) all other fluid properties are constant,
e) no frictional heating,
f) pressure variations are purely hydrostatic,
52
-------
g) changes in density are small enough to be approximated by a linear
equation of state,
h) flow within the jet is axisymmetric,
i) flow within the jet can be approximated as boundary layer flow,
j) the ambient is infinite in extent.
Several of the assumptions are compensated for in the solution. The zone of
flow establishment uses a special entrainment function (see equation (128) of
(24)) which is a function of local Froude number, velocities, port diameter,
spacing, and thickness of the developed flow region.
Zone of Established Flow
The equations of mass, energy, pollutant, and momentum (equations 32-35),
and the additional momentum equations mentioned above are also solved in the
zone of established flow, but are cast into a slightly different form. The
governing equations are written in a cylindrical coordinant system where is
the circumferential angle around the plume and cross section and the independ-
ent variables are r and s. These are evaluated using the assumed 3/2 power
approximation to Gaussian lateral profiles. The angle is the angle between
the centerline projected to the xy plane and the x axis, 02 is the angle be-
tween the centerline and the xy plane. These angles relate the two coordinate
systems. Another entrainment function is used in the zone of established flow
which is a function of the local Froude number, velocities, plume diameter,
and spacing.
Zone of Merging
When adjacent plumes begin to overlap, the discharge is no longer consid-
ered axisymmetric. The distributions of plume properties are superimposed.
-------
Another entrainment function is used which also considers the variable en-
trapment surface during merging. A drag term is also introduced to account
for the additional bending of the plumes after merging.
MODEL DESCRIPTION
An equation of state for sea water has been added to the original program
(subroutine SIGMAT) replacing the linear equation of state. In the zone of
flow establishment, the system of six governing equations are solved simultan-
eously (subroutine SIMQ) and stepped forward in space by a Hamming's modified
predictor-corrector method (subroutine HPCG). This procedure continues until
velocity, temperature and concentration become fully developed. Subroutine
0UTP1 contains the results which are stored as initial conditions for the zone
of established flow.
In the zone of established flow, similar profiles and the integral method
allow solution of the six governing equations for the six unknowns initially
by Runge-Kutta integral approximation and then by the Hamming's modified pre-
dictor-corrector method. At the point where the plumes overlap, the assumed
similarity no longer applies. The merging plumes have axes of symmetry along
the discharge line and normal to it. Only one quadrant of a plume, taken to a
midpoint of the overlap area, is evaluated. The profiles are superimposed in
this region using integral similarity coefficients. Above the point where the
plume and the ambient have equal density, results are obtained by extrapola-
tion.
The program itself contains many comments and explanations which serve as
further documentation.
54
-------
EXAMPLE INPUT, OUTPUT, AND MODEL LISTING
Input
As an example case, suppose a diffuser 1250 m long discharges 12.5 m3/sec.
The ports are 0.178 m in diameter, oriented vertically, and spaced 5 m apart.
The ambient currents are assumed to be 0.05 m/sec and normal to the discharge
line. (Note that the angle of currents to the diffuser should be within about
20° of normal.) The velocity through the discharge ports is 2.0 m/sec. The
effluent is 15.01°C and 1.09 °/oo. (Note that the effluent should have some
finite temperature difference with the ambient.) The ambient at the level of
discharge is 15.0°C and 33.71 °/oo (1025 kg/m3). The ambient temperature
gradient is negligible but must be finite. The ambient salinity gradient is
-0.026 °/oo per m (a density gradient of 0.020 kg/m3/m). Depth of the dis-
charge is 50 m. Computer punch card format for this input follows.
55
-------
HEADER INFORMATION (Nil)
0 DKHPLM EXAMPLE RUN
J3 - Number of runs to follow C0=one run).
-------
-.02600 !
15.01
1.09
15.00
33.71
.00001
F10.6
F10.6
F10.6 i
F10.6
F10.6
F10.6
F10.6
F10.6i
UO - Initial jet velocity, m/sec.
UA - Current speed, m/sec.
TO - Initial jet temperature, °C.
CO - Initial jet salinity, o/oo.
TAO- Ambient temperature at the level of discharge, °C,
CAO- Ambient salinity at the level of discahrge, °C.
TIZ- Temperature gradient, °C/m.
CIZ- Salinity gradient, o/oo/m.
m
o
X
Q>
Q)
-s ^
3
Q. "O
-o
ro 2
(D
HH
3
"O
c
c+
-------
en
CO
r *- ¦
i
konecnr^
j
DO j
i
!'
TH1 |
]
i
TH2
i
1
SF !
SPACE
H
I
0.178
i
i
90.0 j
90.0 L
!
500.0 \
5.0
50.0
F10.6
1
l
F10.6 j
F10.6 !
i .
I
F10.6 j
F10.6
F10.6
5
SPACE
Port diameter, m.
Horizontal discharge angle (DEG) relative to the current (70
-------
Output
Example output and program listing follow. The initial input conditions
are printed along with dimensionless conditions. The first section of data
listing concerns the zone of flow establishment. The columns list: length
along the plume axis (S), horizontal distance parallel with discharge line
(X), horizontal distance normal to the discharge line (Y), vertical distance
from the level discharge (Z), the horizontal angle from the plume's axis to
the discharge line (TH1), the angle the plume's axis makes from the horizontal
(TH2), the radius (B), potential core widths for velocity, temperature and
concentration (RU, RT, and RC), normalized centerline disparities of velocity,
temperature, and concentration with the ambient (DUCL, DTCL, and DCCL), the
ambient density normalized by the density of the discharge (PI), and non-
dimensional time (TIME). In the zone of established flow, width and average
dilution (Q/Qo) are also given. Time is given in seconds. The centerline
dilution before zone of merging is:
Sm = 0.52 Q/Qo (36)
m
after complete merging (which never actually occurs):
Sm = 0.70 Q/Qo (37)
based on the assumed distributions of concentration and velocity.
If stagnant conditions are specified, the dilution should be taken at the
point of buoyant equilibrium with the ambient. The location is announced in
the model output.
59
-------
SOLUTION TO MULTIPLE BUOYANT DISCHARGE PROBLEM WITH AMBIENT CURRENTS AND VERTICAL GRADIENTS
DKHPLM EXAMPLE CASE I.D. = 0
DISCHARGE VELOCITY = 2.00-M/S o TEMP. = 1^,.01-DEG C »» SALINITY = 1.09-PPT
niAMETEW = . lfl-M •» SPACING = 5.00-M «® DEPTH = 50.00-M
AMBIENT CONDITIONS AT DISCHARGE ELEVATION. VELOCITY = ,05-M/S »« TEMP. = 15.O0-DEG C »• SALINITY = 33.71 PPT
AMBIENT STRATIFICATION GRADIENTS »« TEMPERATURE .000010-DEG C/M »« SALINITY •» -0.026000-PPT/M
AMHIENT CONDITIONS AT DISCHARGE - NONDIMENS IONAL
TFMP= .99933 SAL INITY = 30.92661 DENS IT Y= 1.02501 T EMPGR = .000000 SALGR = -0.002123 DENGR=-0.0000018
PORT SPACING L/D = 28.09
FROUDE NO = 9.58 VELOCITY RATIO = .025 STRATIFICATION NO = 14U96.6
ZONE OF
FLOW ESTABLISHMENT
ALL LENGTHS
ARE IN
METERS
S X
Y
Z
THl
TH2
B
RU
RT
RC
DUCL DTCL DCCL
PI
TIME
n o
0
n
90.00
90.00
0
.089
.089
.089
1.000 1.000 1.000
1.02501
0
.09 .00
.00
.09
90.00
99.91
.018
.083
.083
.083
1.000 1.000 1.000
1 .02501
1.0
•in .oo
.00
. l 8
90.00
89.81
.038
.077
.077
.077
1.000 1.000 1.000
1.02501
2.0
.27 .00
.00
.27
90.00
89.70
.059
.070
.070
.070
1.000 1.000 1.000
1.02501
3.0
.3ft .00
.00
. 36
90.00
89.57
.082
.063
.062
.062
1.000 1.000 1.000
1.02500
4.0
.44 .00
.00
. 44
90. 00
89.43
.107
.054
.053
.053
1.000 1.000 1.000
1.02500
5.0
.53 .00
. 0 0
.53
90.00
89.28
.133
.045
.043
.043
1.000 .999 1.000
1.02500
6.0
.6? .00
.00
.62
90.00
89. 1 1
.162
.035
.033
.033
1.000 .999 1.000
1.02500
7.0
.7 1 -0.00
.01
. n
90.00
88.92
.193
.024
.021
.021
1.000 .999 .999
1.02500
8.0
.80 -0.00
.01
.80
90.00
88.72
.226
.011
.00 7
.007
.999 .999 .999
1.02499
9.0
STARTING LFNGTH,
T =
.844
.84 -0.00
.01
.84
90.00
88.62
.243
.004
.999 .999 .999
1.02499
.4
STARTING LFNGTH,
VELOCITY
=
.066
.87 -n.oo
.0 1
.87
90.00
88.57
.252 -
0.000
.999 .930 .930
1 .02499
9.7
ZONE OF
ESTABLISHED FLOW
S X
Y
Z
THl
TH2
WIDTH
DUCL
OTCL
DCCL
PI TIME
Q/UO
.87 -0.00
.01
.87
90.00
88.57
.51
.997
.930
.930
1.02499 .43
2.07
1.5ft .00
.05
1 .58
90.00
P4.59
1 .24
.429
.354
.354
1.02498 1.00
5.48
2.29 .00
. 14
2.2a
90.00
80.81
2.03
.284
. 196
.196
1.02497 2.03
9.99
3.00 .00
.2*
2.98
90.00
77.58
2.82
.223
.126
.127
1.02495 3.43
15.67
3.71 .00
.45
3.67
90.00
74.86.
3.59
.190
.088
.089
1.02494 5.12
22.48
4.43 .00
.65
4.36
90.00
72.50
4.35
.169
.065
• 066
1.02492 7.04
30.37
ZONE OF MERGING
FLOW
5.14 .00
.87
5. 03
90.00
69.97
5.09
.154
.050
.051
1.02491 9.14
39.29
5.65 .00
1.14
5.69
90.00
66.86
5.71
.144
.040
.041
1.02490 11.40
47.81
6.56 .00
1.43
6.34
90.00
64 . 78
6.27
.138
.033
.034
1.02488 13.76
55.90
7.27 .00
1 .74
6.98
90.00
63.27
6.81
.132
.027
.028
1.02487 16.20
63.73
8.01 .00
2.08
7.63
90.00
62.08
7.35
.126
.022
.024
1.02486 18.81
71.60
9.43 .00
2. 77
8. 88
90.00
60.32
8.44
.113
.015
.017
1.02483 24.21
86.36
10.86 .00
3.49
10.11
90. 00
58.89
9.64
.099
.010
.011
1.02481 30.19
100.48
-------
DkmPLM EXAMPLE CASE
1.0. a
0
s
X
Y
2
TH1
TH2
WIDTH
DUCL
DTCL
DCCL
PI
TIME
Q/O
12.2*
.00
4.24
11.3?
90.00
57.55
11.18
.087
.006
.008
1.02479
36.87
114.06
13.75
.00
5.04
12.55
90.00
56.06
12.63
.084
.004
.005
1.02476
44.28
127.61
15.17
.00
5.85
13.72
90.00
54.44
14.13
.080
.002
.004
1.02474
51.69
140.44
16.60
.00
6.70
14.86
90.00
52.58
15.73
.076
-0.000
.002
1.02471
59.38
152.93
in.02
.00
7.58
15.98
90.00
50.38
17.43
.071
-0.002
.000
1.024ts9
67.41
165.00
PLUMF HAS
REACHED
ITS EQUILIBRIUM
HEIGHT -
STRATIFIED
ENVIRONMENT
ia.«2
.00
8.10
16.59
90.00
48.94
18.45
.068
-0.003
-0.000
1.02468
72.11
171.58
19.00
.00
8.22
16.72
90.00
48.60
18.68
.067
-0.003
-0.001
1.02468
73. 17
173.02
19. lfl
.00
8.34
16.86
'90.00
48.25
18.91
.066
-0.003
-0.0 0 1
1.02467
74.24
174.45
19.36
.00
8.46
16.99
90.00
47.89
19.14
.065
-0.004
-0.001
1.02467
75.32
175.88
19.53
.00
8.58
17.12
90.00
47.52
19.38
.065
-0.004
-0.001
1.02467
76.40
177.29
19.71
.00
8.70
17.25
90.00
47.15
19.62
.064
-0.004
-0.001
1.02467
77.50
178.69
19.89
.00
8.82
17.38
90.00
46. 7b
19.86
.063
-0.004
-0.001
1.02466
78.60
180.08
20.07
.00
8.94
17.51
90.00
46 . 36
20.11
.062
-0.004
-0.002
1.02466
79.72
181.46
20.25
.00
9.06
17.64
90.00
45.95
20.35
. 062
-0.005
-0.002
1.02466
80.84
182.84
20.4 2
.00
9.19
17.77
90.00
45.53
20.60
.061
rO.005
-0.002
1.02466
61.97
184.20
20.6 0
.00
9.31
17.89
90.00
45.09
20.85
.060
-0.005
-0.002
1.02465
83.11
185.55
20.7B
.00
9.44
18.0?
90.00
44.65
21.10
.059
-0.005
-0.002
1.02465
84.26
186.89
20 . <36
.00
9.57
18.14
90.00
44. 19
21.36
.058
-0.005
-0.002
1.02465
85.42
188.21
21.14
.00
9.70
18.27
90.00
43.72
21.62
.058
-0.005
-0.003
1.02465
86.59
189.53
21.31
iOO
9.H 2
18.39
90.00
4 3.23
21.87
.057
-0.006
-0.003
1.02464
87.77
190.83
21 .49
.00
9.95
18.51
90.00
42.73
22.14
.056
-0.006
-0.003
1.02464
88.96
192.12
21 .67
.00
10.09
IB.63
90.00
42.21
22.40
.055
-0.006
-0.003
1.02464
90. 17
193.39
21.05
.00
10.22
18.75
90.00
41 .68
22.67
.054
-0.006
-0.003
1 .02464
91.38
194.65
22.03
.00
10.35
18.87
90.00
41.13
22.93
.053
—0.006
-0.003
1.02463
92.61
195.90
22.20
.00
10.49
18.98
90.00
40.57
23.20
.052
-0.007
-0.003
1 .02463
93.85
197.13
22.47
.00
1 0.69
19.16
90.00
39.68
23.60
.051
-0.007
-0.004
1.02463
95.73
198.95
22.83
.00
10.97
19.38
90.00
38. 44
24.15
.049
-0.007
-0.004
1.02462
98.29
201.32
23. 1R
.00
11.25
19.60
90.00
37. 12
24.70
.04 7
-0.007
-0.004
1.02462
100.91
203.62
23.54
.00
11.53
19.81
90.00
35.70
25.27
.045
-0.008
-0.004
1.02462
103.58
205.85
23.89
.00
11.83
20.01
90.00
34. 19
25.84
. 044
-0.008
-0.005
1.02461
106.32
20S.00
24.25
.00
12.12
20.21
90.00
32.58
26.41
.042
-0.008
-0.005
1.02461
109.13
210.06
24,61
.00
12.43
20.40
90.00
30.86
26.98
. 040
-0.009
-0.005
1.02460
112.01
212.02
24.96
.00
12.73
20.57
90.00
29.02
27.55
. 038
-0.009
-0*006
1.02460
114.96
213.89
25.32
.00
13.05
20.74
90.00
27.05
28.12
.036
-0.009
-0.006
1.02460
118.00
215.66
25.67
.00
13.37
20.90
90.00
24.96
28.67
.034
-0.010
-0.006
1.02459
121.11
217.30
26.03
.00
13.69
21 .04
90.00
22.73
29.21
.032
-0.010
-0.006
1.02459
124.30
218.83
26.39
.00
14.03
21.17
90.00
20.36
29.72
.030
-0.010
-0.006
1.02459
127.58
220.22
26.74
.00
14.36
21.29
90.00
17.85
30.20
.029
-0.010
-0.007
1.02459
130.95
221.47
27. 10
.00
14.70
21.39
90.00
15.21
30.64
.027
-0.011
-0.007
1.02458
134.39
222.58
27.45
.00
15.05
21.48
90.00
12.45
31.04
.026
-0.011
-0.007
1.02458
137.91
223.52
DKHPLM example case
I.D.
= 0
S
X
r
Z
THl
TH2
HIDTH .
DUCL
DTCL
DCCL
PI
TIME
Q/O
27.81
.00
15.40
21.54
90.00
9.57
31 .36
. 025
-0.011
-0.007
1.02453
141.49
224.31
28.34
.00
15.93
21.61
90.00
5.09
31.71
.023
-0.011
-0.007
1 .02458
146.96
225.16
28.88
.00
16.46
21.64
90.00
.50
31 .88
. 023
-0.011
-0.007
1 .02458
152.51
225.63
29.06
.00
16.64
21.64
90.00
-1.03
31.91
.023
-0.011
-0.007
1.02458
154.36
225.72
PLOMF HAS REACHED ITS MAXIMUM HEIfiHT - STRATIFIED environment
MAXIMUM Pl.lJMF HEIGHT. EXTRAPOLATED LINEAR = 44.52 PARACOLIC = 43.41
NO OF INTEGRATION STEPS= 898 NO OF HALVINGS = 11 ABSERR= .00100 PRMT<5>= 1.000 FINAL SPACE = 4.000
-------
00001
00002
00003
00004
00005
onoob
D0007
00008
00009
00010
000 11
0001?
00013
00014
II001S
000 If.
00017
000 1H
000 1*
00020
00021
00022
00023
00024
00026
00026
00027
(>0 0?ft
000 29
1)0030
90031
00032
00033
00034
00035
00036
J 0 0 3 7
>0038
00 0 39
0 004 0
10041
00042
.10043
1)0044
11004S
00046
00047
.') 0 0 4 R
;>0049
00 0 SO
00 051
DOO'52
0 0 0 b 3
i)OOS4
0 0 055
00056
00057
U0058
00059
•J 0 0 6 0
)OOt>l
00062
00063
PROGRAM DKHPLM
THIS PROGRAM CALCULATES THE CHARACTFRISTICS OF A LINE UF EOUALLY
SPACFD BUOYANT DISCHARGES INTO A FLOWING STRATIFIED AMBIENT
IT IS ASSUMED THAT BUOYANCY IS DUE TO TEMPERATURE ANL) SALINITY
THE METHOD OF SOL°N INVOLVES 7 ORDINARY DIFFERENTIAL EQUATIONS
WHICH ARE 1. CONSERVATION OF MASS! ?. CONSERVATIONS OF ENERGY?
3. CONSERVATION OF CONCENTRATION (SALINITY)! 4. DENSITY DEFICIENCY
(DEPTVFn FROM 2)3)5 5.6.7. MOMENTUM F'JNS IN THE Z (AXIAL).
K (VERTICAL) ) J(HORIZONTAL AND PARALLEL TO AMBIENT CURRENT)'01MECT1 ON
EON 4 IS USED ONLY IN ALGEBRAIC FORM.
The SIX EQUATIONS are WRITTEN FOR A CONTROL VOLUME WHICH IS FINITE
IN A DIRECTION PERPENDICULAR TO THE JET AND INFINITESIMAL IN THE
DIRECTION OF THE JET AXIS.
««»«»« INPUT VARIABLES °09«»o»
INPUT CONSISTS OF THREE CARDS FOR EACH CASE RUN - SETS OF .THREE
CARDS CAN BE STACKED FOR MULTIPLE RUNS
»»» CARD NO. 1
The FIRST INPUT IS AN INDICATOR AS TO THE NUMBER OF RUNS
REMAINING AFTER PRESENT PUN. ZERO OR BLANK INDICATES LAST
RUN OR A SINGLE RUN. FORMAT (14)
THE REST OF THIS CAPO CONTAINS ANY INFORMATION THE USER WISHES
TO HAVE PRINTED OUT AT Thf TOP OF EACH PAGE.
1.
1.
2.
3.
4.
5.
6.
7.
H.
1.
2.
CARD NO.
FORMAT (8F10.0)
DISCHARGE VELOCITY (M/S)
AMBIENT VELOCITY -
AMBIENT SALINITY AT DISCHARGE LEVEL IPPT)
LINEAR AMBIENT TEMPERATURE GRADIENT (DEG C/M)
LINEAR AMBIENT SALINITY GRADIENT (PPT/M)
~«» CARD NO. 3 «»
FORMAT
(6F10.0)
DISCHARGE DIAMETER (M)
HORIZONTAL 01SCHAriGF ANGLE (DEG) RELATIVE TO THE AMBIENT CURRENT
90 DEG IS FOR CURRENT NORMAL TO LINE OF DISCHARGE PORTS
MODEL IS NOT GOOD FOR LARGE VARIATIONS FROM THIS DIRECTION
DISCHARGE ANGLE (DEG) RELATIVE TO THE HORIZONTAL.
90 DEG IS VERTICAL - ZF.RO IS HORIZONTAL DOWN STREAM -
DISTANCE ALONG PLUMP CENTFRLINF SIMULATION IS TO STOP (M)
DISTANCE BETWEEN PORTS (M) .
USE A LARGE VALUE TO SIMULATE SINGLE PORT DISCHARGES
DEPTH UF DISCHARGE (M)
EXTERNAL DF.RIV.OUTP
PFAL I I 1.112.I 13.I 14
DIMENSION PRMT(S). Al IX(16.6) . F(6). FP(6) .Nll(lH)
COMMON /AMH/ TIZ.CIZ.CIO.TIO.PIO.GA.AL.RE.FR.Rl.FRL
COMMON /OUT/ L.TH1.Th2.H,Z.SNEW.OS.R.K.N11.J3.J2.J5.JJ.SPACE.S1GMA
X•DO.CD.A7.XXX
COMMON /PLTE/ Jb.J7.J8.J9.J10.XlM.YlM.X2M.Y2M.Y3M
62
-------
000& E>KovaJc)
00073
T T 3 = 0.4S«1.07034 « OO O ^"?
I!4 = 0.31558*1.094 1
00074
00075
1
oFAD (5.102) J.I,Nil
0007b
TF(EOF(5).NE.O) GO TO 221
00077
.J5 = 0
00078
RE = 1.4253E-4
00079
GA=-7.6669E-4
00080
cn=2-
00081
ns=4.
00082
Ml = l
00083
jj = l
00084
.11=0
00085
J2 = l
00086
,J6=0
00087
11 = 1
00088
M = 16
00084
C Nil CONTAINS ALPHAMERIC INFO WHICH IS PRINTED IN OUTPUT. J3 IS ID NO.
00090
102
FORMAT (14,16A4.2A2)
00091
220
PRINT 200. N11.J3
0009?
200
FORMAT (1H1, SOLUTION TO MULTIPLE BUOYANT DISCHARGE",
00093
1» PROHLEM WITH AMBIENT CURRENTS AND VERTICAL GRADIENTS",//,
00094
? P0X.18A4," 1.0. = "•14.//)
00095
42
CONTINUE
00096
37
RFAD(5.lOSHJO.UA.TO.CO.TAO.CAO,TIZ,CIZ
00097
105
FORMAT (8F10.6)
00098
RFAll (5,100) DO, TH1,TH?,SF,SPACE,H
00099
PRINT 2000,UO,TO.CO.r>0.SPACE,H
00100
PRINT 200l,UA,TAO,CAO
00101
PRINT 2002,TIZ.C1Z
00102
2000
FORMAT (1M ,"DISCHARGE VELOCITY = "»F5.?."-M/S • TEMP. = ",F5.2
00103
x."-OEG C •« SALINITY = »,F5.2,"-PPT",/,1H ," DIAMETER = ",F5.2,
00104
K.I-M SPACING = "•F5.2• M oo DEPTH = » , F 7 . 2 , "-M")
00105
2001
FORMAT (1H ,"AMBIENT CONDITIONS AT DISCHARGE ELEVATION, VELOCITY
00106
X = "«F 7.2 »"-M/S «» TEMP. = ",F 7.2,"-DEG C •• SALINITY = ",F7.2,"
00107
XPPT")
00 I 08
2002
FORMAT (]H ."AMblENT STRATIFICATION GRADIENTS TEMPERATURE
00109
X •» ",F10.6,"-0EG C/M •• SALINITY •• ",F 1 0 .6,"-RPT/M"//>
00110
h = H/D0«2.
0011 1
CALL SIGMAT(PIO,TO,CO)
001 1?
CALL SIGMAT(PAO,TAO,CAO)
001 13
POVO= 0•5°D0/U0
00114
P=UA/UO
00115
BF=flE°TO/PIO
00116
GA=GA °CO/P10
00117
PI0=ABS(PA0-PI0)
00118
FR = SQRT (PII)/PI0«D0®9.8>
00119
FR=U0/FR
00120
TIZ=TIZ»DO/<2.«T0)
00121
CIZ=CIZ»D0/(2.®CO)
00122
TIO=TAO/TO
00123
CIO=CAO/CO
00124
SF = SF/DO«2.
00125
SPACE = SRACE/DO
00126
100
FORMAT (aFl0.6)
00127
9876
R1 = R
00128
RE = 1.0 E*6
00129
PTO = 1.0 ~(8E«(1.0-T10) *GA«(I.fl-CIO))
63
-------
U.O 1 JU
= -Ht«11/ - GA»CI/
oom
PRINT 207. TI0.CI0.PI0,TIZ,C1Z,PIZ
0 0 132
207
FORMAT (1H AMBIENT CONDITIONS AT DISCHARGE
00133
1" NONI) JMENSIONAL"./, " TEMP=", F H . 5 . 2X , "S AL INI T V = " , F8 . 5»2X,
00 1 34
? •'nENSTTY = ".F8.b«3X.HTFMPGR = ,,,F9.6.2X."SALGR = ",F9.6,2X,
00135
3 '.'DFNGR = ", F 1 0 . 7 ,/)
00136
40
PRINT 202. SPACE
00137
202
FORMAT ( 1H ,¦> PORT SPACING L/D = ",F6.2)
00138
IF (ABS(PIZ) - 1.0E-10) 90.90.91
00139
90
T = 1.0 E*8
00 140
GO TO 92
oom
91
T = (PIO-l.0)/(-P1Z)
001 <.2
92
PRINT 206. FR.R.T
00 1 o
206
FORMAT (1HU." FROUDE NO =•'.F9.2«4X."VELOC1TY RATIO =".
00144
1 F7.3.4*."STRATIFICATION NO = ",F7.1./)
00 145
FR = 2.®FR«FR
on 146
SPACE = 2.«SPACE
00147
nsn = OS
00148
TTT1 = THl
00149
TTT2 = TH2
0(1150
R.3 = R»SIN(0.01745329«TH1)«COS(0.01745329,,TH2)
00151
IF (R3-1.0) 60.61.60
0015?
61
93 = 0.0
00153
60
CALL ZFE
00154
C ZFF TAKES CARE OF COMPUTING IC AT END OF ZFE.
00155
OS = DS11
00156
THl = TTT1«0.01745329
00157
TH2 = TTT2«0.01745329
00158
TIE = T1E°0.01745329
00159
T2F = T2E°0.01745329
001 60
C CONVFRT TH FROM DEGREES TO RADIANS
00161
C TH
ANGLE AT S=0. TE = ANGLE AT SE
00162
SI = SIN(TIE)
00163
C? = COS(T2E)
00164
PS 12 = R'Sl^C?
00165
IF (RE - 1.0H»4) 20,21.21
00166
20
PRINT 214, RE
00167
214
FORMAT (1H FLOW MAY NOT HE FULLY TURBULENT, RE = ",F8.0,/I
00168
21
PRMT(1) = SET
00169
DRMT<2> = SF
00170
°RMT<3> = LIS
00171
PRMT(4)= 0.001
001 72
36
1 = -1
00173
N = 6 .
00174
C L=-
1 USED TO INITIALIZE OUTP. N = NUMBER OF DFFERENT1AL EONS.
00175
«= 0
00176
C K COUNTS NO OF INTEGRATIONS.
00177
no 2 I=l«6
001 78
2
FP(I) = 0.1666667
00179
STOMA = HUE
001 fiO
PRE? = BHd'RBE
001 ft)
F(l) = RRE2°CI I 1~(.5-1 I 1)»RS12)
001R2
F(2)=OTE°t3HE3&( I IP* fl Il-I
001*3
F<3> = l)CEortRE2» (112* ( 1 Il-I I2MRS12)
00184
F(4)=A8E2»II2«(1.-RSI 2)O02*2.<> (1.-RSI2)«RS1?°8HE2»111
00185
1 ~RS12°RS12°RHE2».5
00186
F(5) = TIE
00187
F (6) .= T2E
00 188
C SET INITIAL CONDITIONS AT END OF POTENTIAL CORE.
00 189
PRINT 208
00 1 90
208
FORMAT (1H0.20X,"Z0NF OF ESTAbLISHEO FLOW',/
00191
1 1H • 5x , "S", 8X,"X",7X,'
-------
uoiso
Ih 1 IMLt- -111 1 Ut 1 1 , I 1
00197
1 1
CALL EXTR {S.X,Y,Z.B.PU,DT,5.D0)
00198
10
PR INT 212, K,IhLF.PMMT(4),P«MT(5),PRMT(3)
00199
212
FORMAT (1H NO OF INTEGRATION STEPS=»,14.3X."NO OF HALVINGS"
00200
1 I3.3X,'>ABSEHR = ",F9.5,3X."0PMT (5) =",F9.3." FINAL SPACE ="
30201
2 .F9.3)
00202
GO TO 1
J0203
221
STOP
00204
FND
0020S
SUBROUTINE SIGMAtIrho.T.S)
90206
SIGO=-4.fl2K-4>oS».Rl49>®S-.O93
'10207
C = 1 .E-6«T«( (.01667«T-.8164)<»TM8.03)
00208
n=.001»T»(<.0Ol0843«T-.n9818)»r+4.7867)
00209
SUMT=(T-3.9B)«®2«(T»?83.)/(50 3.57"(T*67.26))
0021 0
RHO =(SItiO*.1324)•(1.-D*C»(SI GO-.1324))-SUMT
30211
DHO =RHO/1000.*1.
00212
RFTIJRN
00213
FND
00214
SUBROUTINE DERIV IS.F.FP)
00215
DIMENSION F ( 6) ,FP16) , Nil(18)
00216
COMMON /OUT/ L.TH1,TH2,H,Z,SOLO.OS.W.K,Nil,J3.J2,Jb.JJ,SPACE.SIGMA
00217
X,00,CD,A7,X XX
00218
S1 =S IN (F < 5) )
00219
S2=SIN(F(6) )
00220
C 1 =COS(F(5))
00221
C2=C0S(F(6))
00222
71 = 1* S2«(S - SOLD)
00223
CALL AMBIEN (Z1,CI.TI.PI.T1Z,C17,PIZ.VTI,VCI»VUI.VV1)
00224
C
AMRTFNT COMPUTES GRADIENTS OF T.C ) TURBULENCE OUANTTIES IN FR STREAM.
00225
CALL PRFLE
00226
C
PRFLF COMPUTES INTEGRAL OUANITITIES WITHIN THE JET
00227
C
F8
IS THE INTEGRAL OF DENSITY DIFFERENCE, AND 8 THE WIDTH OF TE JET
00228
FP<1) = E
00229
FP(2)=-F(1)°TIZ»S2 ~ VTI
00230
FP<3) =-F(1)«CIZ»S2 ~ VCI
00231
FP (4) = E*R°S1«C2 »FP«S2 ~ VUI
00232
F9 = F(4) - 0.25°E0E ~ VVI
00233
TF < BBS(C2) - 0.01) 1.1.2
00234
1
F (5) = F(6)
00235
FP(5) = 0.0
00236
GO TO 3
00237
2
FDY=0.
00238
UNX=-(C2»«2»C1»S1 I
00239
UNY=S2«»2»C2«02*C1°»2
00240
I)NZ = -(S2»C2°S1)
00241
UN=SORT
00247
FP(5)= (E«R°C1)/(C2°F9)*FDY
00248
3
FP(6) = (F8°C2 - (E»M*FD)«Sl«S2)/F9
00249
K =K ~ 1
00250
RFTURN
0 0 251
FND
00252
SUBROUTINE XINT(I,A.F1INT.F21NT.F3INT)
00253
DIMENSION XM25) , X 11 (2b) ,X 12(25)
00254
DATA (XA=1.0E—10,1.,1.1.1.2,1.3.1.4,1.5.1.6,I.7,1.6.1.9,
00255
X2.0,100.),(X I 1 = .4314,.4314..42^5,.4259,.4238,.4221,.4209,
00256
X.4201,.4196,.4194,.4193..4193), (X I 2=.436..3905..3576,
00257
X.3348,.3201,.31 14,.3068..304 7,.304..3038,.30377, .303?7)
0025H
1
TF ( A-XA.I I ) ) 2,3,4
00259
2
I = 1-1
00260
GO TO, 1
00261
3
F1 I NT = X 11 (I )
65
-------
.>0262
10 263
00264
D0 26S
00266
00267
0026B
00269
J0270
J0271
0027,?
:)02 7 3
(10274
00 275
0O27t>
00277
00278
00279
0 028 0
002B1
00282
002H3
00284
00285
00286
00287
0028H
00289
0029 0
00291
0029?
00293
00294
00295
00296
00297
00298
00299
00 300
0030]
00302
00303
00304
00305
00306
00307"
00308
00309
00310
00311
00312
0 0313
00314
00315
00316
0031 7
00318
00319
00320
003? 1
00322
00323
00324
00325
00326
00327
F 2 1 N I = A I 2 ( 1 )
r,0 TO 8
4 IF(A-XA(I*1)>5.6.7
6 FlINT=XI1{I*1)
F2INT = XI2(IM>
T = I ~ 1
r,0 TO 8
7 T = i * 1
RO TO 4
5 clINT = XIl(I)*(XIl(l*l)-XIl(I> )/(XA(1*1)-X A(I) )
-------
00328
00329
1)0330
00331
0033?
00333
00334
00335
00336
00337
II033H
00334
00340
00341
00342
00343
0034 4
00345
00346
00347
0034B
00349
00350
003S1
0035?
00353
00354
00355
00356
00357
00358
00359
00360
00361
0036?
00363
00364
00365
00366
00367
0036H
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
003HO
00381
0038?
00383
00384
00385
00386
00387
00388
00389
00390
00391
0039?
00393
f4=4.•1 1 J* A I
C5=U12»U12»AI3/PI
C6=2.»II4«AI2/P1
C7=2.«II3»U12«AI1/PI
r,0 TO 504
C oooo mfrgeo PLUME
503 CALL XINT(I.A.Al1,Al?,AI 3)
C1=.45«AI3/PI
r?=U12°Al3/Pl
C3=.315SB«AI3/P1
CW.#U12«»C1
C5=U12«>C2
C6=C 3
C7=0.5»C4
504 IF(AHSIU12)-.00]>112,112.113
112 RSUR=F(1)«F(1)«C3//Cl-C3»Z.BC2"F(1)/C1/C1-F(4)
CC=C3»F(1> »F(1)/Cl/Cl
RSQR=-BR-SURT
nSQR=BSQR/(2.«AA>
12 «=SORTIBSQP)
MGMA = B
DLUM=F(1)/(BSQR°C1)-C2/C1
IJCL = DLUM*U12
DTCL = F(2)/(HSQH®(01 UM«C6*C7))
DCCL = F (3)/(HSiiR» (Ol UM«C6*C7) )
PICL = RE°DTCL ~ GA»DCCL
IF I 11/ ((PI 0-1.)®FR>
IF(A.LT.2.)FH=2.»BS0P«IJ3»P1CL«AI1/IPI«FR»(PIO-1.))
IF(A.LE..95)FH=Fa»AI3/AI1/2.14068
C COMPUTE F8 - INTEGRAL OF DENSITY DIFFERENCE
FRLI = PICL®B/(FR«(PIO-1.0)»UCL°UCL)
GO TO 3
2 F8 = 0.0
FRLI = 0.0
3 IF < AB5(FRLI) - 1 0.0«ARS(1.O/FRL)) 10,13,13
13 FRLI = 1. O/FRL
10 CONTINUE
Al = .05
A?2 = 0.0
AA = Al ~ A22«FKLI«SIN(F (6) ).
CC = A7«W«SQHT <1.0-S12®S1?>
RH=R°ABS(UCL-R°S12)
28 TF (A.GE.2. )F = AA« (BB« (1 .-XXX/A) ~CCH)
IF(A.LT.2.)E=AA
RF AL 111.112,113,114
DIMENSION F(6).FP(to),PRMT(5),SO(4),SN(4).N11(18)
COMMON /OUT/ L,TH1,Th2,H,Z,S1,DS»R,K»N11,J3»J2,J5»JJ«SPACE~SIGMA
X,00,CO,A 7,XXX
COMMON /AMB/ T IZ , CIZ ,C10,T 10«P10,GA,AL»BEtFR,R1,FRL
COMMON /PLTE/ Jb,J7,J8.J9,J 10,X1M,Y1M,X2M,Y2M,Y3M
COMMON /ZF9/ SEr,XE«YE»ZE,BBE»DCE,DTE,LCNT»RUE,RTE»WCE, T1E»T2E
COMMON /STR95/ II.M.R3
COMMON /CONS/ 111,112,113,114,111
COMMON/WGM/MERGCK,ROVO
I = I I I
PII=3.14159
67
-------
00394
0039s
00396
00397
0039H
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
004 1 0
U04 1 1
00412
0041 3
00414
00415
004 16
0041 7
0041ft
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
oo43a
00439
00440
0044 1
00442
00443
00444
00445
00446
00447
00448
00449
004SO
0 04b 1
0 0452
00453
00454
00455
00456
00457
00458
00459
IF (L) 6.8.7
6 L= 0
L 19 = 0
Ml =-1
M2 = 0
71 = -OS
SNEW = PWMI(1)
SN (1) = SIN(TIF)
SN<2) = 51N ( T 2 E)
SN(3) = COS(Tit)
SN(4) = COS(T2E)
TIME = SET
x = XE
Y = YE
7. = ZE
urc = l.o
S 1 = PRMT <1>
T 3 I = T 10
C31 = CIO
IF(CIO-1.O.LT..000001)C3I=0.
IF (TIO- 1.0) 63.64.63
64 T 31 = 0.0
63 GO TO 10
8 . I. = l
SNEW = OS ~ SI
SO(1) = SN(1)
SO(2) = SN(2)
SO(3) = SN(3)
SO(4) = SN(4)
7 TF (M2) 18.18.19
19 ns = S-SNEw
TF (S-Sl - 1.0) 83.1,1
18 IF (S - SNEW) 20.1.1
20 IF (46S(S-SNEW) - O.OOl'SNEW) 1.1.83
1 SN(1)=S[N(F(5))
SN(?)=SIN(F(6))
SN(3)=COS(F(5))
SN(4)=COS(F(6))
X= X ~0.5 «S0 < 4 ) ) • ( S-S 1)
Y = Y ~0.5 ~ S0(1)«S0U)->
71 = Z
7 = Z~0•5° (SN ( 2) ~ S0(?) IMS-51)
no 9 1=1.4
9 SO(I)=SN(I)
10 PIZ = -f)E°T 17. - GA°C I Z
TI = TIO ~ TIZ«Z
CI = CIO ~ CIZ«Z
PI = Plo ~ PIZ*Z
IF (F ( 4 ) ) 21.21.55
55 Rl4=R«SN(1)»SN(4)
&=SPACF/SIGMA
IF(A-2.I 501.500.500
501 IF(A-.95)503.503.502
C SINGLE PLUME •«»»
¦500 r1=1 I 1
C?=«14/2.
C3=I 12
C4 = 2.14»11 1
r5=R14»R14/2.
C6=I I 2
C7 = II 1°R14
0,0 TO 504
C MERGING PLUMES 00°®
502 CONTINUE
519 FORMAT(/20X."ZONE OF MERGING FLOW'/).
IF(MEHGCK.EO.O)PRINT 519
68
-------
00460
MFRGCK=1
00461
CALL X1NT(I.A.A 11•A IA I 3)
00462
rl=2.»II3«Al 1/PII
00463
C2=R14«A13/PI 1
00464
C3 = 2.MI4«AI2/PII
00465
C4=4.®II3«AI1«R14/PII
00466
r5=Rl4°Rl4«Al3/PII
00467
C6=2.»II4»AI2/PII
00466
r7=2.»II3»Rl4»Al1/PIT
00469
r,0 TO 504
00470
C
o«o« MERGED PLUMES «•««
0047 1
503
CALL XINT(1,A»AI1,AI?,AI3)
00472
ri=.4S«Al3/PII
00473
C2=W14°A13/P11
00474
r3=.3155H»AI3/PII
00475
C4=2.»R14»C1
00476
C5=R14«C2
00477
C6=C3
00478
07=0.5«C4
00479
504
IF(A0S(RU)-.OO1) 12.12.99
00480
12
BSOR=F<1> »F < 1 >®C3/(F(4)®C1°C1>
00481
GO TO 15
00482
99
AA=C3»C2»C2/C1/C1-C4»C2/C1*C5
00483
RR=C4«F( I)/C1-C3«2.«C2»F(1>/C1/C1-F<4>
00484
CC=C3»F(1)«F(1)/C1/C1
00485
BSOR = -BB-SURT(H8°RB-4.«AA°CC)
U0486
RSOR = BSQW/(2.°AA)
00487
15
B = SQKT(BSQR)
00488
SIGMA = B
00489
S1GMA0 =SI(iMA»DO
00490
DLUM = F(1)/A»OELC» I 1 .O-CIO)
00502
c
DELU= (UCL - UI*51°C2)/UJO - UI> DELT = ( TCL - T 1) ~ ( TO - TI)
00503
F55 = F(5)»57.29578
00504
F 66 = F(6)»57.29578
00505
c
CONVFRT ANGLES BACK TO DFGRFES
00506
DMA = 0.5«(UCL ~ UCCl
00507
TIME = TIME ~ (S-Sll/UMA
00508
TF (FR -0.98E*6) 110.111,111
00509
110
IF (ABS((PI-PCL)°0> -1.0E-6°ABS(UCL»«2#FR«(PIO-1.0)1) 112,113,
00510
113
FRL = FR»(PIO-l.b>»UCL»UCL/
00511
fiO TO 112
00512
111
FRL = FR
00513
112
IF (L19) 70.70,71
00514
70
SOO = S/2.°D0
00515
XO = X/2.»D0
00516
YO = Y/2.»n0
00517
ZO = //2.»00
00518
P0=2.»F(1)
00519
TTMER=TIME»WOVO
00520
PRINT 104,SOO,XO,YO,70.F55.F66,SIGMAO,DELU,DELT,DELC»
00521
•PI,TIMER,HO
00522
104
FORMAT (1H , F9.2.1X.5F8.2, 1X,F8.2,3F9.3,1X,F9.5, 2F9.2)
00523
Ml = Ml ~ 1
00524
M = M« 1
00525
1.19 = L19 ~ 1
69
-------
00526
go iu i d
00S27
71
1.19 = 0
00528
72
PICL = PI-PCL
005?')
PICLO = PIO-l.O
00530
IF C M 21 40*40*4 2
00531
40
IF (SIGN(l.O.PICL) - SIGN<1.O.PICLO)) 41,42.41
00532
41
ns = 0.0
00533
Ml = 0
00534 .
M2 = 1
00535
MFXT=1
0 05 36
PRINT 213
00537
213
FORMAT (1H . 10X."PLUME HAS REACHED ITS EQUILIBRIUM HEIGHT
00538
1 •• STRATIFIED ENVIRONMENT")
00539
42
TF (Ml - 10) 36.32.3?
00540-
32
Ml = 0
00541
OS = 2.0°DS
00542
36
TF (Z) 3.3.31
00543
31
IF (M2) 90,90.91
00544
91
TF (Z-Zl) 21.90.90
00545
21
PRMTI5) = 1.0
00546
TH - 11
00547
PRINT 212
00548
212
format <1H .9X." PLUME HAS REACHED ITS MAXIMUM HEIGHT - " ,
00549
1 "STRATIFIED ENVIRONMENT"./)
00550
IF 83.83.81
00551
90
IF (7.-H) 3.3,4
00552
4
PRMTI5) =1.0
00553
C
thf
CARDS BELOW ARE FOR. USE WHEN THE PROGRAM IS USED FOR
00554
C
DISCHARGE into WATER OF DEPTH H. This MUST HE READ INTO PROGRAM
00555
PRINT 214
00556
214
FORMAT (1H ."PLUME HAS REACHED.WATER SURFACE"./)
00557
C
IF
WF HAVE REACHED THE WATER SURFACE - STOP.
00558
TF (J 6 > 83.83.81
00559
3
SNEW = SNEW ~ DS
00560
SI = s
00561
UCC = UCL
00562
TF
-------
OOSId
00593
00594
00595
00596
00S97
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
Sl(l) = bl t 1 *1>
Xl(I) = X1 (I~1)
rlili = v l (I~l)
71(1) = 21 (1*1)
R1(I) = 81(1*1)
mil (I) = ou 1(1*1)
DT1I1I = 011(1*1)
5 CONTINUE
M=3
1 SI = Z
B1 (M> = B
nu1(m > = ou
DTI(M) = OT
MR = M
RFTURN
4 IF (MR-3) 202.6.202
6 MR2 = MR-2
MR1=MR-1
7MAX = Zl(MRl) * DU1(MR1)«/(2.0®C))°D0
207 PRINT 102.ZMAX. ZM
102 FORMAT 1 1 HO." MAXIMUM PLUME HEIGHT. EXTRAPOLATED
1 "LINEAR ="»F8.2«" PARABOLIC =".F8.2./>
RFTURN
202 PRINT 203
203 FORMAT <1H .ERROR IN HPCG. IHLF = 11''./)
RFTURN
FND
SUHHOUTINE ZFE
THIS PROGRAM COMPUTES ThE SOLUTION TO BUOYANT JET PROBLEMS WITHIN
THF ZONE OF FLOW ESTABLISHMENT. I.E. UP TO XE.
THF SOLUTION CALCULATES RU.HT.RC«B.TH1.TH2 AS FUNCTIONS OF S.
DIMENSION PRMT(5).AUX(16.6).F(6).FP(6).N11(18)
COMMON /AMB/ TIZ.CIZ.CI0.TI0.PI0.GA.AL.BE.FR.R1.FRL
COMMON /OUT/ L.TH1.TH2*H,Z,SNfw,DS»H«K»Nl1.J3»J2»J5.JJ.SPACE»SIGMA
X.DO.CD.A7.XXX
COMMON /ZFsi/ SET.XE.YE.ZE.0HE»r>CE.DTE.LCNT,RUE»RTE.HCE,TlE.T2E
EXTERNAL DERIV1»OUTP1
PIZ = -BE°TIZ - GA'CIZ
TH1 = 0.01745329«TH1
TH2 = 0.01745329«TH2
PRMT(1) = 0.0
PRMT(2) = 100.0
PRMT(3) = 0.01
PRMT(4) = 0.0001
TF (ABS(COS(TH2)) - 0.01)1.1.2
1 TH1 = TH2
2 N=6-
L =0
LCNT = 0
71
-------
no 41=1,6
00659
3
FP(I> = 0.16666667
00660
F(l) = 1.0
0066]
F(2t = 1.0
00662
F(3I = 1.0
00663
F 14) = 0.0
00664
F(5) = TH1
(10665
F (6) = IH2
00666
C Ffl) = HU. F12) = Hi. FI3> = HC • F<4) =H, F<5> = Ihl. F<6) = TH2
U0667
rm_F = 0
00668
PRINT 210
U0669
210
POftMAT 11H0.20Xi"ZONF OF FLOw HSTARLISHMENT — ALL LENGTHS "»
00670
p'ARE IN MEIERS •',/
00671
1 1H .5X.''S",9X.''X",7X. »Y''.7X."Z".6X,"TH1''.5X."TH2',.9X
00672
1 » "R I,»SX."BU"»6)U"WT".6X."RC,,.5X »"OUCL" .4*•"DT CL"» 4X t "ijCCL" »
00673
2 6X."Pl"«7<«"TIME"./l
00674
CAUL 0UTP1 (PHM r(l)«F.FP«lHLf" » N« PKMT )
01)675
TALL HPCG (PHMT,F,FP,N,IHLFil)ERlV1.0UTPl.AUX>
00676
TF . FA =B/AO
006B4
TH1 = F 15)
(10685
TH2 = F(6)
00636
no 9 1=1.6
006H7
9
FPU) = 0.1666667
OO6R0
PRMT11) = SET
00689
PRMT12) = SET ~ 100.0
00690
PRMT 13) = 0.01
0069 1
PHMT(4) = 0.000 i
00692
1HLF = n
00693
("ALL 0UTP1 (PHMT 11) .F.FP.lHLF.N.PRMD
00694
CALL HPCG )
00718
R12 = R«S1*C2
00719
SI 2- = S1«C2
00720
nTO=1.0-TJO - TIZ«Z
00721
OCO=1.0-C10 - CIZaZ
00722
F = El° (0.0204«0.0144,1F (4) 1 o (1.0*E2iiS2/FRI o (A8SC1 ,0-H12> >t"3aH
00723
1 'SORT(1.0-S12«»2>)»(1-E4/SPACE)
72
-------
00724
C
E
IS THE ENTRAPMENT FUNCTION.
007?S
01=.45*.55°R12
00726
011=2.*<.12857*.37143*R12)
00727
02 = .31558*.13442*R1?
00728
n3=2.* t.06b76..06181»R]2)
00729
n4=.31558*.26fl04»R12. .4155««R12«R12
00730
05=2.° (.066 76*. 12362*R12*.3nS)62«R 12«H12)
00731
flf> a DTO»BE«MF«2)*F ~ .4S«F (2> »F U) >
00732
nsaOe* OCb®GA»(Fl3»»FI3)»O.S *0.i2»57«F(4)*F(4» ~»45*Ft3>•fUH
00733
0 = 0.5»F (1) «F (1) ~ F(])«F(4)«[)ft *F(4)«F <4)*0.5»i>5 -0.25®E«E
00 734
A A ] =0.5«F (1 ) «F<] ) ~nl«F<] »«FU) ~ D11'*Q.5*F I 4) u F (4)
00735
A (1.1 > = F (1 ) ~ 01«F(4)
00736
A(1.2>= 0.0
00737
A(1.3) = 0.0
00738
All.4) = D1°F < 1) ~ D1]«FUI
007 3")
C
FIRST EQUATION IS CONTINUITY.
00740
A(2•1) = 0.0
00741
AI2.2) = DT0»F(2) ~ n2°DT0«F(4)
00742
A<2.3> = 0.0
00743
A(2.4) = PT0*(D2»F(2> ~ D3*F14))
00744
c
SECOND EQU&TION IS ENERGY.
00745
A(3 ~ 1> = 0.0
00746
A (3.2) = 0.0
00747
A ( 31 3 > =" DCO* (F ( 3) ~ .D?«F(4»)
00748
A (3.4) = DC0»(D2*F(3> ~ D3°F < 4))
00749
c
THIRD EGUATON IS CONCENTRATION.
00/50
A(4 * 1) = F < 1) '~ D4*F(4)
00751
A ( 4.2) = 0.0
00752
A(4 * 3 > = 0.0
00753
A I 4.4) = 04»F (1 > ~ DS»F(4)
00754
c
FOURTH EQUATION IS MOMENTUM
00/55
FP(1) = E
00 756
TIZ7 = TIZ
00757
FP ( 2) = TIZZ*S2»1-AAI»0.5*F(2)«°? ~ F(2>«FI4)«D2 ~0.5*F(4)*»2»D3)
00758
FP(2)=FP(2)~VTE
00759
FP(3) = CIZ°S2*(-AAI ~0.5*F(3)««2 *T(3)«F(4)«D2 ~0.5*F(4)«*2*D3>
00760
FP(3) = FP< 3> +VCE .
00761
IF (L-2) 3.4.3
00762
4
A(2.2 > = O.S«03»F(41««2
00763
A<2.4) = 03»F(2)«F(4)
00764
A(3.3) = A(2.2)
00 765
A(3.4) = U3*F(3)*F(4)
00766
FP(2J=-TIZZ*S2*AAI*VTE«nTE
00767
FPO) = -CIZ«S2»AAI*VCE"DCE
00768
06 =¦
-------
00790 A?(2.2)=A(4.4)
00791 FP(2)=FP(4)
00792 CALL SIMQ(A2»FP.J«K)
00793 GO TO 7
00794 20 .1 = 3
00795 A3(1.1)=A(1.1)
00796 A3(l«2)=0.
00797 A3I1.3)=A(1.4)
00798 A3(2.1>=0.
00799 A3(2.2)=A(2.2)
00800 A3I2.3)= A(2.4)
OOflOl A3(3.1)=A(4,1)
1)0802 A3 (3.2) =0.
00*03 A3(3.3)=A(4.4)
00804 FP(3)=FP(4)
00805 CALL SIMQ(A3.FP.J.K)
00806 GO TO 7
00807 6 .1 = 3
00808 A3(1.1)=A(1.1)
00809 A3(1.2>=0.
00810 A3(1.3)=A(1,4)
00811 A3(2.1)=0.
00812 A3(2»2)=A(3«3)
00813 A 3 (2. 3) = A (3 . 4 )
00814 A3(3.1)=A(4,1)
00815 A3(3.2)=0.
00816 A3(3.3)=A(4.4)
00817 FP(2)=FP(3)
00818 FP(3)= FP(4)
00819 CALL SIMO(A3.FP.J.K)
00820 7 IF (K) 1.2.1
00821 2 TF(J-3)131.132.136
00822 131 rP(4)=FP(2)
00823 FP(3)=0.
00824 FP(2)=0.
00825 GO TO 136
00826 132 TF(DCO.NE.O.)GO TO 134
00827 FP(4)=FP(3)
00828 FP(3)=0.
00329 GO TO 136
00830 134 FP(4)= FP(3)
00831 FP(3)=FP(2)
00832 FP(2)=0.
00833 136 UNX = -(C2»»2»C1»S1)
00834 l.iNY = S2°»2*C2«»2,tCl«»?
00835 UNZ=- = E®R®C 1 /(0°C2)+FDY
00843 FP (6) = (F8°C2 - (E»R+Fn>°S1«S2>/Q
00844 RFTURN .
00845 1 WRITE(6.100)
00846 100 FORMAT (1HO."SINGULAR MATRIX-EXECUTION TERMINATED AT THIS POINT")
00847 STOP
00848 FNO
00849 SUBROUTINE 0UTP1 (S»F,FP.IH.N,PRmT)
00850 COMMON./AMB/ T 11. C IZ . C 10. T 11),PI0 .GA . AL . BE. FW. R 1. FWL
00851 COMMON /OUT/ L . IH1 , TH2.M. Z . SNEhl. OS. R .K , N1 1 . J3. J2. J5. JJ. SPACE. SIGMA
00852 X.00.CO.A 7.X XX
00853 COMMON /ZF9/ SET.XE»YE.ZE.HBE.DCE.OTE.LCNT.HUE»KTE.RCE.TIE.T2E
00854 COMMON /STR95/ II.M.R3
00855 COMMON/WGM/MERGCK•ROVO
74
-------
00857
0085B
00659
00660
00661
0066?
00863
00864
00865
00866
00867
00666
00669
00870
00871
00672
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00688
00889
00690
00691
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
DIMENSION f (6)lb).PHMI lb).bU14 I»bN(tI.Ml I IItt)
IF(L) 1.1.2
1 SNtl) = SIN(THl)
SN<2) = SIN < TH2)
SNI3) = COS(TH1)
SN(4> = COS(TH2)
TJ = PRMT(l)
r)S = PRMTI3)
BUN s 0.0
RTN=0.0
RCN = 0 • 0
RN=1.0
SSN=0.0
T1N=TH1»57.29576
T?N=TH2»57.29578
XN = 0.0
YN = 0.0
7N = 0.0
X =0 . 0
v = 0 . 0
7 = 0.0
JP = 100
IF 30.31 .30
31 L = 1
GO TO 6
30 L = 2
X = XE
Y = YE
7 = Z E
GO TO 6
2 IF (ABSIS-TJ >-0.01) 3.3.4
3 SN <1> = S1N(F < 5))
SN(2> = SIN(F(6))
SN(3) = COS IF(5))
SN(4) r COS(F(6))
X = X ~ 0 • 5* (SN ( 3) °SN < 4 ) ~ SO I 3>»S0(4))»DS
V = Y ~ 0.5MSN(1)*SN(4) ~ SOI] )«S0(4) )»r)S
7=2* 0.5*(SN(2) ~ SO(2))°DS
6 no 5 1=1.4
5 SOI I) = SN 11)
TIME = S
PTZ = -BE«TIZ - GA°CIZ
CI = CIO *CIZ»Z
TI = TIO*TIZ«Z
PI = PIO ~ PIZ»Z
F55 = F(5)°57.29576
F66 = F (6) <>57.29578
TF (JP-100) 50.51.51
51 DUU = (1.0-R»SN(1)«SN<4)I/I1.0-R3)
DELT = (l.O-Tll/ll.O-TIO)
DELC=0.
nFNOM=l.-CIO
IFIOENOM.NE.O.)DELC=(1,0-CIl/OENOM
C Ft?) = RT OK (TM-TI) /TO. F(3) SAME F"OP C
C DUU = IUO-U12)/(UO - U120 >. DELT = ITM-II)/ITO-T10)
IF (L-2) 22.23.22
23 I1ELT = F I 2)/ 11 • 0-T 10)
n€LC=0.
nENOM=i.o-cio
IF IDEN0M.NE.0.)DELC=FI 3)/DENOH
100 FORMAT (1H .F9.2»lX,SFfl.2.2X»2F8.3.16X.3F8.3«2X» F9.5.F9.1)
22 SOO = S/2.«00
xoo = x/2.«no
YOO = Y/2.°l)0
ZOO = Z/2.«D0
RRO = F(4)/2.°D0
75
-------
00922
RUOO = F( 1)/2."00
00923
riTOO = F 121/2.*00
00924
ncoo = F<3)/2.«DO
00925
TTMER=TJME«R0V0
00 926
TF (1-2)1001,1002,1001
00927
1002
PRTNT 10n.SOO,XOO,YOO.ZOO.F55«Fb6«HRO.RUOO,r>Ull.DELT,
00928
IDFLC.PNTIHfR
009?9
C,0 TO 1004
00930
1001
PRINT 100 31SOO,XOO. YOO, 700,F55,F66,BRO,RUOO.OT 00.OCOO,UUU
00931
lOELT.OELC.PI.TIME
00932
11)04
RPI_OT=F(4)»2.«S0RTIP.)
00933
1 003
FORMAT (1H .F9.2,lX.SFfl.2.2X.7K6.3,2X«F9.5,F9.1>
00934
C
WRITE(7.1000)S»Y,X,Z.BFLOT,F 11>«F<2),F(3>•OUU»OELT tOELC.
00935
C
»TIMEW
00936
C1000
FORMAT(5F8.3.6F5.3.FQ.3)
009 37
M = M ~ 1
00 93ft
TF MYO-YN)/I T 200. SOO
00984
2on
FORMAT (1H STARTING LENGTH, T =",F9.3,)
00985
(. = -2
00986
PRMTI5) = 1.0
0048 7
pfturn
76
-------
0098B
13
TJ = TJ ~ PRMT < 3)
009B9
RETURN
00990
40
TJ = TJ ~ PRMT(3)
00991
RETURN
00992
21
IF (Fill) 10.10.40
00993
10
PRMT < 5) = 1.0
00994
SFT = SSN - RUN"(SSN-SS>/(RUN-RUOI
00995
RUE = RUO ~ (SE T—SS1•(RUO-HUN)/(SS-SSN)
00996
ORE = BO ~ (SE T-SS)•(BO—RN)/ISS-SSN)
00997
me. = WTO ~ (SEr-SS)°(RTO-RTN>/(SS-SSN)
0099ft
nCE = RCO ~ (SE1-5SI«
00999
TIE = T10 ~ (SEr-SSIMTlO«TlN)/»lYO-YN)/(SS-SSN>
01003
7F = ZO ~ (SFT-SS) » (70-7N) / (i,5-SSN»
0 1 004
OUT = <1.0-R«SIN<0.0]745oTlE)»C0S(0.01745FNOM=l .O-CIO
01008
IFIDENOM.NE.O.)OCT =OCE/DENOM
01009
PI = PIO ~ PIZ'ZE
01010
SOO = SET/2.«DO
01011
PRINT 205.500
01012
205
FORMAT <1H STARTING LENGTH, VELOCITY = ". F9.31
01013
*00 = XE/2.°00
01014
YOO = YE/2.°nO
0101S
700 = ZF_/2.»DO
01016
HRO = H8E/2.*D0
01017
RUOO = RU£/?.*D0
01018
DTOO = DTE/2.
01013
OCOO = OCE/2.
01020
PRINT 1'00»SOO.XOO.YOO»ZOO«T1E»T2F»RBO.RUOU»DUTiDTT»
01021
lOCT.PI.TlME
01022
4
RFTURN
01023
FND
01 024
SUBROUTINE SIMQ(A.B.N.KS)
SIMO
01025
C
SIMQ
01026
c
SIMO
01027
0 I MENS I ON A (1) ,d( 1)
SIMQ
01028
c
SIMQ
01029
c
FORWARD solution
SIMO
01 030
c
SIMO
01031
TOL=0.0
SIMO
01032
K S = 0
SIMO
01033
JJ = -N
SIMQ
01034
no 65 J=1,N
SIMO
01035
>(Y = J*1
SIMO
01036
JJ=JJ*N.1
SIMQ
01037
fUGA = 0
SIMO
01038
IT = JJ-J
SIMQ
01039
DO 30 I=J.N
SIMQ
01040
c
SIMO
01041
c
SEARCH FOR MAXIMUM COEFFICIENT IN COLUMN
SIMQ
01042
c
SIMQ
01043
T J=IT» 1
SIMU
01044
IF(ABS(BIGA)-ABS(A(Ij)>) 20.30.30
SIMQ
01045
20
niGA=A(IJ)
SIMQ
01046
TMAX=I ¦
SIMQ
01047
30
CONTINUE
SIMQ
01048
c
SIMQ
01049
c
TEST FOR PIVOT Lt'SS THAN TOLERANCE (SINGULAR MATRIX*
SIMQ
01050
c
SIMQ
01051
IF
-------
01 os4
C
010SS
C
interchange hows tf necessary
01 056
C
010S7
40
T1=J*N«(J-2)
010SH
tt=imax-j
010 59
no SO K = J.N
01060
T1 = 11*N
01 061
12= 11~ I T
0106?
SAVE=A( I 1 )
01 063
A(I 1)=A(12)
01 064
A(12)= S A V E
01 06S
C
01066
c
DIVIDE EQUATION Hy LEADING COEFFICIENT
01 067
c
01068
so
A CI I)=A(I 1)/BIGA
01 069
<;ave = h ( I MA A )
01070
rt(I MAX)< J)
01071
R(J)=SAVE/BIGA
01072
c
01 073
c
ELIMINATE NEXT V4RIARLE
01 074
c
01 07S
TF(J-N) SS.70.5S
01076
5S
TQS=N«(J-l)
01077
no 65 IX=JY.N
01078
TXJ=IQS* IX
01 079
T T = J - I X
01080
no 60 JX=JY.N
01081
TXJX=N»)°a(IXJ)>
01085
c
01 086
c
BACK SOLUTION
01 087
c
01088
70
NY=N-1
01 089
T T = N 0 N
01090
no 80 J=1.NY
01091
I A =IT-J
01092
TR=N-J
01 093
TC = N
010 94
no 80 K=1.J
01095
R(IB)=fl(IBl-A(IA)»H(IC)
01096
T A=IA-N
01097
80
TC=IC-1
0109H
RFTURN
01 099
FNO
01 100
c
01101
01102
c
SURROUTINE HPCG(PRMT.Y.DEWY.NDIM.IHLF.KCT.OUTP,AUX)
01103
c
SUBROUTINE HPCG
0110 4
c
01 10S
c
PURPOSE
01106
c
TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY GENERAL
01107
c
DIFFERENTIAL EOUATIONS WITH GIVEN INITIAL VALUES.
0110S
c
01 109
c
USAGE
01 1 10
c
CALL HPCG (PRMT.Y.DERY.NDIM,IHLF.FCI.OUTR.AUX)
01111
c
PARAMETERS FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT.
01112
c
01113
c
DESCRIPTION OF PARAMETERS
01114
c
PRMT - AN INPUT AND OUTPUT VECTOR WITH DIMENSION GREATER
01115
c
OR EQiial TO 5. WHICH SPECIFIES 1 HE PARAMETERS OF
01116
c
THE INTERVAL AND OF ACCURACY AND WHICH SERVES FOR
01117
c
COMMUNICATION BETWEEN OUTPUT SUrtROUTINE (FURNISHED
01 llfi
c
HY THE USER) AND SURROUTINE HPCG. EXCEPT PRMT(S) 1
01119
c
THE COMPONENTS ARE NOT DESTROYED BY SUBROUTINE 1
MMU
SIMU
S I HQ
SIMQ
SIMO
SIMU
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMO
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMO
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMU
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
SIMU
SIMQ
SIMQ
SIMU
SIMQ
SIMQ
SIMQ
SIMQ
SIMQ
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
HPCG
78
-------
01
I lit)
C
rIPCG ANO 1 HEY ARf
HPCG
01
121
C
PRMT( 1 )-
LOWER BOUND OF THE INTERVAL (INPUT 1,
HPCG
01
122
C
PHMT(2) -
UPPER ROUND OF THE INTERVAL (INPUT),
HPCG
01
123
C
PRMT(3)-
INITIAL INCREMENT OF THE INDEPENDENT VARIABLE
HPCG
01
124
C
(INPUT).
HPCG
01
125
C
PHMT (4)-
UPPER ERROR BOUND (INPUT). IF ABSOLUTE ERROR IS
HPCG
0 1
126
C
GREAfFR Than PRMT(4). INCREMENT .GETS HALVED.
HPCG
01
127
C
if increment is less than prmtoi and absolute
HPCG
01
128
C
ERROR LESS THAN PRMT(41/50. INCREMENT GETS DOUBLED
.HPCG
01
129
C
THE USER MAY CHANGE PRMT(4) BY MEANS OF HIS
HPCG
01
130
C
OUTPUT SUBROUTINE.
HPCG
01
131
C
PRMT<5)-
NO INPUT PARAMETER. SUBROUTINE HPCG INITIALIZES
HPCG
01
132
C
PRMT(5)=0. IF The USER WANTS TO TERMINATE
HPCG
01
133
C
SUBROUTINE HPCG AT ANY OUTPUT POINT, HE HAS TO
HPCG
01
1 34
C
CHANGF P°MT(5) TO NON-ZERO BY MEANS OF SUBROUTINE
HPCG
01
135
C
OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE
HPCG
01
1 3*
C
FEASIBLE IF ITS DIMENSION IS DEFINED GREATER
HPCG
0 I
137
C
THAN 5. HOWEVER SUBROUTINE HPCG DOES NOT REQUIRE
HPCG
01
138
C
ANO CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL
HPCG
0]
139
c
FOR HANDING RESULT VALUES TO THE MAIN PROGRAM
HPCG
01
140
c
(CALLING HPCG) WHICH ARE OUT AINEO BY:SPECIAL
HPCG
01
141
c
MANIPULATIONS WITH OUTPUT OATA IN SUBROUTINE OUTP.
HPCG
01
142
c
Y
INPUT VECTOR OF INITIAL VALUES. (DESTROYED)
HPCG
01
143
c
LATEWON Y IS THE RESULTING VECTOR OF DEPENDENT
HPCG
01
144
c
VARIARLES COMPUTED AT INTERMEDIATE POINTS X.
HPCG
01
145
c
DEMY
INPUT VECTOR OF ERROR WEIGHTS. (OESTROYEO)
HPCG
01
14b
c
THE SUM OF ITS COMPONENTS MUST HE EQUAL TO 1.
HPCG
0 1
147
c
LATERON OERY IS THE VECTOR OF DERIVAT 1 l/ES. WHICH
HPCG
01
148
c
BELONG TO FUNCTION VALUES Y AT A POINT X.
HPCG
01
149
c
NDI M
AN INPUT VALUE» WHICH SPECIFIES THE NUMBER OF
HPCG
01
150
c
EQUATIONS IN THE SYSTEM.
HPCG
01
151
c
IhLF
AN OUTPUT VALUE. WHICH SPECIFIES THE NUMBER OF
HPCG
01
152
c
BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS
HPCG
01
153
c
GREATFR THAN 10. SURPOUTINE HPCG RETURNS WITH
HPCG
01
154
c
F.HROR MESSAGE I HLF = 1 1 INTO MAIN PROGRAM.
HPCG
01
155
c
ERROR MESSAGE IHLF=12 OR IHLF=13 APPEARS IN CASE
HPCG
0 1
156
c
PRMT(3)=0 OR IN CASE SIGN(PRMT(3)).NE.SIGN(PRMT(2)
-HPCG
01
157
c
PRMT(1)) RESPECTIVELY.
HPCG
01
158
c
FCT
THE NAME OF AN EXTERNAL SUBROUTINE USED. IT
HPCG
01
159
c
COMPUTES THE RIGHT HANI) SIDES OERY OF THE SYSTEM
HPCG
01
160
c
TO GIVEN VALUES OF X AND Y. ITS PARAMETER LIST
HPCG
01
161
c
MUST RE X.Y.OERY. THE SUBROUTINE SHOULD NOT
HPCG
01
162
c
DESTROY X AND Y.
HPCG
01
163
c
OUTP
THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED.
HPCG
01
164
c
its parameter list must be x.y.oery,ihlf.ndim.prmt
.HPCG
01
165
c
none OF THESE PARAMETERS (EXCEPT, IF NECESSARY,
HPCG
01
166
c
PRMT(4).PRMT(5)....) SHOULD BE CHANGED BY
HPCG
01
167
c
SUBROUTINE OUTP. IF PRMT(5) IS CHANGED TO NON-ZERO
.HPCG
01
16B
c
SUBROUTINE HPCG IS TERMINATEO.
HPCG
01
169
c
AUX
AN AUXILIARY STORAGE ARRAY WITH 16 ROWS AND NDIM
HPCG
01
170
c
COLUMNS.
HPCG
01
171
c
Hf?CG
01
172
c
REMARKS
HPCG
01
173
c
THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM. IF
HPCG
0 1
174
c
(1) MORE
THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE
HPCG
01
175
c
NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE
HPCG
m
176
c
IHLF
= 11) .
HPCG
01
177
c
(21 INITIAL INCREMENT IS EQUAL TO 0 OR HAS WRONG SIGN
HPCG
01
178
c
(ERROR MESSAGES THLF=12 OR IHLF=13),
HPCG
01
179
c
(31 THE
WHOLE INTEGRATION INTERVAL IS WORKED THROUGH,
HPCG
01
180
c
(4) SUBROUTINE OUTP HAS CHANGED PRMT(5) TO NON-ZERO.
HPCG
0]
181
c
HPCG
01
182
c
SUBROUTINES
AND FUNCTION SUHPROGRAMS REQUIRED
HPCG
01
183
c
THE EXTERNAL SUBROUTINES FCT(X,Y,DERY» ANO
HPCG
01
184
c
OUTP(X.Y
.DERY.IHLF.NDIM.PRMT) MUST HE FURNISHED BY THE USER
.HPCG
01
185
c
HPCG
79
-------
U] 18b
C
METHOD
HPCG
01 107
C
EVALUATION IS DONE RY MEANS OF H AMM IN6S MODIFIED PREDICTOR-
HPCG
01181
C
CORRECTOR METHOD. IT IS A FOURTH ORDER METHOD* USING 4
HPCG
01 1«9
C
PRECEEDI N(> POINTS FOR COMPUTATION OF A NEW VECTOR Y OF THE
HPCG
01 190
C
DEPENDENT VARIABLES.
HPCG
011^1
C
FOURTH ORDER RuNGE-KUTTA METHOD SUGGESTED BY RALSTON IS
HPCG
Ul 192
C
USED FOR ADJUSTMENT OF THE INITIAL INCREMENT AND FOR
HPCG
01 193
C
COMPUTATION OF STARTING VALUES.
HPCG
01 194
C
SUHROUTINE HPCG AUTOMATICALLY.ADJUSTS THE INCREMENT DURING
HPCG
01 19b
C
THE WHOLE COMPUTATION BY HALVING OR DOUBLING.
HPCG
01 196
C
TO GET FULL FLEXIBILITY IN OUTPUT. AN OUTPUT SUBROUTINE
HPCG
01 197
C
MUST HE CODEO RY THE USER.
HPCG
01 19B
C
FOR REFERENCE. SEF
HPCG
01 199
C
(1) RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL
HPCG
01 200
C
COMPUTERS, WILEY. NFW YORK/LONDON, I960. PP.95-109.
HPCG
01201
C
(2) RALSTON, OUNGE— .Y < 1 )»f)ERY(1>,AUX(16.1)
HPCG
01209
c
HPCG
01210
c
HPCG
01211
N=1
HPCG
01212
THI.F = 0
HPCG
01213
X=PRMT(1)
HPCG
01214
h=prm r13)
HPCG
01215
DPMT(5)=0.
HPCG
0121b
no 1 I=1.NUIM.
HPCG
0 1217
AUX(16.I)=0.
HPCG
01 218
All* (15,1) =OERY ( I )
HPCG
01219
1
AUX(1«I)=Y(I)
HPCG
0 1220
IF
-------
01252
N=1
HPCG
0 1 253
T5W = 2
HPCG
01254
r,oro loo
HPCG
01255
C
HPCG
01256
1.3
X = X»H
HPCG
01257
CALL FCT (X.Y.DEKY)
HPCG
01 2S8
N = 2
HPCG
01259
no 14 I = 1t NO IM
HPCG
0 J 260
MIX (2. 1 ) =Y ( I J
HPCG
01 261
14
MJX(9«I)=DE»Y(1)
HPCG
"1262
TSW = 3
HPCG
01263
noro loo
HPCG
01264
C
HPCG
01265
C
COMPUTATION OF TEST VALUE DELT
HPCG
01266
15
OFL T = 0.
HPCG
<> 1 26?
no 16 i=i.NnjM
HPCG
01268
16
DFL T =DELT~ AUX(15.1)«ABS1Y(I)-AUX(4. I))
HPCG
01 269
DFL T = •06666667°DEL T
HPCG
01270
TF«H»(.375»AUX (8.I)~.791666 7®AUX (9.I)
HPCG
0 1 295
I-.208 3333"AIJX(10*1)~.04166667°0EPY(I))
HPCG
01296
23
X = X.H
HPCG
01297
N=N*l
HPCG
01298
CALL FCT(X.Y.DERY)
HPCG
01299
CALL OUTP(X.Y.DERY.IHLF.NDIM.PRMT)
HPCG
01300
T F1PRMT(b) J6.24.6
HPCG
01301
24
TFtN-4)25.200.200
HPCG
01302
25
no 26 1=1.NOIM
HPCG
01303
AUX(N.I)=Y(I)
HPCG
0 1 304
26
4UX(N*7.I)=DERY(I)
HPCG
01 305
IF(N-3)27.29.20U
HPCG
01306
C
HPCG
01307
27
no 28 1=1.NDIM
HPCG
01308
nFLT = AUX(9.1)~ AUX(9.1)
HPCG
01309
dflt=oelt*uelt
HPCG
01310
28
Y(l)= AUX(1.I>~.J333333»H»(AUX(8.I>~DELT*AUX(10.I))
HPCG
01311
ROTO 23
HPCG
01312
C
HPCG
01313
29
no 30 1=1.NDIM
HPCG
01314
nELT = AUH(9.1)~AUX(10.I).
HPCG
01315
nFLT = DELT*DELT ~L)ELT
HPCG
0131b
30
Y(I)=AUX(1.I)».375«H0(AUX(8.I)»OELT.AUX(1 1 . I))
HPCG
01317
ROTO 23
HPCG
81
-------
0
31 ft
C
HPCG
0
314
C
THp FOLLOWING PART OF SUHROUTINF HPCG COMMUTES HY MEANS OF
HPC6
0
3?0
C
RUNGE-KUTTA METHOD STARTING VALUES FOR THE" NOT SELK-START 1NG
HPCG
0
321
C
PREDICTOR-CORRECTOR METHOD.
hpcg
0
322
100
no 101 I = l.Nf)lM .
hpcg
0
323
7=H«AUX =Z
hpcg
u
325
101
Y(I)= AUXIN.T)~.4°Z
hpcg
0
326
C
7 IS AN AUXILIARY STORAGE LOCATION
hpcg
n
327
c
hpcg
0
328
¦7 = X«.4» H
hpcg
u
329
CALL FCT (Z . Y.DEHY)
hpcg
0
330
no 102 i = i•no Ih
hpcg
0
331
7=H«OERY(l>
hpcg
0
332
AUX (6*1) -Z
hpcg
0
333
102
Y ( 1 ) =AUX (N. I) ~ .2969776«AUX (5, 1 ) ~ . I587596°Z
HPCG
0
334
c
HPCG
a
335
7=x ~ .4557372®H
HPCG
a
336
CALL FCT(Z.Y.OERY)
HPCG
0
337
no 103 I = 1»NDIM
HPCG
0
33fl
7=H»DERY(I)
HPCG
0
339
AUX(7 * I ) =Z
HPCG
0
340
103
Y(I)= AUX(N.I)*.2181004«AUX(5,I)-3.050965"AUX(6.1)*3.a32Hb5»Z
HPCG
a
34 1
c
HPCG
0
34 2
7 = X »H
HPCG
i)
34 3
CALL FCT(Z.Y.DERY)
HPCG
a
344
no 104 I=l,NOIM
HPCG
0
345
1040Y ( I ) = AUX (N. I ) ~. 1747603°AUX IS. I ) -.S51 4807°AUX 16. I )
HPCG
0
34 6
1 ~1 .205536®AUX< 7.1)~.171184fl»H«nKRY(I>
HPCG
0
347
r,010(9. 13. 15.21 ) .ISW
HPCG
0
34R
c
HPCG
0
349
c
POSSIBLE HREAK-POINT FOR LINKAGE
HPCG
0
350
c
HPCG
0
351
c
STARTING VALUES ARE COMPUTED.
HPCG
0
35?
c
NOW START HAMMINGS MODIFIED P«EDICTOR-COHRfCTOR METHOD.
HPCG
0
353
200
TSTEP=3
HPCG
0
354
201
T F (N-8) 204. 202.204
HPCG
0
355
c
HPCG
0
356
c
n=« causes the mows of aux to change their storage locations
HPCG
0
357
20?
no 203 N=2.7
HPCG
6
35B
no 203 1=1.ND1M
HPCG
0
359
AUX (N-l . I ) =AUX =AUX(N*7i I)
HPCG
0
361
N = 7
HPCG
0
362
c
HPCG
0
363
c
N LESS THAN 8 CAUSES N*1 TO (SET N
HPCG
0
364
204
N=N*1
HPCG
0
365
c
HPCG
0
366
c
COMPUTATION OF NEXT VECTOR Y
HPCG
0
36 7
DO 205 [=1.NDIM
HPCG
0
368
AUX(N-l.I)= Y(I)
HPCG
0
369
205
aux in»6. i i =r>ERy < I >
HPCG
0
370
X = X + H
HPCG
0
371
206
TSTEP=ISTEP+1
HPCG
0
372
DO 20 7 I=1.NDIM
HPCG
0
373
OnFL T = AUX(N-4.I) ~1.3J3333*H6(AUXIN*6.1)*AIH IN*6.I)-AUXIN*5.1)~
HPCG
0
374
1AUX(N*4,It*AUX(N*4.I 1 )
HPCG
0
3 75
Y(I)=DELT-.925619fl»AUX(16,If
HPCG
0
376
207
AUX(16.1)=OELT
HPCG
0
377
c
PREDICTOR IS NOni GENERATED IN WOW 16 OK AUX. MODIFIED PREDICTOR
HPCG
0
37M
c
IS GENERATED IN Y. DFLT MEANS AN AUXILIARY STORAGE.
HPCG
0
379
c
HPCG
0
3*0
CALL FCT(X.Y.DERY)
HPCG
0
3fll
c
DERIVATIVE OF MODIFIED PREDICTOR IS GENERATED IN DERY
HPCG
0
3fl2
c
HPCG
0
3H3
110 208 1 = 1. ND IM
HPCG
82
-------
U1 384 ODFL I = . 12b«(9.®AUX (N-l,l)-AUX*J.«H®(UEkY* HPCG
01305 1AUX(N*6.I)-AUX HPCG
01388 C HPCG
01389 C TEST WHETHEP H MUST RE HALVED OR DOUBLED HPCG
01390 OFLT=0. HPCG
01391 no 209 1 = 1.NDIM HPCG
01392 209 DFLT=DELT*AUX(15.1)®ARS(AUX(16,I>) HPCG
01393 TF(0£LT-PRMT<4)1210.222.222 HPCG
01 39<» C HPCG
01395 C H MUST NOT BE HALVED. THAT MEANS Yd) AWE GOOD. HPCG
0J396 210 CALL FCT(X.Y.DEHY) HPCG
01397 CALL OUTP(X.Y.DERY.IHLF.NDIM.PRMT) HPCG
01398 TF(PRMT(5)>212.211.212 HPCG
01399 211 IF(IHLF-11)213.212.212 HPCG
01400 212 RETURN HPCG
01401 213 TFIH®(X-PRMT(2)))214.212.212 HPCG
01402 214 IF(ABS(X-PRMT(2))-.1«ARS(H>)212,215.215 HPCG
01403 215 IF(DELT-.02»PHMT(4))216.216.201 HPCG
01404 C HPCG
01405 C HPCG
01406 C H COULD HE DOUBLED IF ALL NECESSARY PRECEEDING VALUES ARE HPCG
01407 C AVAILABLE HPCG
01408 216 TF(IHLF)201»201»217 HPCG
01409 217 fF=AUX(N-2, I > HPCG
01418 AUX(N —2 ,1) = AUX(N—4. I ) HPCG
01419 Al)X(N-3.I)=AUX(N-6.I) HPCG
01420 AUX(N*6.I)=AUX(N*5.I) HPCG
01421 AUX(N*5.I)=AUX(N*3,1 ) HPCG
01422 AtJX (N* 4 • I ) = AUX (N* 1 . I ) HPCG
01423 nFLT=AUX(N*b.I)*AUX(N*5.I) HPCG
01424 DFLT=OELT»DELT*OELT HPCG
01425 22 inAUX(16. I)=8.962963*(Y(1)-AUX(N-3.I)>-3.361111®H®(DERY(I) ~DELT HPCG
01426 1*AUX(N*4.1)> HPCG
01427 ROTO'201 HPCG
01428 C HPCG
01429 C HPCG
01430 C H MUST 8E HALVED HPCG
01431 222 IHLF = IHLF~ 1 HPCG
01432 IF(IHLF-10)223.223.210 HPCG
01433 223 H=.5«H HPCG
01434 TSTEP=0 HPCG
01435 DO 224 I = 1.NO IM HPCG
01436 OY(I)=.00390625®(80.®AUX(N-l.I)~135.®AUX(N-2.I)*40.°AUX(N-3. I) ~ HPCG
01437 1AUX(N-4.I)>-. 1 171875®(AUX(N*b.I)-6.«AUX(N*S.I)-AUX(N+4,I))»H HPCG
01438 r>AUX (N—4. I) =.00390625® (12.®AUX (N-l. 1 ) *135.MUX (N-2. I) ~ , HPCG
01439 1108.®AUX(N-3.I)~AUX(N-4,I) >-.0234375®(AUX ~18.®AUX(N*5.I>- HPCG
01440 29. ®AUX(N* 4.I))"H HPCG
01441 AUX(N-3.I)=AUX(N—2.1) HPCG
01442 224 AUX(N»4.I)=AUX(N.b.I) HPCG
01443 X=X-H HPCG
01444 OELT=X-(H*H) HPCG
01445 CALL FCT(DELT.Y.DERY) HPCG
01446 HO 225 I = 1»NDIM HPCG
01447 AUX(N—2.I)=Y(I) HPCG
01448 AUX(N*5.I)=DEWY(I> HPCG
01449 225 Y(I>=AUX(N-4,I) HPCG
83
-------
014b0 DF-'L r =UtL I - (rl»H I nrtu
01451 CALL FC T (DEL T • Y »DEH Y) HPCG
0145? 00 226 I = 1»NO IM HPCG
01453 DFLT=AUX HPCG
01454 r>FLT = OELT*DELT*OtLT HPCG
01 455 niUJX ( 16. I > = 0.962*63° UUX (N- 1 , I) -Y ( I) ) -3. 361 111«H» ( AUX HPCG
01457 226 AIJX (N» j, 1 ) =DEKY ( I ) HPCG
01458 r,oT0 206 HPCG
0)45
-------
REFERENCES
1. U.S. Environmental Protection Agency. April 25, 1978. Modification of
Secondary Treatment Requirements for Discharges into Marine Waters.
Federal Register, Vol. 43, No. 80, p. 17484-17508. Washington, D.C.
2. 33 USC 1251, et sec). PL 92-500 as amended December 28, 1977.
3. California State Water Resources Control Board. 1978. Water Quality
Control Plan for Ocean Waters of California. Sacramento, California.
4. EPA. 1976. Quality Criteria for Water. U.S. Environmental Protection
Agency, Washington, D.C. EPA-440/9-76-023.
5. Baumgartner, D. J. and D. S. Trent. 1970. Ocean Outfall Design: Part
I. Literature Review and Theoretical Development. U.S. Dept. of Int.,
Federal Water Quality Adm. (NTIS No. PB-203 749).
6. Briggs, G. A. 1969. Plume Rise. U.S. Atomic Energy Commission, Oak
Ridge Tennessee. (NTIS No. TID-25075).
7. Davis, L. R. and M. A. Shirazi. 1978. A Review of Thermal Plume Model-
ing. Keynote address. In: Proceedings of the Sixth International Heat
Transfer Conf., ASME, Aug. 6-11, 1978, Toronto, Canada.
8. Rouse, H., C. S. Yih and H. W. Humphreys. 1952. Gravitational Convec-
tion from a Boundary Source. Tellus, Vol. 4, pp. 201-210.
9. Priestley, C. H:. B. and F. K. Ball. 1955. Continuous Convection from an
Isolated Source of Heat. Quarterly Journal of the Royal Meteor. Soc.,
Vol., 81,"p. 144-157.
85
-------
10. Morton, B. R., G. I. Taylor and J. S. Turner. 1956. Turbulent Gravita-
tional Convection from Maintained and Instantaneous Sources. Proc. of
the Royal Soc. of London, Vol. A234, pp. 1-23.
11. Morton, B. R. 1959. Forced Plumes. Journal of Fluid Mechanics, Vol. 5,
pp. 151-163.
12. Abraham, G. 1963. Jet Diffusion in Stagnant Ambient Fluid. Delft Hy-
draulics Publication No. 29. Delft, Netherlands.
13. Fan, L. H. 1967. Turbulent Buoyant Jets into Stratified and Flowing
Ambient Fluids. Cal. Inst. Tech., Keck Hydraulics Lab., Rep. No. KH-R-15.
Pasadena, California.
14. Abraham, G. 1971. The Flow of Round Buoyant Jets Issuing Vertically
into Ambient Fluid Flowing in a Horizontal Direction. Delft Hydraulics
Publication No. 81, Delft, Netherlands.
15. Cederwall, K. 1971. Buoyant Slot Jets into Stagnant or Flowing Environ-
ment. Cal, Inst, of Techn., Keck Lab. Rep. No. KH-R-25. Pasadena,
California.
16. Roberts, P. J. W. 1977. Dispersion of Buoyant Waste Water Discharged
from Outfall Diffusers of Finite Length. Cal. Inst, of Techn., Keck Lab.
Report No. KH-R-35. Pasadena, California.
17. Sotil, C. A. 1971. Computer Program for Slot Buoyant Jets into Strati-
fied Ambient Environments. Cal. Inst, of Techn., Keck Lab. Tech. Memo
71-2. Pasadena, California.
18. Koh, R. C. and L. N. Fan. 1970. Mathematical Models for the Prediction
of Temperature Distribution Resulting from the Discharge of Heated Water
in Large Bodies of Water. U.S. Environmental Protection Agency, Water
Poll. Cont. Res. Series Rep. 1613 0DW0/70.
-------
19. Baumgartner, D. J., D. S. Trent, and K. V. Byram. 1971. User's Guide
and Documentation for Outfall Plume Model. U.S. Environmental Protection
Agency, Pac. N.W. Water Lab. Working Paper No. 80. Corvallis, Oregon.
(NTIS No. PB 204-577/BA).
20.; Winiarski, L. D. and W. E. Frick. 1976. Cooling Tower Plume Model.
U.S. Environmental Protection Agency, Corvallis Environ. Res. Lab.
EPA-600/3-76-100. Corvallis, Oregon.
21\ Davis, L. R. 1975. Analysis of Multiple Cell Mechanical Draft Cooling
Towers. U.S. Environmental Protection Agency, Corvallis Environ. Res.
Lab. EPA-660/3-75-039. Corvallis, Oregon.
22. Brooks, N. H. 1973. Dispersion in Hydrologic and Coastal Environments.
U.S. Environmental Protection Agency, Corvallis Environ. Res. Lab.
EPA-660/3-73-010. Corvallis, Oregon.
23: Winiarski, L. D. and W. E. Frick. 1978. Methods of improving Plume
Models. Presentation at Cooling Tower Environ.—1978, Univ. of Maryland,
College Park, Maryland.
24! Kannberg, L. D. and L. R. Davis. 1976. An Experimental/Analytical In-
vestigation of Deep Submerged Multiple Buoyant Jets. U.S. Environmental
Protection Agency, Corvallis Environ. Res. Lab. EPA-600/3-76-001. Cor-
vallis, Oregon.
25. Hirst, E. A. 1971. Analysis of Round, Turbulent, Buoyant Jets Dis-
charged into Flowing Stratified Ambients. U.S. Atomic Energy Commission,
Oak Ridge Nat. Lab., Rep. ORNL-4685. Oak Ridge, Tennessee.
26. Hirst, E. A. 1971. Analysis of Buoyant Jets within the Zone of Flow
Establishment. U.S. Atomic Energy Commission, Oak Ridge Nat. Lab., Rep.
N. 0RNL-TM-3470. Oak Ridge, Tennessee.
-------
Guthrie, D. L. 1975. Documentation for OUTFALL: A Computer Program for
the Calculation of Outfall Length Based on Dilution Requirement. USEPA,
Region II. New York, New York.
88
------- |