|
A COMPUTERIZED METHOD OF EVALUATING TRACE
PARTICLE RECORDS OF TURBULENT FLOW FIELDS
By
GARY RALPH JACKMAN
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1976
ACKNOWLEDGMENTS
I wish to express my gratitude to Professor E. Rune Lindgren,
my supervisory chairman, for his guidance throughout this work. His en-
thusiastic support and technical assistance were invaluable. Special
thanks for the many hours devoted to the written work of an inexperienced
author. I wish to thank Dr. Kurzweg for his stimulating support and
advice on a near daily basis. I would like to thank Dr. Boykin for his
overall interest and in particular for his help with digital filtering.
I would also like to express my thanks to Drs. Keefer and Roan for their
encouragement and assistance.
My thanks go to Dr. Elkins for his help and friendly criticism
during most phases of this work. I thank Dr. Yoo for his work with the
Pictorial Data Acquisition Computer and associated computer programs.
My thanks go to Randall Crowe who spent much time helping apply the data
tablet method of digitizing to this work and helping me while I learned
to program. Special thanks to Terry Lapinsky, who typed this dissertation.
Words fail to express my appreciation to these people who helped
so much.
This work was supported in part by the National Science Foundation
under Grant No. ENG75-03470 and in part by the Office of Naval Research
under Grant No. N00014-75-C-1090.
TABLE OF CONTENTS
ACKNOWLEDGMENTS...; .......................................
LIST OF FIGURES...........................................
ABSTRACT........... .....................................
CHAPTER
1 INTRODUCTION.......................................
1.0 Techniques for Measuring Turbulence..........
1.1 The Problem and the Proposed Solution.........
2 EXPERIMENTAL ARRANGEMENT..........................
2.0 Introduction................................
2.1 Pipe Flow Apparatus..........................
2.2 Data Acquisition Equipment...................
3 METHOD OF DATA ACQUISITION.........................
3.0 Introduction.................................
3.1 Trace Particle Method of Flow Visualization...
3.2 Prism Analysis and Calibration................
3.3 Operation of the Experimental System..........
4 METHOD OF DATA REDUCTION..........................
4.0 Introduction.................................
4.1 Pictorial Data Acquisition Computer Method....
4.2 Data Tablet/Digitizer Method.................
4.3 Method of Following Particle Images and
Matching Image Pairs........................
Page
ii
v
vii
4.4 Method for Determination of Three Di-
mensional Positions and Velocities............ 57
4.5 Referencing Data to a Grid................... 75
4.6 Method for Displaying Reconstructed
Flow Fields.................................... 81
4.7 Comments............................................. 97
5 THE STRUCTURE OF TURBULENCE AT LOW REYNOLDS
NUMBERS............ .............................. 103
5.0 Introduction................................. 103
5.1 Digital Filter Cutoff Frequency Determination. 103
5.2 Observations of Pipe Flow in the Transition
Zone........................ .................. 111
5.3 The Experimental Results of the Pipe Flow at
a Reynolds Number of 3500...................... 115
6 SUMMARY AND CONCLUSIONS............................ 125
6.0 Summary of the Data Acquisition and
Evaluation................................... 125
6.1 Conclusions of the Data Acquisition and
Evaluation................................... 133
6.2 Summary of the Investigation of Turbulence
at Low Reynolds Numbers....................... 134
6.3 Conclusions of the Investigation of Turbulence
at Low Reynolds Numbers....................... 137
APPENDIX A DERIVATION OF THE RELATIONSHIP BETWEEN THE
IMAGE POSITIONS AS VIEWED IN TWO FACES OF A
PRISM AND THE PARTICLE POSITION IN A PIPE..... 139
BIBLIOGRAPHY.............................................. 145
BIOGRAPHICAL SKETCH...................................... 149
Chapter
Page
LIST OF FIGURES
Figure Page
2.1 Experimental pipe flow apparatus................... 15
2.2 Data acquisition section of pipe (test section).... 16
2.3 Relationship of the apex of the viewing volume
to prism size...................................... 17
3.1 Prism, pipe, and camera............................ 27
3.2 Prism arrangement with field of view in flow....... 28
3.3 Prism and pipe................................. 29
4.1 Flowchart of the method of data reduction.......... 32
4.2 Computer program for converting data tablet coded
coordinates to referenced X-Y coordinates.......... 39
4.3 The particle image series startup process........... 48
4.4 Computer program for following image positions and
matching image pairs............................... 49
4.5 Computer program to calculate filtered positions
and velocities.................................... 65
4.6 Computer program to rearrange data by frame......... 74
4.7 Computer program to reference velocity field to
known locations................................... 78
4.8 Three dimensional velocity field................... 85
4.9 R-6 plane segment velocity fields.................. 86
4.10 X-R plane segment velocity fields.................. 87
4.11 Three dimensional particle path display............. 88
4.12 Computer program to construct velocity fields...... 89
4.13 Computer program to construct particle path
display............................................ 95
Figure Page
5.1 Position error vs. frequency........................ 109
5.2 Velocity error vs. frequency....................... 109
5.3 Flow RMS velocity vs. cutoff frequency.............. 110
5.4 Change in radial RMS velocity vs. cutoff frequency.. 110
5.5 Average axial velocity vs. distance from pipe wall.. 122
5.6 Root-mean-square velocities vs. distance from
pipe wall........................... .............. 123
5.7 Mean-velocity distribution.......................... 124
A-1 Prism, pipe, and camera............................. 143
A-2 Prism and pipe...................................... 144
Abtract of Dissertation Presented to the Graduate Council
of the University of Florida in Partial Fulfillment of the Requirements
for the Degree of Doctor of Philosophy
A COMPUTERIZED METHOD OF EVALUATING TRACE
PARTICLE RECORDS OF TURBULENT FLOW FIELDS
By
Gary Ralph Jackman
December 1976
Chairman: E. Rune Lindgren
Major Department: Engineering Sciences
In a Ph.D. dissertation at the University of Florida in 1974,
Johnson used a prism to view particles tracing the motion of fluid in
a pipe. The particles viewed in the prism were recorded with a cine-
matographic camera. The prism produces two images for each particle,
and geometric relationship uniquely determines the true particle position
from the two image positions. Starting from Johnson's manual method,
computerized evaluation techniques have been developed to digitize
the cinegraphic records of the flow entrained trace particles. The
digitized data are processed through computer programs developed to
follow the recorded particle images as they traverse the prism faces
and pair the image paths that correspond to particles traversing the
studied flow volume. These data are further processed to yield three
dimensional positions and velocities. Computer-graphic methods are
developed and used to display the instantaneous flow in two or three
dimensional views.
CHAPTER I
INTRODUCTION
1.0 Techniques for Measuring Turbulence
The study of the structure of turbulence is a problem which has
attracted the collected efforts of the scientific community for more
than 100 years. The difficulties have stimulated the development of
a variety of techniques which have been used with varying degrees of
success and range from the simple qualitative but very informative
trace-dye observations by Reynolds, 1883, to the quantitative, sophis-
ticated laser measurements of today.
Anemometer probes which are inserted into flow fields to measure
average and fluctuating velocities record changes in electrical,
chemical, or mechanical probe parameters caused by the fluid flowing
next to the probe. Besides the requirement that the measuring probe
must cause a minimum of interference with the flow field, there are a
number of specific restrictions for this type of equipment to be used
successfully in measuring turbulence. Some of the crucial requirements
have been listed by Hinze, 1975. The probe must be as small as possible
in order to minimize flow field distortions induced by it. There must
not be a velocity gradient in the region of the probe sensor therefore
the sensor must be smaller than the dimension of the smallest turbulent
"eddies". (The smallest length scale of the turbulence is 0.25 mm at
Re = 8000 and 0.8 mm at Re = 2200 in a 12.7 cm diameter pipe. The
Reynolds number is Re = U D/v where U is the average pipe flow
velocity, D is the diameter of the pipe, and v is the kinematic viscosity.)
The inertia of the instrument must be low enough for it to respond to
the highest frequencies encountered in the turbulent flow. The sensing ele-
ment must be sufficiently rigid to prevent flow induced vibrations while
being sensitive enough to detect the very small fluctuations of velocity
which are typically a few percent of the average flow velocity. The
probe anemometer system must be stable enough to exclude calibration
shifts during a test run. All these restrictions are very difficult for
probe anemometer instruments to satisfy since they require the probe to
be large for rigidity but small for minimizing flow interference and
avoiding velocity gradient distortions.
Up to the present, the hot-wire and hot-film anemometer's develop-
ment and application for quantitative measurements of turbulent flow
have by far surpassed all other probe anemometers for use in gases and
liquids. The vast bulk of quantitative information about the structure
of turbulence has been obtained using these anemometers. The sensing
elements of the hot-wire and hot-film anemometer are heated by electric
current and cooled by the fluid flowing around them. The total heat
transfer depends on 1) the flow velocity, 2) the difference in
temperature between the fluid and the wire, 3) the physical properties
of the fluid, and 4) the dimensions and physical properties of the wire.
At present the constant-temperature anemometer is used almost exclusively.
Here the electric resistance of the sensing element and, consequently,
the temperature are kept constant with the cooling by the flow being
compensated for by an increase in electric current. The operation of
the hot-wire anemometer is described in detail by Hinze, 1975. The
turbulent intensity must be small compared to the average flow velocity
in order for a linearization of the hot-wire response to turbulent
fluctuations to be valid. Errors of up to 10% may be found for tur-
bulent intensities of 20%. Wygnanski, Sokolov, and Friedman, 1975,
recorded turbulent intensities exceeding 20% of the average flow
velocity within the transition range of Reynolds numbers for pipe flow.
Johnson, 1974, recorded axial turbulent intensities 15% of the average
flow velocity at Reynolds numbers of 4000, 6500, and 9000 in a 12.7 cm
diameter pipe. According to Perry and Morrison, 1971 a and b, errors
up to 20% may be associated with the usual static calibration procedures
and up to 20% difference in turbulent intensity may be recorded by
different anemometer systems in the same flow. Even with the calibration
complete and the above mentioned errors known, there are additional
problems associated with the hot-wire and hot-film systems. The cali-
bration can shift with impact or deposit of contaminants on the sensing
element. Changes in the wire properties can occur due to excessive
temperatures. According to observations by Johnson, 1974, hot-film
calibration data indicate the existence of a critical velocity value
below which the hot-film probe is insensitive to velocity change close
to the flow boundaries. The probe size and thermal interference
restrict the proximity of the probe placement to a boundary. It is
notable that the hot-wire and hot-film anemometers can only measure
one or two components of velocity at a single location for each probe.
A more sophisticated system used for flow studies is the Laser
Doppler Anemometer (LDA) which was first used by Yeh and Cummins, 1964.
Here the fluid velocity is determined by recording the Doppler shift
of the scattered laser light beam. The scattering is caused by particles
entrained in the flow which traverse the measuring section. The anemometer
consists of a laser, a beam splitter, a light collector, a photodetector,
and a signal processor. The advantages are that there is no probe and
no flow calibration is required. The Laser Doppler Anemometer measures
the flow velocity within a volume which is usually somewhat larger than
the linear dimensions of a hot-film probe. Noise is generated by a mean
velocity gradient in the measuring volume when particles traverse dif-
ferent sections of the measuring volume at different velocities. Noise
is also generated by phase and amplitude variations in the Doppler shift
from multiple particles traversing the measuring volume at the same time.
Thus, low particle densities are necessary. Laser systems are very
expensive. It is a difficult task to focus the crossed laser beams
and to determine the location of the measuring volume in a pipe. Two
crossed beams allow only one component of velocity to be determined and
at only one location at a time.
Some other techniques for studying turbulence utilize various means
of making the flow visible to a detecting apparatus outside of the flow
field. These techniques run the spectrum from qualitative through quan-
titative observations. The discussion is limited to methods which make
the flow of liquids visible. The methods utilize discrete particles
called trace particles or"continuous"materials such as dyes, powders, or
birefringent solutions.
The transition experiments done by Reynolds, 1883, utilized trace-
dye to visualize the intermittency of the onset of turbulence in pipe
flows. The first experiments using aluminum powder on free surfaces of
water for studying channel flows were done by Ahlborn, 1902. Prandtl
and Tietjens, 1934, used this same technique in their experiments on
flow around cylinders and airfoils. Prandtl, 1904, suspended small
5
flakes of mica in water to visualize regular motions since a large
number of flakes would orient themselves to the flow direction. Binnie
and Fowler, 1947, and Lindgren, 1954-1969, studied the transition process
in pipe flows of water by addition of birefringent colloidal suspensions,
viewing the flow by means of crossed polarization sheets. This tech-
nique is very sensitive for the discrimination between laminar and tur-
bulent flow but does not allow any quantitative measurements of the
structure of turbulence. Similar observations can be made by mixing
aluminum powder in the flow as done by Coles, 1965. These are but a few
of the many continuous or near continuous flow visualization methods
that have been used to visualize the flow of liquids. They are all
limited to qualitative information or, at best, some gross quantitative
information.
Another flow visualization technique is the hydrogen bubble method
which was first described by Kolin, 1953. Clutter, Smith, and Brazier,
1961, applied this technique to the study of flow around aircraft models.
Schraub, Kline, Henry, Runstadler and Littell, 1965, achieved a major
advancement in the use of the hydrogen bubble technique for quantitative
studies by pulsing the voltage and insulating spanwise portions of the
wire to produce combined time-streak markers. Schraub et al., 1965,
give a detailed description of the measurement of velocity using the
hydrogen bubble technique and an error analysis of the velocity. Kline,
Reynolds, Schraub, and Runstadler, 1967, and Kim, Kline, and Reynolds,
1971, made use of this technique for obtaining qualitative and quantitative
information on turbulence in a boundary layer. The hydrogen bubble
technique is limited to low speeds and low frequency turbulent fluc-
tuations. Only one component of velocity at several locations along
the electrode can be determined at a time. Kline and his coworkers
calculated the fluctuations in velocity by subtracting the hydrogen
bubble determined instantaneous velocity from long time averaged hot-
wire anemometer readings of velocities. Kim et al., 1971, state that
the accuracy of the data and particularly of the fluctuations is low
for this technique. A serious limitation of this method is the bouyancy
of the bubbles which will cause them to rise relative to the flow
pattern they supposedly should record.
The use of discrete particles for flow visualization is also a
well known technique for flow visualization. A method for making spheres
with a specific gravity very close to water was reported by Marey, 1893.
The spheres were a mixture of wax and rosin for use as trace particles
in water. Eden, 1912, injected oil into the fluid by means of an
atomizer and illuminated a thin sheet for visualization of the fluid
flow. Fage and Townend, 1932, used their ultramicroscope to measure
motions close to the wall of a pipe by observing dust particles naturally
abundant in the water. Van Meel and Vermij, 1961, obtained three dimen-
sional position components by using a number of layers of different
colored light stacked normal to the photographic plane. A color photo-
graph was taken and the component of position normal to the plane was
known to the limit of the thickness of the color layer. Nedderman, 1961,
used two cameras placed at an angle to each other to simultaneously
photograph small air bubbles to obtain three dimensional positions and
velocities near the wall in turbulent pipe flow. Nieuwenhuizen, 1964,
used a set of mirrors to put a front and side view of the particles in
the flow on one camera frame. Corino and Brodkey, 1969, used colloidal
sized particles of magnesium oxide suspended in trichloroethylene which
was the experimental fluid. High speed motion pictures were taken of
a very thin planer section as the camera was transported downstream.
Nychas, Hershey, and Brodkey, 1973, used pliolite trace particles in
water in the same apparatus as Corino and Brodkey, 1969.
Johnson, 1974, introduced a prism for studying the motion of
pliolite trace particles in the flows of water through pipes. This
technique permits determination of the three dimensional position of
each particle in the flow field viewed through the prism. The prism
produces two images for each particle and a geometric relationship
uniquely determines the true particle position from the two image
positions. There is no probe to disturb the flow field. The particles
are almost neutrally bouyant and are chosen small enough to follow the
fluid motion very closely. Recording the image positions by means of
a frequency monitored cinematographic camera then permits a quantitative
as well as a qualitative three dimensional reconstruction of the studied
flow field and its variation in time. The changes in particle position
can be progressively more precisely followed by using a higher film
rate (frequency) to record them. No calibration is required and the
data recording cannot shift if the equipment is undisturbed. The prism/
trace particle method seems potentially to be superior to any other
technique presently in use to describe the instantaneous characteristics
of flow fields and seems to imply a major step forward in the field of
experimental fluid dynamics. As far as this author knows, this is the
only method of acquiring a quantitative record of an entire flow volume
simultaneously through time.
1.1 The Problem and the Proposed Solution
The problem now is to implement the data in the record of observations
in order to utilize the potential of Johnson's prism/trace particle
method. This necessitates some kind of automation of the processing
of the cinematographic records. The sheer volume of stored data prevents
any manual evaluation of the records should all the information stored
be utilized for analysis of the three dimensional flow field. The
present work is limited to the development of techniques for numerical
reconstruction of the three dimensional flow field being studied, in-
cluding automated processes for its graphical display in a variety of
fashions. It should be mentioned that the present method of data
reduction can be used with very little modification to process the data
contained in the records of the experiments by Nychas, Hershey, and
Brodkey, 1973. Because of the magnitude of manual work required they
could not previously evaluate these data quantitatively.
The method of data acquisition described by Johnson, 1974, has
been used with only slight modifications. The cinematographic records
of the flow entrained trace particles have been translated into digital
data by two methods. In the first method, a Pictorial Data Acquisition
Computer (PIDAC) is utilized. It automatically scans light which passes
through the motion picture film and records the light intensity level at
various locations on the film. This information is thereafter processed
to yield the digitized location of each particle image. In the second
method, a Scriptographics data tablet/digitizer is utilized. Each cine-
matographic frame is enlarged and displayed on the tablet for digitization
one by one. The translated data from both methods has been processed
through computer programs that developed to follow particle images as
they traverse the experimental area on the prism faces and pair the
image paths that correspond to particles traversing the studied flow
9
volume. These data were further processed to yield true positions
and velocities. Computer graphics methods have been developed and used
to display the instantaneous flow field in two and three dimensional views.
CHAPTER II
EXPERIMENTAL ARRANGEMENT
2.0 Introduction
The experimental equipment used for this work is located in the
Fluid Mechanics Laboratory of the Engineering Sciences Department at
the University of Florida. It consists of two parts: the constant-
head, closed circuit pipe flow apparatus which uses water as its ex-
perimental fluid; and the prism arrangement or the data acquisition
apparatus which consists of a prism, a light source, and a movie camera.
2.1 Pipe Flow Apparatus
The flow apparatus, equipped for flow studies by means of the
prism/trace particle method,is shown in Figure 2.1. The first of its
components, the constant-head overflow tank (1), is simply an open
topped container within a container. The fluid is routed to the con-
tainers through a polyvinyl chloride pipe from the pump (7). The flow
must be slightly greater than the flow through the test pipe (3) in
order to ensure a constant head on the system. The inner container
is thus kept full with the excess fluid overflowing into the outer
container which drains into the mixing tank (8). The fluid in the inner
container travels down the third polyvinyl chloride pipe to the settling
chamber (2) where it passes through a honeycomb rectifier. The fluid
then passes around a flat plate perpendicular to the flow with approxi-
mately two millimeters of clearance thus ensuring a thorough mixing of
the fluid at the entrance of the test pipe (3).
The test pipe (3) is a 23.6 meter long plexiglass pipe.of 127 mm
inside diameter. The test pipe was carefully aligned both longi-
tudinally and horizontally to ensure a straight horizontal test section.
It has eleven permanent joints and five flanged, o-ring-sealed joints
which were carefully matched. The inside diameter of the test pipe
was checked and found to vary by less than + 0.5 mm. The data ac-
quisition equipment is located in the section of the test pipe farthest
downstream, labeled (3a) in Figure 2.1. The fluid exits the test pipe
through a plexiglass access chamber (4), whose dimensions are 0.45
meters in length and .203 meters in diameter. A mount for a hot-film
anemometer probe and a mirror to reflect light from a light source up
the test pipe parallel to the longitudinal axis are installed in the
access chamber. There is also a thermometer to monitor the fluid
temperature.
The flow exits the access chamber, passes through the return line
(5) which is a 75 mm inside diameter plexiglass pipe, and reverts to
the pump to be recycled through the system. Two Ramapo Mark V tur-
bulent drag type flowmeters (models V-3-SP and V-1-SP) each upstream of
a regulating valve are located in the return line. The two flowmeters
and regulating valve systems (6) are arranged in parallel, which allows
the use of either system. The Mark V-3-SP has a range of 0.0013 to
0.013 m3/s for high Reynolds number flows and the Mark V-1-SP has a
range of 0.00013 to 0.0013 m3/s for low Reynolds number flows. This
arrangement allows accurate measurements over a wide range of flow rates.
The strain gauge in the Mark V-1-SP flowmeter is connected to a Wheatstone
bridge and a Sola DC power supply. The output is connected to a Honeywell
Electronik 194 potentiometric strip chart recorder and a Hewlett Packard
Model 412 DC vacuum tube voltmeter. The flow measuring unit was calibrated
by the bucket, scale, and stopwatch method with a minimum of four
replications of every reading with less than 0.5% error. The scale was
a Detecto floor scale which was checked by three different methods for
accuracy. Downstream of the flowmeters, the fluid exits the return
line into the pump. The pump (7) is a glass lined Gould centrifugal
machine with a maximum flow rate of 0.013 m3/s at 10 m water overhead
pressure. It has a flexible connection with the return line in order
to insulate the flow in the pipe from mechanical vibrations of the pump.
The mixing tank (8) is connected to a storage tank (9) and the overflow
from the constant-head tank (1). The fluid in the mixing tank
maintains a constant pressure on the upstream side of the pump (7) which
helps prevent cavitation in the pump. The constant pressure i n the
mixing tank also maintains a constant pressure on the return line so
as to ensure steady flow conditions in the test pipe (3).
Two cooling elements are located in the system to maintain a con-
stant temperature throughout the operating time for an experimental run
since the pump and the ambient conditions add energy to the system thus
attempting to elevate the temperature of the fluid. One is a Precision
Scientific Company portable cooler with its cooling element located in
the constant-head tank (1). This cooler is not fluid temperature re-
gulated. The other is a Blue M 250 W cooler with its cooling element
and sensor located in the mixing tank (8). This is a variable setting
fluid temperature regulated cooler.
The water used as the experimental fluid is from the city water
supply. In order to ensure that the water is clean, an Aqua-Pure water
filter model P10-1 was installed on the city water tap. Six deaerating
taps allowed for the elimination of entrapped air in the experimental
system. Of the six, four taps are manually regulated. Two of these
manually regulated taps are located on the test pipe (3); one is situated
on the access chamber (4); and one tap is on the return line (5). The
two remaining taps are open bleed taps. One open bleed tap is located
in the settling chamber (2), and the other is situated in the mixing
chamber (8). The chlorine in the city water supply cannot attack the
system since all metal parts are fiberglass coated except the Ramapo
flowmeters which are stainless steel.
2.2 Data Acquisition Equipment
The data acquisition equipment (located at 3a in Figure 2.1) is
shown in Figure 2.2. The light source is a Beseler slide projector
with a 750 W bulb. The collimated beam from the projector is reflected by
a mirror aligned 450 to both the projector and the flow. This allows
the beam to pass through the prism viewing volume and parallel to the
flow. The end of the test pipe was painted black to prevent the light
from traveling in the plexiglass pipe which would raise the level of
the background light intensity. A 35 mm model 71 Bell and Howell movie
camera was used to film the pliolite trace particles in the fluid as
seen through the two faces of the glass prism.
Two glass prisms mounted side-by-side were used in this work. The
first prism is the same one used by Johnson, 1974, in his dissertation
work. It is a 68 mm base, 900 prism, 56 mm long. Figure 2.3 shows that
the apex of the viewing volume is just prior to the centerline of the
pipe. This is satisfactory for turbulent flow at low Reynolds numbers
but it was felt that for transition flow, the apex of the viewing volume
should be such that a large area in the vicinity of the centerline of the
pipe is in view. For this reason, another prism was mounted on the test
pipe. The larger prism is a 113 mm base, 160 mm long, 900 crown glass
prism made by Ealing Optics Company. This prism has the apex of its
viewing volume greater than two pipe diameters from the pipe wall/
prism in water. A complete analysis of the prism and pliolite trace
particle method will be developed in Chapter III.
To have a complete acquisition of data, the movie frames must
have a reference point to locate the axial coordinate origin and two
reference points to line up the prism vertical. The frame rate or
film speed must also be known. The referencing is accomplished by using
two wheat grain bulbs of low light intensity which are reflected off
the prism to give two reference spots aligned with the pipe axis. One
spot will define the horizontal position of the frame and the two to-
gether will define the rotation or misalignment of the frame. The frame
rate is acquired by filming a Hewlett Packard 5243L electronic counter.
Li]
rJLL'
1 Constant-head overflow tank 5 Return pipe
2 Settling chamber 6 Flowmeters and
-Test pipe regulating valves
3 Test pipe
7 Pump
3a Data acquisition sectionump
8 Mixing tank
4 Access chamber 8 Mixing tank
9 Storage tank
FIGURE 2.1 Experimental pipe flow apparatus
Prism ,
Castor
oil
Pipe
FIGURE 2.2 Data acquisition section of pipe (test section)
113 mm
68 mm
Large 4
Prism
A I
I
Lr
I
I
I
r
I .
I
V
Apex of the viewing
volume for
small prism
Apex of the viewing
volume for
large prism
FIGURE 2.3 Relationship of the apex of the viewing volume to prism size
/ \
Small \
Prism
/
I
P
CHAPTER III
METHOD OF DATA ACQUISITION
3.0 Introduction
The data acquisition phase consists of all steps required to
cinematographically record the flow field pattern. This is accomplished
by adding trace particles to the fluid and photographing the flow using
a method that will relate the cinematographic film to the three dimen-
sional flow field. The method of data acquisition used in this work was
first developed by Johnson, 1974. It employs pliolite particles as the
flow visualization material and a prism to relate the two dimensional
cinematographic film to the three dimensional pipe flow. The method
will be developed in three parts. The first part will develop the
trace particle method flow visualization. This will be followed by
the prism analysis and calibration. Lastly, the operation of the data
acquisition apparatus and the experimental flow apparatus will be discussed.
3.1 Trace Particle Method of Flow Visualization
The trace particles used in this method of flow visualization are
pliolite supplied by Goodyear Chemicals, Akron, Ohio. The pliolite,
pebble sized when received from the supplier, was ground in a colloidal
mill and sieved to yield particles of 48 to 63 pm in size. These particles
were then mixed in water and left for 24 hours. The particles that settled
and the particles that floated were removed. The water with the neutrally
buoyant particles was added to the system through the constant-head tank.
It was decided to try adding the particles to the system without going
through the tedious and time consuming task of mixing and letting the
mixture stand. Within one hour, the particles from both methods acted
in similar fashion in the flow system. The density of pliolite is
1.026 g/cm3. The particles responded the same whether mixed in water
or dry when added to the system. The particles selected were small
enough to follow closely the fluid motion and large enough to be photo-
graphed with a movie camera under the available lighting. Johnson
justifies the use of the pliolite particles in a lengthy argument
(Johnson, 1974, pp. 38-39) which will not be repeated here, although a
comment on this analysis seems appropriate. Johnson uses the equation,
U' (aT + 82)
U' (aTL + 1)
00
where a = 2p+ ) d2 3p ) ,TL = Lagrangian time scale =RL(-r)dT,
=(2p+ p) d2 =(2p + p) L
P P a
U' is the root-mean-square value of the particle motion, U' is the root-
p
mean-square of the fluid motion, p and p are the densities of the particle
and the fluid respectively, p is the fluid viscosity, and d is the diameter
of the particle. This equation is used to show that the pliolite trace
particles indeed follow the fluid motion. Johnson notes that aTL is very
large even for small values of TL. He mentions also that B is very close
to one. Thus, with aTL very large and 8 approximately one, the pliolite
particles follow the fluid flow. However, it was found that the actual
value of a for pliolite and water at 20C is 0.98 which shows his use
of 8 equal to one is very close. Even if the quantity aT was equal to
L
zero, U' would still be 96% of U' which is still reasonably close to the
fluid motion.
fluid motion.
3.2 Prism Analysis and Calibration
The two prisms were mounted on the plexiglass test pipe with castor
oil filling the gap between the flat face of the prism and the curved
pipe wall as shown in Figure 2.2. The intersection of the two smaller
faces of the prism was parallel to the longitudinal axis of the test
pipe. Castor oil was chosen to fill the gap because its index of re-
fraction (1.48) is very close to the plexiglass (=1.49) and crown glass
(1.52). A common index of refraction of 1.50 was chosen for the prism,
castor oil, and pipe unit for ease of analysis since very little error
was introduced.
Viewing the trace particles in the fluid through the prism yields
two images for each particle as seen in Figure 3.1 from the camera position.
The distance between these images is related to the Y coordinate while
the difference in positions of these images from the two face inter-
section is related to the Z coordinate. The X coordinate is measured
relative to any fixed reference location along the longitudinal axis of
the pipe. The general field of view in the fluid as seen through the
prism is sketched in Figure 3.2. This flow field volume cross section is
a circle sector with its apex from the pipe wall determined by the prism
dimensions as detailed in Chapter II and partly shown in Figure 2.3.
For a prism mounted on the test pipe illustrated in Figures 3.1 and
3.2, the complete derivation of equations for relationship between
the image positions and the particle position is developed in Appendix A.
The origin of the Cartesian coordinate system is at the center of the pipe
(Figures 3.1(A-l)and3.3(A-2)). The equations are written for two paths
from point A to point P (ABCP and AEDP) by writing the equations of the
individual lines and intersecting them as shown in Figures 3.1 and 3.3.
Introduction of the angle < allows the derivation to compensate for
a misaligned camera. Angles i and j are angles of incidence. Since
the pipe wall is curved, the normal to the pipe wall associated with
angle j changes with respect to the Cartesian coordinate system as the
image position in the prism moves away from the Y axis (as B or E gets
further from the center of the prism).
The results of the derivation in Appendix A are equations which
locate the particle position in the pipe based on ZB, ZE, the indices
of refraction, and the geometry of the pipe-prism unit. The equations
are
Zc Cot(ec + v) + z Cot(e + Y) /R2 2 /R2 Z
P Cot(OC + YC) + Cot(eD + YD)
(3.1)
Y = /R2 Z + ( -C) Cot(C + ) (3.2)
where
C =Tan ZC = Tan ZD
2_I 2 / 7-2_
RC ( ) D (
N N
YC =Sin-[ Sin( n e Sin-l[ a in( 0)])]
w g
N N
Y = Sin-l[ N- Sin( 0 Sin[ -a Sin( 5)])]
w g
-b+kA2 4ac
ZC 2a
with
a = Tan2( + 6C) + 1
22
b = 2[h ZB Tan( 2)] Tan(n + 6C) 2ZB Tan2(+ + 6)
B_ T C B C2
c = [h ZB Tan( )]2 R2 2[h ZC Tan( Q- -)]ZB Tan(Q2 + 6 C)
+ Z Tan2 ( + 6C
S -e /e2 4df
ZD =
2d
with
d = Tan2(( + 6 ) + 1
D
e = 2[h ZE Tan( 0)] Tan(Q + 6D) -2ZE Tan2(Q + 6D)
T 2 2 _
f = [h ZE Tan( )] 2[h ZE Tan( )]ZE Tan(S 6D)
+ ZE Tan(S + 6D)
E D
6C = Sin-l[
N
-1 a
6 = Sin [ -
g
6 Sin1 a
D N
g
Sin( 1- + ()]
Sin( 0 )]
2
These equations are valid for any position of a particle in the
prism field of view. The equations relating the image positions to the
flow space particle position were then written into a computer program
which was checked for validity by placing a grid with line intersections
at known locations in the fluid behind the prism. This validity check
and
was done for both the large and small prism. Each was checked in two
ways. In both cases the first way wasbyamanual graphics method. The
second way was by use of the method of data reduction utilized to reduce
the data acquired through that prism. The second method for the large
prism was done by using the Pictorial Data Acquisition Computer, while
the second method for the small prism utilized the data tablet/digitizer.
Both second methods will be discussed in detail in Chapter IV. The first
method (manual graphics) was the same method used by Johnson, 1974, to
check the validity of the equations for the small prism. This method
consisted of filming the prism with the grid in the fluid behind the
prism, mounting the film in a slide projector, projecting and tracing
the grid from the film to graph paper, measuring the positions of the
intersections, and processing the data through the prism equations com-
puter program with the appropriate scale factor. The scale factor for
the present work was determined for the change in the X coordinate, which
was unaffected by the prism, and was checked by relating the graph paper
prism size to the actual prism size. Table 3.1 shows the results. The
calculated positions agree with known positions although there were some
noticeable differences for small values of the Y coordinate (farthest from
the prism) and for large absolute values of the Z coordinate (vertical
distance from the XY plane). These discrepancies can be attributed to
the use of a single index of refraction for the prism-pipe unit and the
fact that small errors in reading become slightly larger errors in true
position with a decreasing Y coordinate and an increasing Z coordinate.
3.3 Operation of the Experimental System
This section describes the operation of the experimental flow
apparatus and the data acquisition equipment as an integrated unit.
The discussion will encompass the general operation of this equipment.
It will discuss the start up and running, the data acquisition, and
the shutdown. Figure 2.1 shows the experimental flow apparatus and
Figure 2.2 shows the data acquisition equipment.
When using the experimental pipe flow apparatus, the fluid level
tube on the side of the storage tank (9, Figure 2.1) should be checked
before the pump is turned on. This will ensure adequate fluid for
operation of the pipe flow apparatus. The proper fluid level is marked
on the fluid level tube. The electronic equipment must be warmed up
prior to use for 15 to 20 minutes (vacuum tube voltmeter, strip chart
recorder, and flowmeter). Both cooling units (cooling elements located
in 1 and 8, Figure 2.1), have to be turned on and the desired temperature
set on the Blue M unit rheostat. The pump (7, Figure 2.1) is then turned
on with both control valves (6, Figure 2.1) closed. When the fluid is
seen in the sight glass on the overflow line from the constant-head
tank (between 1 and 8, Figure 2.1), the large line control valve is
opened quickly to pick up any pliolite particles that have settled out
of the fluid from previous experiments. The valve is then closed and
the strip chart recorder and vacuum tube voltmeter are "zeroed". At
this point, the temperature of the fluid is noted (4, Figure 2.1) and
the flow rate is adjusted for the desired Reynolds number using a graph
plotted for this purpose. As precautionary measures during each experi-
mental run,check the strip chart potentimetric recorder for a constant
reading, the sight glass for fluid (constant head), and the thermometer
for a constant temperature.
Caution in the data acquisition procedure can reduce or eliminate many
problems in the data reduction phase. For maximum light intensity, the projector
light beam must be collimated in such a fashion that at the prism
location it illuminates only the convergent flow region viewed through
the prism as discussed in Chapter II. If the beam illuminates any
volume outside of the flow region to be mapped, particles will be
illuminated which can be seen through only one face of this prism.
This will produce series of single view positions of particles that
cannot be matched to images in the other prism face rendering the
matching of image pairs much more difficult. The limited illumination
can be accomplished by using a shaped hole in a black slide in the
projector. With the proper illumination accomplished, the movie
camera is set up. For the small prism, a 152 mm telephoto lens with
a close up adaptor, and for the large prism, a 100 mm telephoto lens
with a close up adaptor were found to fit adequately the vertical
dimension of the prism onto the movie film. For the proper depth
of field, an aperture setting of 11 was used for the small prism and
an aperture setting of 16 was used for the large prism. The electronic
counter was filmed before and after the film record was made in order
to check its rate and continuing constancy.
The shut down of the experimental pipe flow apparatus is accom-
plished by turning off the pump and coolers. This was done without
closing any valves before the fluid level tube on the storage tank
indicated that the quiescent water level is again attained.
TABLE 3.1 CALIBRATION VALUES FOR PRISM ANALYSIS CHECK
Data for large prism:
True (Z) = 0.00
Calculated (Y)
63.5
63.4
59.9
54.9
49.4
45.1
39.5
34.8
30.0
20.0
15.3
11.1
7.2
0.1
0.0
0.0
0.0
0.0
0.0
0.0
7.05 14.10 21.15 28.20
Calculated (Z)
14.2
14.3
14.1
21.2
21.5
True (Y)
63.5
63.4
60.0
55.0
50.0
45.0
40.0
35.0
30.0
20.0
15.0
10.0
5.0
0.0
28.2
28.6
28.0
Data for small prism:
63.5
60.0
54.8
50.0
44.2
40.0
30.0
25.0
20.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
7.0
7.2
7.1
7.1
7.0
6.9
6.7
14.3
14.0
By tablet:
63.5
59.8
55.0
50.0
45.0
40.0
34.8
29.9
24.6
19.8
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
7.0
7.2
7.2
7.1
7.1
7.1
7.0
6.9
7.0
6.8
14.3
14.3
14.2
14.1
14.1
14.1
14.4
21.0
63.5
60.0
55.0
50.0
45.0
40.0
30.0
25.0
20.0
63.5
60.0
55.0
50.0
45.0
40.0
35.0
30.0
25.0
20.0
FIGURE 3.1 Prism, pipe, and camera
Field of view
Shell section for velocity averaging -
FIGURE 3.2 Prism arrangement with field of view in flow
FIGURE 3.3
Prism and pipe
CHAPTER IV
METHOD OF DATA REDUCTION
4.0 Introduction
Data reduction begins by translating the cinematographic record
of trace particle images to digitized data for the IBM 370 Computer.
This has been accomplished using the Pictorial Data Acquisition Computer
and the Scriptographics data tablet/digitizer. The output of both of
these methods is a set of two dimensional particle image locations
as recorded on the film frames. These image locations grouped by frame,
served as the input data for a program which followed the images from
frame to frame as they traversed the fluid volume visible through the
prism, and matched each pair of the two particle image series that
corresponded to a particle entrained in the flow. The data were then
processed through a program that solved for the three dimensional posi-
tions of the particles, filtered out noise associated with the conversion
of data from cinematographic records to computer recognizable locations,
solved for velocities, and interpreted the data statistically. Part of
the output was the positions and velocities of the particles as they
moved through the prism visible region. These data arranged by particle
path, were used in a program that plotted a three dimensional perspective
view of the reconstructed particle paths and in another program that re-
arranged the data by frame (constant time grouping) for use in a display
of instantaneous velocity fields. The input to the velocity field display
program was either directly from the rearranged data or the output of a
program that interpreted the velocities at known grid pLnts from the
rearranged particle position data. The general flow;of data through
the method of data reduction is flowcharted in Figure 4.1. The number
beside each step refers to the corresponding section in Chapter IV.
4.1 Pictorial Data Acquisition Computer Method
The first method of converting cinematographically recorded data
to digitized data recognizable by the computer was a method using the
Pictorial Data Acquisition Computer (PIDAC) located in the Picture Pro-
cessing Laboratory at the University of Florida's Center for Information
Research. The PIDAC is a flying spot scanner system interfaced with an
IBM 7094-11 computer. The PIDAC was designed for biomedical research
but it possibly seemed to be adaptable for processing trace
particle image records into digital information for use in the IBM 370
computer operated by the Northeast Regional Data Center located at the
University of Florida.
The flying spot scanner consists of a high intensity light source,
a converted Nikon still camera, and an image scanner. The high intensity
light was directed through the lens of the Nikon camera onto the film
frame which was placed in the camera. The light that passed through the
2
film was picked up by the image scanner. As it was set up, the 35 x 24 mm
still camera slide was scanned in 648 lines along the 35 mm side and 464
lines along the 24 mm side. Since the 35 mm cinematographic camera
takes half still camera frame sized frames, two frames could be processed
at one time. This meant that the effective scanning of one movie frame
was 464 lines by 324 lines. This corresponds to a spatial sampling of
approximately 0.05 mm on the film or approximately 0.27 mm on the large
prism which was used to acquire data for this method.
DATA ACQUISITION
(Chapter III)
PICTORIAL DATA
ACQUISITION COMPUTER (1)
I I
PARTICLE IMAGE FOLLOWING AND IMAGE
SERIES MATCHING (3)
3-D POSITION AND VELOCITY (4)
DATA TABLET/
DIGITIZER (2)
1
OUTPUT: average
velocities and RMS
values
PARTICLE PATHS (6)
OUTPUT: 3-D perspective
view of particle paths
.REPUNCH BY
FRAME (4)
) i --
GRID REFERENCING (5)
OUTPUT: grid position
and velocities
VELOCITY FIELD (6)
OUTPUT: 3-D perspective
velocity fields, 2-D R-e
2-D X-R views
view of
views, and
Flowchart of the method of data reduction
FIGURE 4.1
The output of the scanner could be viewed on a storage cathode ray
scope before it was processed through the IBM 7094-11 computer and re-
corded on magnetic tape. This visibility was useful in determining the
optimum intensity level of the light source to yield the highest quality
particle image output possible. Since the storage scope has a memory,
several frames could be displayed simultaneously to ensure that the
particle image could be seen from frame to frame as it traversed the
viewing area. In this step, an eight level (0-7) light intensity output
was associated with each scan intersection.
This was a tremendous amount of information since there were more
than 300,000 locations each with eight computer bits of information re-
quired for every two cinematographic frames scanned. The data could be
reduced by a factor of eight by putting this information in binary form
since only one bit would be required for each scanned intersection.
This was acceptable because the only information necessary was whether
or not aparticle image was at the intersection, not its intensity level.
This was accomplished by processing the data through a binary conversion
program using a threshold intensity level. If the intensity number was
at or below the threshold level (6), the number "1" would be associated
with that location, while above that level the location associated number
would be "O". The output of this process was put on computer tape.
The final part of the PIDAC phase of data reduction was to process
the binary data through a computer program that ultimately yielded Cartesian
coordinate locations for the center of each particle image. The program
encoded the boundary of each particle and solved for its center location
in Cartesian coordinates. Two prism fixed reference marks were included
on each cinematographic frame. These marks were used in the next section
of the program to orient the images of each frame to an XY coordinate
system which was nonvarying with respect to the prism. This was ac-
complished by designating the left reference mark center as the origin,
and rotating all the data around this origin so that the right reference
mark center was on the X axis. This operation corrected the data for
mispositioned and misaligned film in the camera film mount of the flying
spot scanner or a misaligned cinematographic camera. The output, then
referenced and rotated, was recorded on computer cards. These computer
cards formed one of the input data sets for the program described in
the method of following particle images and matching image pairs.
The major advantage of this system is the speed and ease with which
the data were processed. The method utilized by Johnson, 1974, required
manually mounting the frames in slideholders, projecting the slides,
marking the particle image centers, following the images, and matching
the image series that corresponded to particles. The automatic scanning
process accomplished the first three steps several orders of magnitude
faster and with considerable less effort. However, there are two major
trouble spots in the PIDAC method which made the method only marginally
capable of yielding useful data. The first problem was the scan interval
since it could not be altered appreciably. There were 424 by 324 scan
positions for each cinematographic frame record corresponding to an in-
terval of 0.05 mm on the film or 0.27 mm in the viewed 113 mm base prism.
The particles used were from 48 to 80 pm in size. The cinematographically
recorded particle image diameters were three to five times smaller than
the scan width. This resulted in the possibility that particle image
center locations could be in error by as much as + 0.1 mm in the X and
Y locations. It was also possible to miss an image position entirely if
the recorded image was less dense than most of the others. A recorded
image could also be lost if it was located near the outer limit of the
scan line width since the sensitivity of the scanner is diminished con-
siderably near the line limits. The missing of a particle image would
interrupt the image series and cause the loss of much data. The accuracy
and data loss problems could be considerably reduced by using a unit
with higher resolution. The second problem with the system was that the
lensing of the flying spot scanner did not allow uniform illumination of
the film. The center of the scanned area received the most illumination
with the intensity decreasing radially toward the boundaries. This
caused a radial decrease in particle image density with progressive loss
of data toward the scan area boundaries. This is a compounding problem
since the particles farthest from the prism in the test pipe, which appear
smaller and less intense on the cine film, were closest to the upper and
lower boundaries of the cinematographic frames and, consequently, the
upper and lower boundaries of the scan area.
4.2 Data Tablet/Digitizer Method
The second mode of converting cinematographic records to digitized
data recognizable by the IBM 370 computer was to utilize the Scripto-
graphic data tablet/digitizer. The movie film was enlarged using a stan-
dard photographic enlarger and projected onto a standard 0.279 by 0.279 m
Scriptographic data tablet/digitizer manufactured by Summagraphics. The
data tablet was interfaced with the IBM 370 computer through a terminal
control program (TCP). The location of each particle image was digitized
by touching the data tablet at the center of each projected image with
the recording stylus. The coded tablet coordinates of these positions
were then transmitted to the computer. The physical principles of the
data tablet's operation are explained in Summagraphics Corporation, 1973, p.2,
Response to a questionnaire from the University of California on April
30, 1973,
Magnetostriction and inverse magnetostriction. The substrate
under the vinyl surface is a magnetostrictive metal. Currents
are pulsed along the X and Y axis, setting up plane strain waves
in the metal. As these waves pass under the small pick-up coil
in the stylus (cursor), the inverse magnetostriction produces an
electrical signal. Digital data is derived from counters which
begin counting when the 'start' currents are pulsed and end the
counting when the 'stop' signal is received from the stylus (cursor).
The relative coordinates of the image centers were verified by viewing
each frame on a Tektronix 4013 cathode ray tube storage scope. The data
were next stored on computer discs and computer cards. The last step of
this operation was to process the raw output (symbols) through a computer
program that converted the symbols to useful coordinate data.
The program took as input the symbols from the data tablet/digitizer
and converted them to locations related to test pipe fixed coordinate
system so that it would not vary from frame to frame. The input symbols
were first processed through the subroutine called "XLT" which translated
them to numbers. The first two data points enteredwere the left and right
reference marks from the film. The left mark was used to reference the data
to the tablet origin and both marks were used to reference the data to the
tablet horizontal. At this point each frame was referenced to the same
origin with the same orientation. This operation corrected the data for
mispositioned and misaligned film in the enlarger or a nonaxially aligned
data acquisition motion picture camera. The data were next recorded on
computer cards and became the second input data set for the program des-
cribed in the method of following particle images and matching image pairs.
The complete program is shown in Figure 4.2.
Although this method of converting visual data to digital data was
tedious and time consuming when compared to the PIDAC method, it was still
several orders of magnitude faster than the manual method applied by
Johnson, 1974. It did not have the lighting and scan width problems asso-
ciated with the PIDAC method. The high intensity light from the enlarger
made the particle images very easily seen on the tablet surface so data
were not lost in the process of digitizing.
The errors in position associated with digitizing data by this
method were the result of the angle of tilt of the stylus from vertical,
the accuracy of the tablet, and the human error of location. The error
caused by tilting the stylus was a line count shift which was two counts
for 220 and four counts for 450. The shift was in the direction of the
tilt only and was a shift not a random error. If, as was the case,
reasonable care was taken to ensure that the sylus was held at approxi-
mately the same angle, no relative shift would be induced between particle
image locations. The absolute shift was compensated for by the fact that
the shift was also in the two referencing marks. The manufacturer
specifies an accuracy of + 0.127 mm for location specification on the
tablet. The data processed with the data tablet/digitizer were obtained
using the small prism. The maximum errors in the specification of location
in the flow region that can possibly result from the data tablet accuracy are
+ 0.035 mm in the longitudinal coordinate, + 0.075 mm in the radial co-
ordinate, and + 0.025 mm in the circumferential coordinate. It should
be noted that the largest error was in the radial direction as would be
expected from observing the image rays in the derivation of the prism
equations in Appendix A since small deviations in image position result
in large radial and small azimuthal deviations in particle position.
The radial root-mean-square velocity i s the primary component used
to confirm the maximum frequency retained after noise filtering reported
in Chapter V, since it contains the largest maximum error and i s the
38
easiest to distinguish. Taking into account the particle image size
as projected on the tablet, the estimated human error in particle center
location was less than one-third of the radius of a particle or + 0.01 mm
in true position in the test pipe. Both machine and human errors
could be reduced to less than one-third by using the now available 0.914
by 1.219 m tablet and enlarging the frame to fit the larger tablet which
has the same + 0.127 mm resolution associated with it.
651-FOLLCW -FORT -SYSPRINT
LEVEL 2.1 ( JAN 75 ) 05/360 FORTRAN H EXTENDED DATE 76.239/20.
REQUESTED OPTIONS: NCOECK.NOLCADCPT=0
OPTIONS IN EFFECT: NAME(MAIN) NOOPTIMIZE LINECOUNT(74) SIZE(MAX) .ATOCEL(NCNE)
SOURCE EeCDIC NOLIST NODECK NOOBJECT MAP NOFORMAT GOSTMT NCXREF ALC NOANSF NOTER
C THIS PROGRAM ANC THE DATA INPUT METHOD MAS DEVELCPEC EY
C GARY JACKMAN AND RANDY CROWE
C THE INPLT IS BY MEANS CF A X-Y TABLET
C THE PROGRAM USES REFERENCE MARKS CN THE FILM 7C ROTATE ALL FRAMES
C TC TABLET HORIZONTAL AND TRANSLATE ALL FRAME BOUNDARIES TC THE
C SAME X-Y COORDINATES
ISN 0002 INTEGER DATA(SOC,4),XTYTXS(500),YS(500)
ISN C003 INTEGER XI.Yl X2,Y2,XYD(4)
ISN 0004 LOGICAL PNCH
ISN 0005 REAL*4 NAME(5)
ISN C006 DATA STOP/'STCP'/
ISN 0007 READ(5.1001) XT,YT
ISN C008 1001 FORMAT(214)
ISN 0009 READ (5.1002) PNCH
ISN 0010 1002 FORMAT(L1)
ISN C011 1 READ(5,1003) NPTSNAME
ISN 0012 1002 FORMAT(14,5A4)
ISN 0013 [F(NPTS.EQO.) GO TO 50
ISN C015 IF (NAME(1).EG.STGP) GC TC 50
ISN 0017 READ(5,1000) ((DATA(I.J).J=1,4)9I=1,NPTS)
ISN C018 IC CONTINUE
TSN 0019 11 CONTINUE
ISN 0020 WRITE(6,1010) ((DATA(I.J),J=1.4).I=1,NPTS)
ISN 0021 101C FORMAT(' 1*//,1OX,'RAW DATA*,/(* -,',20(4A1).'-,*,/))
ISN 0022 100C FORMAT(20(4A1))
ISN 0023 00 21 J=1,4
ISN C024 21 O(J)=DATA(1,J)
C XLT(D,X,Y) DECODES THE TABLET SYMBOLS TO X-Y COORDINATES
ISN C025 CALL XLT(D,X.Y)
ISN 0025 Xil=X
ISN C027 Y1=Y
ISN C028 *RITE(6.9100) X1IYI
ISN 0029 9100 FORMAT(///,10XI'POINT 1: ',2110)
ISN C030 920C FORMAT(1OX,'PCINT 2: ',2110)
ISN 0031 DO 22 J=1,4
ISN 0032 22 C(J)=DATA(2,J)
ISN CC33 CALL XLT(D.X,Y)
ISN 0034 X2=X
!SN C035 Y?=Y
ISN C036 ARITE(5,92O0) X2,Y2 o
ISN 0037 ANG=FLOAT(Y2-Y1)/FLOAT{X2-X1)
.FIGURE 4.2 Computer program for converting data tablet coded coordinates to referenced X-Y coordinates
ISN 0038 ANG=ATAN(ANG)
ISN 0039 %RITE(6,9000) ANG
ISN 0040 90CC FORMATI1OX,'ANGLE: ',F10.4)
ISN 0041 C=CCS (ANG)
ISN 0042 S=SIN(ANG)
ISN 0043 CO 20 I=3,NPTS
ISN 0044 00 25 J=1,4
ISN 0045 25 D(J)=DATA(IJ)
ISN 0046 CALL XLT(DX,Y)
C THIS SECTION ROTATES AND TRANSLATES THE DATA
ISN 0047 XS(I-2)=-((C*(X-Xl)+S*(Y-Y1))-XT)
ISN 0048 YS(I-2)=((-S*(X-X1))+C*(Y-Y1))+YT
ISN C049 20 CONTINUE
ISN 0050 NPTS=NPTS-2
ISN 0051 HRITE{6,1100) NAME,NPTS,((XS(I),YS(I))J,I=1NPTS)
ISN 0052 1100 FORMAT ('1'e//, IX,'FRAME ID = *'5A4,/,
.10X.'NUMBER CF PARTICLES = ',14,
.//,(' 9(2159' ),/))
ISN 0053 WRITE(7. 1201) NPTS,NAME
ISN 0054 1201 FORNAT(I4,5A4)
ISN 0055 WRITE(7.12CC) ((XS(I).YS(I)),I=1, NPTS)
ISN 0056 1200 FORMAT(10(214))
ISN 0057 4C CONTINUE
ISN 0058 GO TO 1
ISN C059 50 CONTINUE
ISN 0060 STOP
-
o
FIGURE 4.2 continued
REQUESTED OPTIONS: NCCECK,NCLOAD,CPT=0
OPTIONS IN EFFECT: NAt'E(MAIh) NOOPTINIZE LINECOUNT(74) SIZE(MAX) ALTODBL(NCNE)
SOURCE EBCDIC NOLIST NODECK NOOBJECT MAP NOFORMAT GOSTMT NCXREF ALC NCANSF NOTER
ISN 0002 SUBROUTINE XLT(DATAX,Y)
ISN 0003 INTEGER DATA(4),CHR(32),XY
ISN 0004 DIMENSION ND(4)
ISN CC05 DATA CHR/ *** ', ', $@*'%' *&** *''**9 )9**** '
.*: *,*. ,* *<. *=*,*>*, *7
ISN 00 006 DO 2 J=1,4
ISN 0007 00 1 I=1,32
ISN 0008 IF (CHR(I).EQ.DATA(J)) GO TO 2
ISN 0010 1 CONTINUE
ISN CO11 2 ND(J)=I-1
ISN 0012 X=ND(3)*32+ND(4)
ISN 0013 Y=NO(1)*32+ND(2)
ISN 0014 RETURN
ISN 0015 END
I-
i-
FIGURE 4.2 continued
OS/360 FORTRAN H EXTENDED
LEVEL 2.1 ( JAN 75 )
DATE 76.239/20.
4.3 Method of Following Particle Images and Matching Image Pairs
As discussed earlier, the particle position in the test pipe is re-
lated to the positions of the two particle images as viewed in pairs
through the faces of the prism. One image is on the upper face and the
other image is on the lower face. Each frame of the cinematographic
record specifies the relative positions of all images of particles
within the test volume of the pipe flow viewed through the prism at a
given instant in time. Identifing each particle image and recording its
position from frame to frame as it traverses the viewed flow field behind
the prism is the first step towards a numerical reconstruction of the
particle path in the test volume. The next step is to pair the series
of particle image positions in the lower prism face with the series of
positions of the particle image in the upper prism face that correspond
to the same individual particles in the fluid. The prism equations
relate the image positions to the position of the particles in the pipe
test volume. The data needed to accomplish the following of the image
from frame to frame and matching image pairs are the particle image
positions from the PIDAC or the Scriptographic data tablet records.
The algorithm for this program imitates the procedure of a human
operator in following the images and matching related pairs. This is
a first generation program so it necessarily imitates only the simplest
human procedures. In succeeding generations, it should be possible to
develop more sophisticated programs. As it stands now, the program
follows the position of the image from frame to frame by using the pre-
ceding two positions to construct a vector to the expected next position.
It then uses a tolerance in both the X and Y coordinates to allow for
deviations from the expected positions. This is done repeatedly until
the image is no longer visible in the prism. The human eye uses this
procedure to follow a moving object. This method is quite successful if
the images are sparse enough and the image paths cross each other infre-
quently. The start up of a series is more difficult since there is no
reference for displacements between frames. This is compensated for by
choosing a distance that the image could not have exceeded between frames,
using the vector construction procedure described above on each image in
this area and finding the second frame image that has a third frame image
within the specified tolerance in X and Y. Next the program scans the
image series in.the upper prism face for that series which corresponds
to the same particle. This is easy to do since it only entails comparing
the X coordinates of the series in the same frame over a sufficient
number of frames to be sure that they are traveling identically in X.
It is possible that if two series should cross, the computer
following series would jump from one particle image series to the other.
If two series or image paths were to come close enough to each other, it
would also be possible to have a series jump. This would result in
losing both series and their matches in the other prism face, or, in
other words, losing two particle paths in the fluid. The higher the
particle density, the more likely it is that this will happen. The
human operator solves this problem by following the series as far as
possible, skips the doubtful section, and constructs new series after
the doubtful area. Next, he constructs the series from the other prism
face record and matches it to the original image path fragments to sort
out the confusion. The second generation program should use a similar
procedure. While the computer matches image pairs it also follows the
series using a trace back method to resolve cross-over problems. This
program should thus be able to handle seven or eight times the data of
the present program. The approximate amount of particle images that
can be handled by this method is at present 150 to 200 which corres-
ponds to 75 to 100 particles in the flow region at any one time. As we
will see later, four or five times this much data are the minimum amount
necessary to give a good description of the instantaneous characteristics
of the flow regions used in these experiments.
Since the program is lengthy and several modifications may be
necessary, the program and modifications are discussed in a brief verbal flow
chart below. The complete program is listed in Figure 4.4. All references
to the program are in sequential order. The first section of the program
follows the image centers from frame to frame as they traverse the flow
field. The first part of the first section identifies an image and
starts the tracking process. The second part tracks the image until it
disappears. The second section pairs the two image series that are
associated with the same particle.
In order to follow a particle image, it must first be identified.
The method of identification or start up for a particle image series
consists of choosing an image position in the first frame in which it
appears and constructing a digital window through the X coordinate of
the first frame image position and the maximum distance the image could
be expected to travel in + X and + Y. This is referred to as "WINDOW 1"
in the program in Figure 4.4. All images in the second frame that
appear in this window are identified as possible second frame positions
for the series. Vectors are calculated between the first frame image
and each of the second frame images. The vectors are then extended from
the second frame positions with the same magnitude and direction. A
second window, "WINDOW 2" in Figure 4.4, is constructed around the end
point of each vector extension. Figure 4.3 shows an example of the
windows and particle image positions. The data from the third frame
are searched for an image position in any "WINDOW 2". If no image
position is found in any of the second windows ("WINDOW 2"), the image
position of the first frame is considered noise and discarded. If one
image is found in a second window, the image position in the first frame
is considered the start of a series, and the image positions in frames
two and three that were located in successive windows are added to the
series and removed from the data. If more than one image position is
found in the second window, the first one tested by the computer will
be chosen. If more than one image position is located in the second
window, there is a high likelihood of error. The computer print-out
will indicated how often two or more images were found in the second
window and which frames. If this happens too many times, the size of
"WINDOW 2" must be decreased. The window sizes must be changed by the
program user to reflect the changes in flow conditions. The higher the
flow velocity the larger the window widths but the smaller the window
heights required. "WINDOW 1" is varied by changing "WIDTH 1" and "HEIGHT
1", and "WINDOW 2" is varied by changing "WIDTH 2" and "HEIGHT 2". Figure
4.3 illustrates the start up procedure. The computer tracks the particle
image through the flow field by constructing a vector from the last
position by an equal magnitude. "WINDOW 2" is digitally constructed
around this expected location. The data of the next frame are searched
for an image position in this area. The image position identified is
added to the series and removed from the data. If no image position is
located in the window, the series is terminated. When all the image
positions of a frame have been added to the proper series and removed
from the data, the remaining image positions are treated as first
positions for start up as previously described.
After the particle image series are constructed, they must be paired
to completely define each particle path in the flow field. This is accom-
plished in the last section of the program which is labeled, "Try to match
a pair of identified series". The series in the lower prism face are
taken one at a time and tested against each series in the upper prism
face. In virtually all cases, the initial frame number of each two
series that successfully paired were within two of each other. Conse-
quently, this was checked first to eliminate series from the upper
prism face that could not be matched, thus saving much computer time.
The X coordinates were compared until a series in the upper prism face
coincided with the series in the lower prism face being tested 15 times.
This was considered to be a sufficient number of matches to determine
the image pairs of a particle. The errors discussed earlier made it
necessary to put a tolerance on the X coordinates of the upper prism
face series. The tolerance used in this case was six data tablet units
which was much larger than the errors but small enough not to cause
pairing mistakes.
Employment of this program may necessitate change in various parts
for compatibility with the data. In the section concerned with following
images, it may be requisite to change the window sizes, the maximum length
of the series, the maximum number of series, and the number of frames of
data to be considered. In that section of program pertaining to image
series matching, it may be essential to alter the number of matches to
consider two series paired, and to change the tolerance on the X coordinates
of the upper prism face series.
The output from this method is the paired series that define the
particle paths in the flow field. The information for each series pair
is the particle number, the initial frame number, and the list of image
positions by frame on computer cards. This information is the input
data necessary for themethod for determination of three dimensional
positions and velocities. The printed output is the image series from
the first section of the program and the paired series output with both
of the X and Y values from the second section. This permits a manual
check of paired series for accuracy and makes it possible to manually
match more series if necessary.
Wi ndow2
-I
r----I
I I
I I
I I
I I
L.---------
Wi ndowl
*Particle position in frame 1
SParticle position in frame 2
0 Particle position in frame 3
Window2
r-----i
S0
Wi ndow2
r ----J
L ..
FIGURE 4.3 The particle image series startup process
r-:- -- -
230-FOLLOW -FORT -SYSPRINT
LEVEL 2.1 ( JAN 75 )
OS/360 =0RTRAN H EXTENDED
DATE 76.251/15.
REQUESTED OPTIONS: NODECKNOLOAD.OPT=0LC(48)
OPTIONS IN EFFECT: NAME(MAIN) NOOPTIMIZE LINECOUNT(48) SIZE(MAX) AUTODBL(NONE)
SOURCE EBCDIC NOLIST NODECK NOOBJECT MAP 4OFORMAT GOSTMT NOXREF ALC NOANSF NOTER
ISN 0002
ISN 0003
ISN 0004
ISN 0005
ISN 0006
ISN 0007
ISN 0008
ISN 0009
ISN 0010
ISN 0011
ISN 0012
ISN 0013
ISN 0014
ISN 0015
ISN 0016
ISN 0017
ISN 0018
ISN 0019
INTEGER WINDY1,WINDY2,WINDX1,WINDX2 WIDTH1,WIDTH2,HI GHTIHIGHT2,
1 DELTAX,DELTAYCOUNT
INTEGER IMAGE1(2),WINDXWINDY.IMAGE2(2).SERIES( 800,63) ,FRAME(
13,400),NP(;3) FLAG(3,400),FLAG2.FLAG3,IMAGE3(2) FLAGI
C ORIGINAL FLOW CHART BY OR R E ELKINS
C ANALYSIS FOR APPLICATION BY GARY R JACKMAN
C PARTICLE FOLLOWING PROGRAMED BY DR J K YOO
C SERIES MATCHUP PROGRAM BY GARY R JACKMAN
C PARTICLE IMAGE IDENTIFICATION PROGRAM
C THIS PROGRAM FIRST IDENTIFIES A SERIES OF A PARTICLE IMAGES
C THEN MATCHES PAIRS OF SERIES THAT DEFINE A PARTICLES MOTION
C ***** **********************
C PARTICLE IDENTIFIER
C FLAGS =0:NOT TESTED OR NOT IDENTIFIED
C FLAGS =1:TESTEO IN FRAME AND IDENTIFICATION IN PROGRESS(FRAME 2,3)
C FLAGS =2:SELECTED AS A SERIES OF A FLOW IMAGE
C FLAGS =3:TESTED AND CONSIDERED AS NOISE
C SERIES(1):PARTICLE ID
C SERIES(2): INITIAL FRAME
C SERIES(3): NO* OF PARTICLE IDENTIFIED IN SERIES
C INITALIZATION OF FLAGS AND SERIES
DO 1000 KOOK=1,400
DO 1000 JAE=1,3
1000 FLAG(JAEKOOKJ=O
DO 1001 KOOK=163
00 1001 JAE=1, 800
1001 SERIES(JAE,KOOK)=0
C SET INITIAL FRAME NUMBER=1
IFRAME=1
C DEFINE WINDOW SIZES
WIDTH1=30
HIGHT1=16
WI DTH2=20
HIGHT2=16
C READ THREE FRAMES FROM CARDS, WITH FRAME=J
COUNT=1
DO 10 J=1,3
READ(5.2001) NP(J)
2001 FORMAT(I4)
NPTWO=NP(J )*2
FIGURE 4.4 Computer program for following image positions and matching image pairs
MAIN OS/360 FORTRAN H EXTENDED
ISN 0020 READ(5.200'2) (FRAME(J,K).K=1,NPTWO)
ISN 0021 2002 FORMAT(2014)
ISN 0022 10 CONTINUE
C REINITIATE FRAME#
ISN 0023 101 NPl=NP(1)
ISN 0024 NP2=NP(2)
ISN 0025 NP3=NP(3)
ISN 0026 9 FLAG=0O
ISN 0027 J=l
C SELECT AN UNIDENTIFIED IMAGE FROM FRAME1) TO BE SET AS AN IMAGE
ISN 0028 2 DO 20 L1=1.,NP1
ISN 0029 JLI=L1*2-1
ISN 0030 IF(=LAG(JL1).EQ.O) GO TO 1
ISN 0032 20 CONTINUE
C IF ALL THE PARTICLES IN FRAME 1 ARE CHECKED, GO TO 4
ISN 0033 21 IF(=LAG1.EQ.O) GO TO 4
ISN 0035 STOP
C AN UNIDENTIFIED PARTICLE IS FOUND, RECORD IT AS IMAGES
ISN 0036 1 =LASl=1
ISN 0037 FLAG(J*L1)=1
ISN 0038 SERIES(COWNT,1)=COUNT
ISN 0039 IMAGEl(1)=FRAME(Jv.JL1)
ISN 0040 LLI=L1*2
ISN 0041 IMAGEI(2)=FRAME{JLL1)
C SET WINDOW AREA
ISN 0042 11 WIN3X=IMAGE1(1)+WIDTH1
ISN 0043 WIND'YI=IMAGE (2)+HIGHT1/2
ISN 0044 WINDY2=IMAGE1(2)-HIGHTI/2
ISN 0045 111 FLAG2=0
ISN 0046 J=2
C ANY IMAGE IN FRAME(2) IN WINDOW AREA OF IMAGES
ISN 0047 DO 30 L2=1,NP2
ISN 0048 JL2=L2*2-1
ISN 0049 IF(FLAG(JL2).NE.0) GO TO 30
ISN 0051 LL2=L2*2
ISN 0052 IF(FRAME(JJL2).GEAGIMAGEI(1).AND.FRAME(JJL2).LE.WINDX.ANDFRAME(J
1.LL2).GE.WINDY2.AND.FRAME(JLL2.).LE.WINDY1) G3 TO 31
ISN 0054 30 CONTINUE
C IF NO PARTICLE' IN WINDOW1.SET IMAGEl AS NOISE AND FIND ANOTHER PARTICLE IN
C FRAME(1)
ISN 0055 32 IF(FLAG2.NE*0) STOP
C NO PARTICLE IN WINDOW
ISN 0057 FLAG(1 L1)=3
C RESET FLAGS FOR THE UNACCEPTED PARTICLES IN FRAME2 WINDOWl
ISN 0058 DO 321 K=1,NP2
ISN 0059 321 IF(FLAG(2,K).EQ.1) FLAG(2,K)=C
FIGURE 4.4 continued
LEVEL 2.1 ( JAN 75 )
DATE 76.251/15.
MAIN OS/360 FORTRAN H EXTENDED
ISN 0061
[SN 0062
ISN 0063
ISN 0064
ISN 0065
ISN 0066
ISN 0067
ISN 0068
ISN 0069
ISN .0070
ISN 0071
ISN 0072
ISN 0073
ISN 0074
ISN 0075
ISN 0076
ISN 0078
ISN 0079
ISN 0081
ISN 0082
ISN 0083
ISN 0085
ISN 0086
ISN 0087
ISN 0088
ISN 0089
ISN 0090
ISN 0091
ISN 0092
ISN 0093
ISN 0094
ISN 0095
GO TO 9
C YES, THERE FOUND A PARTICLE IN WINDOWl. RECORD IT AS
31 FLAG2=1
FLAG(JL2)=1
IMAGE2( 1)FRAME(J.JL2)
IMA6E2(2)=FRAME(J,-LL2)
IMAGE
C EXTRAPOLATE THE THIRD POSITION BY VECTOR MANIPULATION
C FROM IMAGES TB IMAGE2 INSIDE THE THIRD FRAME WINDOW AREA
C
C CALCULATE VECTOR
33 DELTAX=IMAGE2(1)-IMAGE1(1)
DELTAY=IMAGE2(2)-IMAGE1(2)
J=3
FLAG3=C
C SET WINDOW AREA
WINDXI=IMAGF2(:1)+DELTAX-WIDTH2/2
WIN3X2=IMAGE2(1)+DELTAX+WIDTH2/2
WIN3YI=IMAGE2(2)+ DOLTAY-HIGHT2/2
WIN)Y2=IMAGE2(2)+DELTAY+HIGHT2/2
C SEARCH FOR A POINT IN THE THIRD FRAME
DO 40 L3=I,NP3
JL3=L3*2-1
IF(FLAG(JL3).NE.O0 GO TO 40
LL3=L3*2
IF(FRAME(JJL3).GE.WINDX1.AND.FRAME(JJL3).LE.WINDX2.AND.FRAME(J*
1LL3).GE.WINDY1.AND.FRAME(J,LL3). LE.WINDY2) GO TO 41
GO TO 40
C YES, THERE IS A PARTICLE IN WINDOW
41 FLAG3=FLAG3+1
IF(=LAG3.GE.2) GO TO 44
C ASSUME THIS PARTICLE IS A PART OF A SERIES CURRENTLY SEARCHING
FLAG(J L3)=1
IMAGE3(1 )=FRAME(J.-JL3)
IMAGE3(2)=FRAME{JLL3)
C STORE INITIAL FRAME NUMBER AND NJMBER OF PARTICLES IDENTIFIED IN THIS
C SERIES AND IMAGES
42 SERIES(COUNT,2)=IFRAME
C STORE SERIES OF THREE IMAGES IN ARRAY SERIES(800,63)
SERIES(COUNT 3)=SERIESICOUNT,3)+3
SERIES(COUNT 4)=IMAGE(11)
SERIES(COUNTs5)=IMAGE1(2)
SERIES(COUNT,6)=IMAGE2(1)
SERIES(COUNT7 )=IMAGE2(2)
SERIES(COUNT,8)=IMAGE3(1)
SERIES(COUNT,9)=IMAGE3(2)
C SET FLAGS FOR THE IDENTIFIED SERIES
FIGURE 4.4 continued
LEVEL 2.1 ( JAN 75 )
DATE 76.251/15o
MAIN OS/360 FORTRAN H EXTENDED
ISN 0096
ISN 0097
ISN 0098
ISN 0099
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
0100
0102
0104
0105
0106
0107
0108
0110
0111
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
FLA6(3 L31=2
FLAG(2.L2)=2
FLAG( 1 L 1=2
40 CONTINUE
C IF THERE WAS NO MATCH IN FRAME3, GO TO FRAME TO FIND ANOTHER
C PARTICLE IN WINDOW
IF(cLAG3.EQ.0) GO TO 11
C IF THERE IS AN EXACTLY' ONE PARTICLE IN WINDOW, THAT IS A SERIES
IF(FLAG3.EQ.11 GO TO 422
STOP
C IF MORE THAN I IN WINDOW2, SET 1ST ONE AS SERIES AND WRITE MESSAGE
44 WRITE(6,441) IMAGE1
441 FORMAT('**** MORE THAN 1 PARTICLE IN WINDOW 2 *** FOR IMAGE 2('*2I
15,')')
C RESET FLAGS FOR THE UNACCEPTED PARTICLES IN FRAME WINDOWi
422 DO 421 K=1INP2
421 IF(FLAG(2,KJ).EQ.l FLAG(2,K)=0
C INCREASE PARTICLE IDENTIFICATION # IE PARTICLE NAME(COUNT).
COUNT=COUNT+1
IF(COUNT.GE. 800) GO TO 1003
C GO TO FIND ANOTHER PARTICLE IN FRAME 1
GO TO 9
C ADJUST FRAMES 2 AND 3 INTO 1 AND 2
4 NPNP=NP2*2
NP(1)=NP2
COUNT=COUNT-1
DO 50 K=1,NPNP
FRAME(1,K)=FRAME(2.K)
50 FLA( I,K)=FLAG(2,K)
NPNP=NP3*2
NP(2)=NP3
DO 51 K=1,NPNP
FRAME( 2,K)=FRAME (3 K)
51 FLAG(2,K)=FLAG(3.K)
DO 52 K=1,400
52 FLAG(3.K)=0
C NOW READ NEW FRAME3
READ(5,2001) NP(3)
NPTWO=NP(3)*2
READ(5 2002) (FRAME(3,K),K=1,NPTO )
C INCREASE INITIAL FRAME NO.
IFRAME=IFRAME+1
IFRAM=IFRAME+2
WRITE(6.59) IFRAM
59 FORMAT (FR.AMEE' 5. IS ENTERED')-
C NUMBER OF FRAMES TO BE WORKED
FIGURE 4.4 continued
LEVEL 2.1 ( JAN 75 )
DATE 76.251/15.
MAIN OS/360 FORTRAN H EXTENDED
ISN 0134 IF(IFRAM.EGO168) GO TO 1003
C EXTRAPOLATE A POINT BY VECTOR FROM THE LAST TWO POINTS.
ISN 0136 DO 50 K=1-COUNT
ISN 0137 IF((IFRAM-(SERIES(r(P3)+SERIES(K,2))).GE*I) GO TO 60
ISN 0139 IF(SERIES(K,3J*GE.30) GO TO 60
ISN 0141 M8=SERIES( K,3)*2+2
ISN 0142 M9=M8+1
ISN 0143 M7=M8-1
ISN 0144 M6=M8-2
ISN 0145 K8=M8
ISN 0146 K9=M9
ISN 0147 DELTAX=SER.IES(K M8 )-SERIES(K M6J
ISN 0148 DELTAY=SER.IES(KM9)-SERIES(K*M7)
ISN 0149 WINDXI=SER;IES(K.M8)+DELTAX-WIDTH2/2
ISN 0150 WINDX2=SEPIES(KMS)+DELTAX+WIDTH2/2
ISN 0151 WIN3Y1=SER;IES tK, Mg)+DELTAY-HI GHT2/2
ISN 0152 WINDY2=SERIES(KM9)+DELTAY+HIGHT2/2
ISN 0153 FLAG3=0
ISN 0154 NP3=NP(3)
ISN 0155 00 61 L3=I1NP3
ISN 0156 JL3=L3*2-1.
ISN 0157 IF(=LAG(3,L3).NE.0) GO TO 61
ISN 0159 LL3=L3*2
ISN 0160 IF(=RAME(3 JL3).GE.WINDXI.AND.FRAME(3eJL3).LE*WINDX2*AND.FRAME(3,
1LL3)JGE*WINDY1.AND.FRAME(3,LL3).LE*WIN3Y2) GO TO 62
ISN 0162 GO TO 61
C YES, A ARTICLE IN WINDOW 2
ISN 0163 62 FLAG3=FLAG3+1
ISN 0164 IF(.LAG3.GE*2) GO TO 63
C ASSUME THIS PARTICLE IS A PART OF SERIES
ISN 0166 FLAG(3.L3)=1
ISN 0167 M9=M9+2
ISN 0168 M8=M8+2
ISN 0169 IMAGE3(1)J FRAME(3,JL3)
ISN 0170 IMAGE3(2)=FRAME(3.LL3)
ISN 0171 64 SERIES(KM8)=IMAGE3(1)
ISN 0172 SERIES(K.M9)=IMAGE3(2)
ISN 0173 SERIES(K.3)=SERIES(K.3)+1
ISN 0174 FLAG(3.L3)=2
ISN 0175 61 CONTINUE
ISN 0176 IF(FLAG3.LE*1) GO TO 60
ISN 0178 STOP
ISN 0179 63 WRITE(5,441) SERIES(KK8),SERIES(KK9)
ISN 0180 60 CONTINUE
C NOW TAKE CARE OF PARTICLES NEWLY APPEARED IN FRAME 1 AT THIS TIME
ISN 0181 COUNT=COUNT+1
FIGURE 4.4 continued
DATE 76.251/15,
LEVEL 2.1e ( JAN 75 )
MAIN OS/360 FORTRAN H EXTENDED
ISN 0182 IF(COUNT.GE* 800) GO TO 1003
ISN 0184 GO TO 9
ISN 0185 1003 CONTINUE
C PRINT TEST OUTPUT
ISN 0186 WRITE(6,1006)
ISN 0187 1006 FORMAT('THE FOLLOWING SERIES ARE FOUND' )
ISN 0188 DO 1005 J1f,COUNT
ISN 0189 ITOTAL=SERIIES(J,3)*2+3
ISN 0190 1005 WRITE(6.1004) (SERIES(JK),K=1,ITOTAL);
ISN 0191 1004 FORMAT(/,3I6,/,10('(',215,'J"))
C ******** ************** ***** ** ********* **** ** t 4t
C TRY TO MATCH A PAIR OF IDENTIFIED SERIES
C ****************** *******#*#**'#***t***********#*(ii^*)te(ct *
ISN 0192 JJJ=C
ISN 0193 DO 700 M=1,COUNT
ISN 0194 ISM3=SERIES(4M,3)
ISN 0195 JFRAME=SER;IES(.M.2)
ISN 0196 IF(ISM3.LT.5) GO TO 700
ISN 0198 IF(SERIES(fM,5).GT.500) GO TO 700
C TRY TO COMPARE' TWO SERIES M,N
ISN 0200 DO 70 N= 1 COUNT
ISN 0201 IF(M.EQ.N) GO TO 70
ISN 0203 ISN3=SERIES(N,3)
ISN 0204 IF(ISN3.LT.5) GO TO 70
ISN 0206 IF(SERIES4N,5).LT.500) GO TO 70
ISN 0208 IFRAME=SERIES(N.2)
ISN 0209 IF(JFRAMEGE*IFRAME) GO TO 722
C INITIAL FRAME NO. IFRAME OF SECOND SERIES IS LARGER THAN JFRAME OF FIRSTSERIES
ISN 0211 IF((IFRAME-JFRAME).GT.2) GO TO 70
ISN 0213 JK=1
ISN 0214 MK=IFRAME-JFRAME+1
ISN 0215 JACK-=JK-1
ISN 0216 MEET=0
ISN 0217 MDE==ISM3
ISN 0218 MKL=MDEF-MK+1
ISN 0219 NASC=ISN3-MK+1
ISN 0220 IF(NASC*LT.MKL) MDiF=MK+NABC-1
ISN 0222 DO 18 MANTMK,MDEF
ISN 0223 MX=(MAN+1)*2
ISN 0224 JACK=JACK+
ISN 0225 JX=(JACK+I)*2
ISN 0226 KFRAME=IFRAME
ISN 0227 IF(IABS(SE-RIES(M,MX)-SERIES(NJX)).LE.6) GO TO 17
ISN 0229 GO TO 70
C FOUND ONE MATCH
ISN 0230 17 MEET=MEET+1
FIGURE 4.4 continued
LEVEL 2.1 ( JAN 75 )
DATE 76.251/15
MAIN OS/360 FORTRAN H EXTENDED
ISN 0231
ISN 0233
ISN 0234
ISN 0235
ISN 0237
ISN 0238
ISN 0239
ISN 0240
ISN 0241
ISN 0242
ISN 0243
ISN 0244
ISN 0246
ISN 0247
ISN 0248
ISN 0249
ISN 0250
ISN 0251
ISN 0253
ISN 0254
ISN 0255
ISN 0257
ISN 0258
ISN 0259
ISN 0260
ISN 0261
ISN 0262
ISN 0263
ISN 0264
ISN 0265
ISN 0266
ISN 0267
ISN 0268
ISN 0269
ISN 0270
ISN 0271
ISN 0272
ISN 0273
ISN 0274
IF(MEET*GE.15) GO TO 71
18 CONTINUE
GO TO 70
C INITIAL FRAME NO. OF 2ND SERIES (IFRAME) IS LESS THAN OR EQUAL TO
C JFRAME OF FIRST SERIES
722 IF((IJFRAME-IFRAME).GT.2) GO TO 70
MK=1
JK=JFRAME-IRRAME+1
726 JACK=JK-1
MEET=0
MDEF=ISM3
MKL=MDEF-MK+1
NABC=I SN3- JK+
IF(NABC*LT .MKL) MDEF=MK+NABC-1
00 723 MAN=MKMDEF
MX=(MAN+1 ) *2
JACK=JACK+1
JX=(JACK+1)*2
KFRAME=JFR.AME
IF(IABS(SERIES(MMX)-SERIES(NJX)).LE*6) GO TO 725
GO TO 70
C FOUND ONE MATCH
725 MEET=MEET+.l
IF(MEET.GE.15) GO TO 71
723 CONTINUE
GO TO 70
C FOUND A PAIR OF SERIES M.N
71 WRITE(6.711) SERIES(M.)11SERIES(NI1
711 FORMAT(' SERIES ',14, AND I*4,' ARE A PAIR')
C (PARTICLE ID, INITIAL FRAME), ((XY.,Z), (XYZ),. (0,00))
C PUNCH OUT THE OUTPUT IN POSITION INPUT FORMAT
JJJ=JJJ+1
WRITE(6.72) JJ-JKFRAME
72 FORMAT(215)
WRITE(7,75) JJJ.KFRAME
75 FORWAT{I4.2X.I4)
DO 73 LNG:MK,MDEF
LONG=(LNG+,1)*2
LONG1=LONG+l
MALL=( JK+1 )*2+1
MALL1=MALL-1
JJ=K+1
WRITE(7,76) SERIES(MLONG),SERIES(N.MALL),SERIES(MLONG1)
76 FORMAT( I4,13.13.2X)
73 WRITE(6,731) SERIES(M ,LONG).SERIES(M ,LONG1).SERIES(N .MALL)
1.SERIES(N *MALL1)
FIGURE 4.4 continued
DATE 76.251/154
LEVEL 2-.1 ( JAN 75 )
ISN 0275 731 FORMAT(1X,3I15,17)
ISN 0276 WRITE(7,77)
ISN 0277 77 FORMAT(14HO000000CO0000)
ISN 0278 GO TO 700
ISN 0279 70 CONTINUE
ISN 0280 700 CONTINUE
ISN 0281 701 STOP
ISN 0282 END
ln
FIGURE 4.4 continued
MAIN OS/360 FORTRAN H EXTENDED
DATE 76*251/15*
LEVEL 2*1 ( JAN 75 )
4.4 Method for Determination of Three Dimensional Positions and Velocities
The records from the method of following particle images and matching
image pairs contain the particle number, the initial frame number, and
the listing of coordinates for each image pair of the series for each
particle. This information is reduced to three dimensional coordinate
series or reconstructed particle paths. The data must be filtered to
remove digital noise obscuring the true image positions. The filtered
position data are used to calculate the velocity components at each data
point and the statistical quantities for the flow field. The results
of this method can be used to describe the instantaneous flow field as
well as the time and space averaged flow field. This method and its
computer program will be discussed in four sections. The first subject
addressed will be the conversion of the input series pairs to three di-
mensional coordinates. The second topic will be the digital filtering
out of noise induced by the digitization process. The calculation of
instantaneous particle velocities at the particle path data points will
be the third area of study. The final subject will be the method of
acquiring statistical information about the flow field.
The first step of this method is converting image pairs to three
dimensional cylindrical coordinates. The input data are of the form of
a series of pairs of two dimensional particle image coordinates referenced
to the prism. These pairs consists of three coordinate values: The
distance of the image pair along the horizontal from the prism edge,
labeled "IX" in Figure 4.5; the vertical distance of the image in the
upper prism face from the prism axis (two face intersection), labeled
"IZB" in Figure 4.5; and the vertical distance of the image in the lower
prism face from the same prism axis (two face intersection), labeled
"IZE" in Figure 4.5. The integers are converted to real numbers and
referenced to the prism axis for use in the computer program.
These values are processed through the prism equations (3.1) and (3.2)
to yield Cartesian coordinates which are converted into cylindrical
coordinates for use in the remaining parts of the program. The height
of the apex of the prism above the centerline of the test pipe ( "h"
in Figure A-2) must be changed when changing prisms for data acquisition.
The angle 4 referred to in the prism analysis (Appendix A) is the
angle "ANG", in degrees, of this computer program. The easiest method
of determining this angle is to locate two particle image paired series,
one for a particle in the lower section of the pipe and one for a
particle in the upper part, which are located on the pipe wall and
vary the angle until the program locates both on the pipe wall. Mis-
alignments of up to 50 are very difficult to detect and are nearly
impossible to measure accurately. At the completion of this section
of the program, the reconstructed particle paths are expressed in
cylindrical coordinates with noise included.
The position error or noise from the digitization process is random
in amplitude up to the coordinate maximum errors discussed in section 4.1
and 4.2 and in all frequencies through half the cinematographic recording
rate. The position errors at higher frequencies are eliminated from the
reconstructed particle paths by using a digital filter. First, the
particle path position data are transformed into Fourier coefficients in
the frequency domain. The transformation is accomplished by using a
Fourier series of sine functions and varying the amplitudes of their
coefficients until the resultant series exactly matches the reconstructed
particle path. This accomplished by using a Fast Fourier Transform
(FFT) for sampled particle paths. The frequency domain is a plot of
the coefficients against their corresponding frequencies. The coefficients
above a predetermined frequency are set equal to zero. The remaining
terms are transformed back into the time domain using an inverse FFT
to yield reconstructed particle paths with no frequencies present above
the predetermined frequency which is the cutoff frequency for the data.
An interesting observation is that data can be totally obscured by noise
but providing the noise is of a different frequency than the data, it
can be removed or digitally filtered out leaving only noise free data.
This procedure of transforming the reconstructed particle paths into the
frequency domain, filtering out the frequencies above the cutoff frequency,
and transforming back into the time domain is the type of digital filtering
that is used in this program.
This digital filtering is a simple procedure to eliminate the higher
frequencies of noise but before using it, the pitfalls must be known.
The sample time increments must be equal. This is automatically accom-
plished since the cinematographic recording rate is constant and, thus,
the time interval between frames (sample time increment) is constant.
The highest frequency in the data must be at or belowhalf the cinemato-
graphic recording rate (folding frequency) in order to eliminate aliasingg"
which occurs when high frequency components of a time function match low
frequencies if the sampling rate is too low. The highest frequency noise
that can be induced in frame by frame data entry is at the folding fre-
quency so noise cannot be "aliased". The cinematographic recording rate
must be high enough to insure that the highest frequency of the turbulence
is below the folding frequency to eliminate turbulence frequency aliasingg".
Thus, sampling must be at a rate at least twice as high as the highest
frequency present in the turbulent fluctuations. This digital filter is equi-
valent to multiplying the Fourier coefficients by a rectangular function
of unit amplitude up to the cutoff frequency. This rectangular function
multiplication in the frequency domain is equivalent to performing a
convolution with the function (sin t)/t in the time domain. This function
does not go to zero but has side lobes which cause oscillatory variations
on the ends of the filtered particle path reconstruction. The oscilla-
tions are a small percentage of the amplitude of each end of the recon-
structed particle path that is processed through the FFT program and
rectangular filter. This truncation error was reduced to a minimum for
the rectangular filter by fitting a line through the end points of each
sampled particle path, and processing only deviations from that line in
order to reduce the amplitude of the time domain data. A plot of selected
reconstructed particle paths with and without filtering shows that phase
shift is not a problem. When the filtered frequencies are the higher
frequencies in the spectrum, the filtered particle path reconstruction
appears smoothed much the same as if it were subjected to a polynomial
smoothing formula.
In the second section of the computer program the higher frequency
noise is removed from the data using a digital filter in conjunction with
the Fast Fourier Transform (FFT) program written by Robinson, 1967. The
digital filter is entered by calling the subroutine "FILTER" in line 84
of the main program in Figure 4.5. The data areprocessed one reconstructed
particle path at a time with the filter using the cylindrical coordinates
of the reconstructed particle path data points ( x, rad, theta), the
number of data points in the particle path (JTOT), and the cinematographic
recording rate (RATE) as its working data. Before the data are trans-
formed to the frequency domain for filtering, the linear trend in-
formation is removed by constructing a straight line between the un-
filtered particle path end points and processing only the deviations
from this line through the FFT program and the digital filter. By
eliminating the trend information from consideration in this manner,
the amplitudes of the end points of the modified particle path data
processed through the filter are zero which reduces the filter error
at the ends of the filtered particle path reconstructions that results
from truncation of the Fourier series when using a square filter. The
FFT program is entered at this point with the trend removed data by
calling the subroutine "NLOGN" for each of the three coordinates in
lines 41, 42, and 43 of the subroutine "FILTER". The transformed
data are returned to the subroutine "FILTER". The highest retained
frequency is set by altering line 44 of this subroutine. For this
work, 6 Hz was the maximum retained frequency (see section 5.1). The
data, now in the frequency domain, have the amplitudes of all frequencies
above 6 Hz set equal to zero which eliminates the higher frequencies
in the data. The data are transformed back to the time domain by
calling the subroutine "NLOGN" in lines 60, 61, and 62 which are in-
verse transforms. The digitally filtered fluctuating data are recom-
bined with the trend information and returned to the main program. The
now filtered particle path data can be used with a much higher degree of
confidence since the noise above 6 Hz has been eliminated. If tempo-
rarily.by-passing the digital filter becomes necessary, it can be ac-
complished by replacing "go to 117" with "go to 2000" in line 35 of
the main program.
62
The velocities of the particles can be calculated at each data
point in the reconstructed particle paths. The three components of
velocity are calculated for both the cylindrical and Cartesian coordi-
nate systems since both will be utilized in later programs. A three
data point calculation of velocity is used since it approximates the
reconstructed particle path by a parabola from time tl to time t3. This
method approximates a parabola only if the time increments are equal as
is the case in this work. The velocity components of point (i) are
calculated by using the equations,
( +1- x. )
Si+1l 1-1
i = ( t t i)
Ri i(i+1 e1)
U
i (t i+ t )
where Uxi is the instantaneous axial velocity, URi is the instantaneous
radial velocity, U i is the instantaneous circumferential velocity, and
(Xi, R 0i, ti) are the cylindrical coordinates and time at reconstructed
particle position (i). With the positions and velocities determined at
the data points, the average velocities and root-mean-square velocities
can be calculated.
In the final section of the computer program, the average and root-
mean-square velocities are calculated for the time and space-fixed volume
observed in this work. These quantities are determined at radial incre-
ments. The assumption is made that these turbulent flow quantities only
have a radial dependency. The flow volume viewed through the prism is
63
divided into shells which have a radial thickness of 2% of the radius.
One of the shells is shown in Figure 3.2. The velocities are summed for
all reconstructed particle positions in each shell for a spatial total,
and over time fora space-time total which is then divided by the number
of data points in the space-time shell to yield space-time average
velocities. An identical summing procedure is used for the squared
velocities which are needed to calculate the turbulent intensities or
root-mean-square velocities. The average velocities and turbulent in-
tensities are computed using the equations,
K
z Uxi
S i=l
Ux =
K
E U R"
i=l
R K
K
E U0i
i=l
U =K
K
K U2
E U-
U =
UXRMS K U
K 2
i=lRi
RRMS K
K
/ K 2
E U
i=l ei
ORMS K
where ( UX, UR, U ) are the cylindrical components of the average velocity
for a shell, (UXRMS, URRMS, UeRMS) are the turbulent intensities or root-
mean-square velocities for a shell, and K is the total number of data
points utilized in the shell calculation.
The output from this program is average and root-mean-square
velocities, and instantaneous reconstructed particle path positions
and velocities. The instantaneous positions and velocities are recorded
on computer print-out and cards at the reconstructed particle path
data points. The information is needed in this order for use in con-
structing a visual display of the reconstructed particle paths, but
particle path data positions must be grouped by frame for use in con-
structing visual displays of instantaneous velocity fields and referencing
the data to a fixed grid into which the studied flow volume may be
referenced. A copy of the computer cards arranged by frame is made
using a simple computer program. This program is listed in Figure 4.6
for continuity of the method of data reduction.
224-FOLLOW -FORT -SYSPRINT
LEVEL 2.1 ( JAN 75 )
OS/360 FORTRAN H EXTENDED
DATE 76.251/15.
REQUESTED OPTIONS: NODECKNOLOAD.OPT=0,LC(48)
OPTIONS IN EFFECT: NAME(MAIN) NOOPTIMIZE LINECOUNT(484 SIZELMAX) AUTODBL(NONE)
SOURCE EBCDIC NOLIST NODECK NOOBJECT MAP NOFORMAT GOSTMT NOXREF ALC NOANSF NOTER
ISN 0002
ISN 0003
ISN 0004
ISN 0005
ISN 0006
ISN 0007
ISN 0008
ISN 0009
ISN 0010
ISN 0011
ISN 0012
ISN 0013
ISN 0014
ISN 0015
ISN 0016
ISN 0017
ISN 0018
ISN 0019
ISN 0020
ISN 0021
ISN 0022
ISN 0023
ISN 0024
ISN 0025
ISN 0026
ISN 0027
ISN 0029
!SN 0030
C THIS PROGRAM FINDS SMOOTHED POSITION AND VELOCITY WITH AVERAGE
C VELOCITIES AND RMS VALUES BY SHELLS
C READ. DATA INTO COMPUTER AND SET UP HEADINGS FOR OUTPUT
REAL NA.NG*ZB(64),ZE(64),X(64),Z(64),Y(64) RAD(64),THETA(64),
IVELX(64),VELR(64).VELT(64), NW
DIMENSION IX (64) ,IZB(64).IZE(64),VELY(64),VELZ(64) D(64).E(64)
IF(64)
COMPLEX A(*64).B(64),C(64)
DIMENSION SUMU(28)-SUMV(28),SUM (28 )SUMUU(28),SUMVV(28).SUMWW(28)
1,RATRO(28).K1(28)AUAVG(28).VAVG428).WAVG(28),URMS(28).VRMS(28),
IWRMS(28) RATUUA(29).RATVUA(29),RATWUA(29) RURMS(29 )RVRMS(29),
IRWRMS(29)
C SET INITIAL VALUES TO ZERO
DO 1-0 1=1,57
SUMU(I)=0
SUMV(I)=0
SUMW(I)=0
SUMUU(I)=0
SUMVV( I)=0
SUMWW(I)=0
K1(I)=0
10 CONTINUE
SUO=O
SVO=0
SWO=0
TURMS=0.
TVRMS=0.
TWRMS=0.
READ(5,100) JRE,RATEsSCALEHoSCALEVTEMP.IDATEXZEROZBZEROZEZERO
100 FORMAT( 14.2X.F4.142X.F6.4.2XF6.4.2XF4.1.2X.I6.3(2XF5.1)3
WRITE(6,101) JRE, TEMP. IDATE
101 FORMAT(1H1. 1OX,16HREYNOLDS NUMBER=.14,2HAT,F4. 12HON 16,/////)
102 READ(5,103) IPART,.J
103 FORMAT(14.2XI4)
IF(IPART.EQ.0)GO TO 109
C JRE=REYNOLDS NUMBER. RATE=FRAMES PER SECOND, SCALE=MAGNIFICATION
C RATIO OF FILM PRISM TO TRUE PRISMs IPART=PARTICLE, ZBZERO ZEZERO
C XZERO ARE INITIAL VALUES OF ZB ZE X,- J=FRAMES
WR.ITE(6,104)
104 FORMAT(9X.5HPRTCL,2X 1HI 2X.5HFRAME 8X IHY,8X1 HZ.8XIHX.8X 1HR,
I8X,5HTHETA.8SX.2HUX,8X,2HVR.8X.2HWT,8X,2HVY,8X.2HVZ,///)
FIGURE 4.5 Computer program to calculate filtered positions and velocities
MAIN OS/360 =ORTRA4 H EXTENDED
ISN 0031 JTOT=0
ISN 0032 00 108 1=1,65
ISN 0033 READ(5,105) IX(I ),IZB(I),IE( I )
ISN 0034 105 FORMAT ( I4, 13 I3,2X)
ISN 0035 IF(IX(I).EQ.O.AND*IZB(I).EQ.O.AND*IZE(1).EQ.0.) GO TO 117
ISN 0037 JTOT=JTOT+1
ISN 0038 IZB(I)=IZB(I)-500
ISN 0039 IZEII)=500-IZE(I)
ISN 0040 X(I=IX(I)
ISN 0041 ZB(I)=IZB l)
ISN 0042 ZEI )=IZE(-I)
C CALCULATION OF TRUE PARTICLE POSITION
ISN 0043 R=63.5
ISN 0044 NA=I.OO
ISN 0045 NG=1.50
ISN 0046 NW=1.33
ISN 0047 H=110.2
ISN 0048 PI=0.78540
ISN 0049 ANG=O.0
ISN 0050 ZB(I)=(Z8(aI)-Z8ZERO)*SCALEV
ISN 0051 ZE(I)=(ZE(I)-ZEZERO)*SCALEV
ISN 0052 X(I)=(X(I)-XZERO)*SCALEH
ISN 0053 DELTAC=ARSIN(NA/NG*SINIPI+ANG*PI/45 ))
ISN 0054 DELTAO=ARSIN(NA/NG*SIN(PI-ANG*PI/45 ).)
ISN 0055 TC=TAN(PI+OELTAC)
ISN 0056 TD=TAN(PI*DELTAD)
ISN 0057 U=TAN(PI)
ISN 0058 AA-TC**21 e
ISN 0059 BB=2**TC*(H-ZB(I)*U)-2**ZB(I*TC**2
ISN 0060 CC=(H-ZB(I)*'J**2-R**2-2.*ZB(I) TC*(H-ZB(I)*U)+Z( I)**2*TC**2
ISN 0061 ZC=('-BB+SQRT(BB*BB-*4.*AA*CC))/(2.AA)
ISN 0062 OD=TD**2+1.
ISN 0063 EE=-2**TD*(H-ZE(I)*U)+2.*ZE(I)*TD**2
ISN 0064 FF=(H-ZE(I)*U)**2-R**2-2 *ZE(I )TD*(H-ZE(I)*U+ZE( )**2*TD#2
ISN 0065 ZD=(-EE-SQRT(EE*EE-4.*DD*FF))/(2**DD)
ISN 0066 VC=SQRT(R*R-ZC*ZC)
ISN 0067 VD=SQRT(R*R-ZD*ZDJ
ISN 0068 THETAC=ATAN2(ZC*VC)
ISN 0069 THETAD=ATAN2(-ZD VD)
ISN 0070 GAMMAC=ARSIN(NG/NSSIN(PI-DELTTAC-THETAC))
ISN 0071 GAMMAD=ARS[IN(NG/NIwSIN(PI-DELTAD-THETAD )
ISN 0072 WC=TAN(2.*PI-THETAC-GAMMAC)
ISN 0073 WD=TAN(2.*PI-THETAD-GAMMAD)
ISN 0074 Z( I=(ZC*WC+ZD*W0-VC+VD)/{WC+WD)
ISN 0075 Y(I)=VC+(Z(I)-ZC)*WC
C CONVERT FROM CARTESIAN TO CYLINDRICAL
FIGURE 4.5 continued
LEVEL 2.1 ( JAN 75 )
DATE 76.251/15.
OS/360 FORTRAN H EXTENDED
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
0076
0077
0079
0080
0081
OC82
0083
0084
0085
0086
0087
0088
0089
0090
0091
0093
0094
0095
0096
0097
0098
0100
0102
0103
0104
0105
0106
0107
0108
0110
0111
0112
0113
ISN 0114
ISN 0116
RAD( I )=SQRT( Y4 I )*Y(I)+Z(I)*Z(I)I
IF(Z(I).EQ.0.0.AND.Y(I).EQ.0.0) GO TO 113
THETA(I)=(ATAN2(Z(I),Y(I)))*45./PI
GO TO 115
113 THETA(I)=0.0
115 CONTINUE
108 CONTINUE
C FILTER NOISE OUT
117 CALL FILTER(X,RADTHETA.JTOT.RATE)
2000 CONTINUE
C CONVERT FROM CYLINDRICAL TO CARTESIAN
DO 112 I=1.JTOT
Y(IJ=RAD(I)*COS(THETA(I)*PI/45.3
Z(I)=RAD(I)*SIN(THETA(I)*PI/45.i
112 CONTINUE
DO 118 I=1IJTOT
C CALCULATION OF VELOCITIES
IF(I.LE.2.)GO TO 114
VELX(I-1)=(X(IJ-X(I-2))*RATE/2.
VELR(I-1)=(RAD(I)-RAD(I-2))*RATE/2.
VELT(I-1)=(THETALI)-THETA(1-2))*RAD(I-1)*PI/45.*RATE/2.
VELY(I-1)=(Y(I)-Y(II-2))*RATE/2.
VELZ(I-1)=(ZIZ(-Z(I-2))*RATE/2.
C THIS WILL CORRECT VELT ACROSS THETA= 180.0
IF((THETA( I)-THETA(I-2)).GT.180.0) VELT(I-1)=-VELT(I-1)+
1(6.28319*RAD(I-1)*RATE/2.,
IF((:THETA(
1(6.28319*R;AD I-1J)RATE/2. )
114 CONTINUE
C VELOCITIES ON FIRST POSITION SET EQUAL TO ZERO BUT DO NOT USE
VEL(1I)=0.
VELR(1)-=0
VELT(1)=0.
VELY(1)=0b
VELZ 1)=0.
IF(I.EQ.1)GO TO 107
K=I-1
L=J-1
WRITE(6,106) IPART K, L, Y(K), Z(K), X(K). RAD(K). THETA(K),
IVELK(K), VELR(K), VELT(K),VELY(K), VELZ(K)
106 FORMAT(113, 14l7,F12.2,F9.2,F9.2.F9.2,F11.2,5(FI0.2))
C THIS SECTION WILL YIELD AVERAGE VELOCITIES AND RMS VELOCITIES
C RATIO 0= RADIUS OF DATA POINT TO PIPE RADIUS TO YIELD FRACTION OF
C RADIUS 3F PIPE THE DATA POINT IS IN FROM THE WALL
IF(I.EQ.2) GO TO 107
RATIOR=1.-RAD(K)/63.5
FIGURE 4.5 continued
MAIN
DATE 76.251/15.
LEVEL 2.1 ( JAN 75 )
MAIN 05/360 FORTRAN H EXTENDED
ISN 0117
ISN 0118
ISN 0119
ISN 0121
ISN 0122
ISN 0123
ISN 0124
ISN 0125
ISN 0126
ISN 0127
ISN 0128
ISN 0129
ISN 0130
ISN 0131
ISN 0132
ISN 0133
ISN 0134
ISN 0135
ISN 0136
ISN 0137
ISN 0138
ISN 0139
ISN 0140
ISN 0141
ISN 0142
ISN 0143
ISN 0144
ISN 0145
ISN 0146
ISN 0147
ISN 0148
ISN 0149
ISN 0150
ISN 0151
ISN 0152
ISN 0153
ISN 0154
ISN 0155
ISN 0156
ISN 0157
ISN 0158
DO 20 M=1,28
C ESTABLISH SHELL CENTERS AT 2% INCREMENTS AND CHECK DATA POINT IN SHELL
RAT:O(M)=(M-1)*0.02+0.010
IF(RATIOR.LT.,(RATRO(M)-0.010). OR.RATIORGT.(RATRO(M)+0.010))GOT020
SUMU(M)=SUMU(M)+VELX(K)
SUMV(M)=SUMV(M)+VELR(K)
SUMW(M)=SUMW(M)+VELT(K)
SUMUU(M)=SUMUU(r)+VELX(K)*VELX(K)
SUMvV(M)=SUMVV(M).VELR(K) *VELR(K)
SUMWW(M)=SUMWW( M)+VELT(K)*VELT(K)
K1(M)=K1(M)+1
20 CONTINUE
107 J=J 1
118 CONTINUE
110 WRITE(6,111)
111 FORMAT(/)
GO TO 102
109 DO 30 J=1,28
C AVERAGE VELOCITIES AT THE SHELL LOCATIONS
IF(Kl(J)) 11,11,12
11 UAVS(J)=0
VAVG(J)=0
WAVG(J)=0
URMS(J)=0
VRMS(J)=0
WRMS(J)=0
GO TO 13
12 UAVG(J)=SUMU(J)/K1(J)
VAVG(J)=SUMV(J)/K (J)
WAVG( J )=SUMW(J)/K (J)
C RMS VELOCITIES AT SHELL LOCATIONS
URMS(J)=SQRT(SUMUU(J)/Kl(J)-UAVG(J2*UAVG(J))
VRMS(J)=SQRT(SUMVV(J)/K1(J))
WRMS(J)=SQRT(SUMWW(J)/K1(J))
C
13 SUO=SUO+UAVG(J)
SVO=SVO+VAVG(JJ
SWO=SWO+WAVG(J)
TURMS=TURMS+URMS(J)
TVRMS=TVRMS+VRMS(J)
TWRMS=TwRMS+WRMS(J)
C UA IS PIPE AVERAGE VELOCITY IN MM/SEC
UA=29.1
RATUUA(J)=UAVG(J)/UA
RATVUA(J)=VAVG(J)/UA
RATWUA(J)=WAV( J)/UA
FIGURE 4.5 continued
LEVEL 2.1 ( JAN 75 )
DATE 76.251/15
MAIN 0S/360 FORTRAN H EXTENDED
ISN 0159 RURMS(J)=URMS(J)/UA
ISN 0160 RVRMS(J)=VRMS(JJ/UA
ISN 0161 RWRMS(J)=WRMS(J)/UA
ISN 0162 30 CONTINUE
C CALCULATE AREA AVERAGE VELOCITIES (UCoVOWO)-28.IS NJ. OF SHELLS
ISN 0163 UO=SUO/28.
ISN 0164 VO=SIVO/28.
ISN 0165 WO=SWO/28.
C CALCULATE AREA AVERAGE RMS VALUES
ISN 0166 AURMS=TURMS/28.
ISN 0167 AVRMS=TVRMS/28.
ISN 0168 AWRMS=TWRMS/28.
C
ISN 0169 WRITE(6,200) UOVOWO
ISN 0170 200 FORMAT(1H1.1OX,*35HTHE AREA AVERAGE VELOCITIES ARE UO=,F6.2.4X,3HVO
1=,F6.2,4Xv3HWO=,F6.2,//)
ISN 0171 WRITE(6.201)
ISN 0172 201 FORMAT(5X,3HY/R 5X,5HUXAVG,5495HVRAVG,5X.5HWTAVGe6X,5HUXRMS X, 5HV
1RRIMS,5X5HWTRMS10OX.3HNO../)
ISN 0173 WRITE(6,202) (RATRO(J),UAVG(J) VAVG(J)*WAVG(J) URMS(J) VRMS(J).
1WRMS(J)(KI(J),RATUUA(J),RATVUAJ)RATAJRJRM(JRVRMS(J
1RWRMS(J) J=1,28)
ISN 0174 202 FORMAT(4X,F5*3,4XF5.2.SX,F5.2,5XF5.2,5XF6.2,4XF6.2,4X.F6.28X,
1I4,4X,6F8.-3)
ISN 0175 WRITE(6,203) AURMSAVRMSAWRMS
ISN 0176 203 FORMAT(////,10X,38HTHE AREA AVERAGE RMS VALUES ARE AURMS=,F6.2,4X.
16HAVRMS=,F6.2*4X,&HAWRMS=.F6.2)
ISN 0177 STOP
ISN 0178 END
FIGURE 4.5 continued
DATE 76*251/15.
LEVEL 2.1 ( JAN 75 )
REQUESTED OPTIONS: NODECK,NOLOAD:rOPT=0,LC(48)
OPTIONS IN EFFECT: NAME(MAIN) NOOPTIMIZE LINECOUNT(48) SIZE(MAX) AUTODBL(NONE)
SOURCE EBCDIC NOLIST NODECK NOOBJECT MAP NOFORMAT GOSTMT NOXREF ALC NOANSF NOTE
ISN 0002 SUBROUTINE FILTER(XRTHETA.JRATE)
ISN 0003 DIMENSION X(64),R(64),TrIETA(64).D(64)JE(64).F(64)
ISN 0004 COMPLEX A(64),.(6 (6C(64),WKHOLD*Q
C THIS SECTION FILTERS HIGHER FREQUENCY NOISE OUT OF DATA USING FAST-
C FOURIER TRANSFORM METHOD
ISN 0005 3 IF(JULE.2) N=1
ISN 0007 IF(J.GT.2.AND.J.LE.4) N=2
ISN 0009 IF(J.GT.4.ANDOJ LE .8) N=3
ISN 0011 IF(U.GT*8.ANDO.JLE.16) N=4
ISN 0013 IF(J.GT.16.AND.J.LE.32) N=5
ISN 0015 IF(J*GT.32.ANDOJ.LE.64) N=6
ISN 0017 K=J+1
ISN 0018 L=2**N
C USE ENDPOINTS OF SERIES DATA FOR STRAIGHT LINE TO ELIMINATE SIGNAL
C OF GROSS MOVEMENT OF PARTICLE FROM FOURIER CONSIDERATION
ISN 0019 XSLOPE=(X(J)-X(1))/(J-1)
ISN 0020 RSLOPE=(R(J)-R(1))/(J-1)
ISN 0021 THETAS=(THETA{J)-THETA(1))/(J-1)
ISN 0022 XBEG=X(1)
ISN 0023 RBEG=R(1)
ISN 0024 THETAB=THETA(1)
ISN 0025 DO 5 I=1,J
ISN 0026 XO=XSLOPE(1I-1 )+XBEG
ISN 0027 RO=RSLOPE*(I-I1+RBEG
ISN 0028 THETAO=THETAS*(I-1)+THETAB
ISN 0029 A(IU=X(I)-XO
ISN 0030 B(I)=R(I)-RO
ISN 0031 C(I)=THETA(I)-THETAC
ISN 0032 5 CONTINUE
ISN 0033 DO 1 I=KL
ISN 0034 X(I)=0.
ISN 0035 R(I)=0.
ISN 0036 THETA(I)=0.
ISN 0037 A(I)=0.
ISN 0033 B(I)=0.
ISN 0039 C(IJ=0.
ISN 0040 1 CONTINUE
ISN 0041 CALL NLOGN(N,A,1.0)
ISN 0042 CALL NLOGN(N,8B1.0) 0
ISN 0043 CALL NLOGN(NC,1.0)
C FILTER OUT FREQUENCIES ABOVE 'FREQ'
FIGURE 4.5 continued
S0/360 PORTRAY H EXTENDED
DATE 76.251/15
LEVEL 2.1 ( JAN 75 )
LEVEL 2.1
ISN 0044
ISN 0045
ISN 0046
ISN 0047
ISN 0048
ISN 0049
ISN 0051
ISN 0052
ISN 0053
ISN 0054
ISN 0055
ISN 0056
ISN 0057
ISN 0058
ISN 0059
ISN 0060
ISN 0061
ISN 0062
ISN 0063
ISN 0064
ISN 0065
ISN 0066
ISN 0067
ISN 0068
ISN 0069
ISN 0070
ISN 0071
ISN 0072
ISN 0073
ISN 0074
ISN 0075
ISN 0076
ISN 0077
OS/360 =ORTRAN H EXTENDED
FREQ=6.
IDELTA=L/2*(1.-FREQ*2./RATE)
LLIM=L/2+1-IDELTA
ULIM=L/2+1+IDELTA
DO 4 I=1lL
IF(I.LE.LLIM.OR.I.GE.ULIM) GO TO 4
A( I)=0.
B( I=0.
C( I=0.
4 CONTINUE
DO 5. I=1 L
D(I)=CABS(A(I))**2
E(I)=CABS(B( I))**2
F(I)=CABS(C( 1))**2
6 CONTINUE
CALL NLOGN(NA,-1.0)
CALL NLOGN(NB,-1.0)
CALL NLOGN(NC-1..0)
DO 7 I=1 J
XO=KSLOPE (I-1)+XBEG
RO=RSLOPE*(I-1)+RBEG
THETAO=THETAS*(I-1)+THETA3
X( I)=A(Ii+.XO
R(I)=B(I)+RO
THETA(I)=C(I)+THETAO
7 CONTINUE
DO 8 I=K,L
X(I)=A(I)
R(I)=B(I)
THETA(I)=C(<)
8 CONTINUE
RETURN
END
FIGURE 4.5 continued
DATE 76.251/15
( JAN 75 )
FILTER
REQUESTED OPTIONS: NODECK*NOLOAD.OPT=O,LC(48)
OPTIONS IN EFFECT: NAME(MAIN) NOOPTIMIZE LINECOUNT(481 SIZE(MAX) AUTODBL(NONE)
SOURCE EBCDIC NOLIST NODECK NOOBJECT MAP MOFORMAT GOSTMT NOXREF ALC NOANSF NOTE
ISN 0002
ISN 0003
C
ISN 0004
ISN 0005
ISN 0006
ISN 0007
ISN 0008
ISN 0009
ISN 0010
ISN 0011
ISN 0012
ISN 0013
ISN 0014
ISN 0015
ISN 0016
ISN 0017
ISN 0018
ISN 0019
ISN 0020
ISN 0021
ISN 0022
ISN 0023
ISN 0024
ISN 0025
ISN 0026
ISN 0027
ISN 0028
ISN 0029
ISN 0031
ISN 0032
ISN 0033
ISN 0034
ISN 0035
ISN 0037
ISN 0038
ISN 0039
FIGURE 4.5 continued
SUBROUTINE NLOGN(NX.SIGN)
FFT ROUTINE BY ROBINSON-1967
.NMAX=LARGEST VALUE OF N TO BE PROCESSED
NONDUMMY DIMENSION M(NMAX.)
FOR EXAMPLE, IF NMAX=6 THEN
DIMENSION M(6)
DIMENSION X(2**N)
DIMENSION X(64)
COMPLEX X,WKHOLD,Q
LX=2**N
DO 1 I=1 N
1 M(I)=2**(N-I)
DO 4 L=1,N
NBLOCK=2**(L-1)
LBLOCK=LX/NBLOCK
LBHALF=L6LOCK/2
K=0
DO 4 IBLOCK=1,NBLOCK
FK=KI
FLX=LX
V=SIGN*6.2831853*FK/FLX
WK=CMPLX(COS(V),SIN(V))
ISTART=LBLOCK*(IBLOCK-1)
DO 2 I=1,LBHALF
J=ISTART+I
JH=J+LBHALF
Q=X(JH)*WK
X( J )=X(J)-Q
X(J)=X(J)+Q
2 CONTINUE
DO 3 I=2,N
II=1
IF(K.LT.M(I)) GO TO 4
3 K=K-M(I)
4 K=K+M(II)
K=0
DO 7 J=l,LX
IF(K.LT.J) GC TO 5
HOL)=X(J)
X( J=X(K+1)
X(K+1)=HOLD
OS/360 =ORTRAN H EXTENDED
DATE 76.251/15
LEVEL 2.1 ( JAN 75 )
LEVEL 2.1 ( JAN 75 )
ISN 0040 5
ISN 0041
ISN 0042
ISN 0044 6
ISN 0045 7
ISN 0046
ISN 0048
ISN 0049 8
ISN 0050
ISN 0051
OS/350 FORTRAN H EXTENDED
DATE 76.251/15.
DO 6 I=1.N
II=I
IF(K.LT*M(I)) GO TO 7
K=K-M( I )
K=K+M(II)
IF(SISN.LT.O.0) RETURN
DO 8 I=1,LX
X( I=X(I)/FLX
RETURN
END
FIGURE 4.5 continued
NLOGN
641-FOLLOW -FORT -SYSPRINT
LEVEL 2.1 ( JAN 75 ) CS/360 FORTRAN H EXTENDED DATE 76.239/20.
FECUESTED OPTIONS: NCDECK,NOLOADCPT=0
OPTIONS IN EFFECT: NAME(AIh) NCCPTIMIZE LINECOLNT{74) SIZE(MAX) ALTODBL(NONE)
SOURCE EBCDIC NOLIST NODECK NOOBJECT MAP NOFORMAT GGSTMT NCXREF ALC NOANSF NOTER
C PROGRAM TO PUNCH CARDS EY FRAME
C MUST BE "COOC0" CARD AT END OF DATA INPUT
C INPUT CARDS ARRANGED BY PARTICLE, OUTPUT CARDS ARRANGED 2Y FRAME
C I IS ATLEAST NO. OF CARDS INPUT, N IS ACTUAL NC. OF CARDS INPUT
C IFRAME IS NO. OF FRAMES USED IN ANALYSIS
ISN 0002 DIMENSION DATA(4000.20),IDATA(4000,3)
ISN 0003 N=0
ISN 0004 DO 1 I=1,4000
ISN 0005 READ(5.116) (IDATA(IJ),J=13),(DATA(I.K),K=1,2C)
ISN 0006 116 FORMAT(14,I2,I2,20A3)
ISN 0007 IF(IDATA(I,3).EQ.C) GO TO 2
ISN C009 N=N+1
ISN 0010 1 CONTINUE
ISN 0011 2 DO 10 IFRAME=1,55
ISN 0012 WRITE(6,101) IFRAME
ISN 0013 101 FORMAT(1HI,//.30X.6HFRAME=.14,//)
ISN C014 WIRITE(6,102)
ISN C015 102 FORMAT(9X.5HPRTCL,2X 1HI,2X,5HFRAME98X,1HY,8X,1HZ,eX,1FX,8X,1HR,
18X SHTHETA,8X.2HUX,8X.2HVR,8X,2HWT 8X 2HVY X,2HVZ ///)
ISN 0016 DO 20 I=1,N
ISN 0017 IF(ICATA(I,3).NE.IFRAME) GO TO 20
ISN 0019 IF(IDATA(I,2).EC.1) GO TC 20
ISN 0021 NW ITE(5, 1C') (I ATA(I J),J=1,3), (CATA(I,K),K<=1,20)
ISN 0022 106 FORMAT(I13,14,I7,6X,2A3,3X,2A3,3X,2A3,3X,2A3*5X,2A3,5(4X,2A3))
ISN C023 WRITE(7.116) (ICATA(I1J),J=1,3),(DATA(IIK),K=1.20)
ISN 0024 20 CONTINUE
ISN 0025 10 CONTINUE
ISN 0025 STOP
ISN 0027 END
F-e
4N
FIGURE 4.6 Computer program to rearrange data by frame
S4.5 Referencing Data to a Grid
The three dimensional positions and velocities computed in the
previous section can be used for analysis and display or can be processed
through a position referencing procedure to yield velocities at regular
intervals in the flow field. The computer program to reference the
velocity field to specified fixed geometric grid locations in the flow
field is listed in Figure 4.7. Plane views can be constructed since the
reference locations lie on regularly spaced planes. The number of re-
ference or grid locations can be altered for clarity of display with
relatively little effort. Nonuniform data locations that change from
frame to frame must be studied much more closely to understand the
changing velocity field. Caution should be exercised when using any
system which interprets velocities at random locations to predict
velocities at grid locations since it necessarily smooths the velocity
field. This type of interpretive system should only be used in con-
junction with nonreferenced data displays in order to see what effect
the procedure has on the appearance of the velocity field. Any grid
referencing system must necessarily apply some kind of a weighted average
of the velocities over a predetermined volume of influence. This has
the tendency to smooth the velocities, suppressing the velocity gradients
and shear stresses. If large strain variations over short distances are
present in the flow field, these variations will also be reduced. The
result can be that a highly irregular velocity field may appear as a
relatively regular and smooth velocity field.
The referencing grid boundaries for this work range from 63.5 to
35.5 mm in the radial direction, 280 (+ 140) in the circumferential
direction, and 10 to 38 mm in the axial direction. The grid reference
positions are constructed at the intersections between radial planes
every 1.27 mm from each other, azimuthal planes every 20 of arc, and
axial planes every 4 mm of axial length. A 4 mm radius spherical volume
is digitally constructed around each reference position. The velocities
of all data points located within this volume are considered to be
correlated to the velocity of the fluid at the reference position, and
inversely, the velocity of the fluid at the reference position is cor-
related to the velocities of the data points in the volume. The amount
of influence each of these velocities has on the reference position
fluid velocity is based on the distance of the data point from the
center of the sphere which is the reference position. The reference
position velocities are determined by the equations,
N
E (U WTFi)
i=l
U =
X N
E WTF.
i=l
N
E (U Ri WTF )
S i=l
R N
E WTF.
i=l
N
E (U WTFi)
i=l
8 N
Z WTF.
i=l
where WTF is the weight factor and N is the number of data points within
4 mm of the reference position. The assumption has been made that the
velocity of a data point at the reference position would have the most
influence on the grid or reference position fluid velocity, while a
point 4 mm away would have no influence, and that there was a linear
decrease in influence. Thus, the weight factor was determined by
the expression,
WFT= (4-radius)
4
where "radius" is the distance between the data point and reference
point in millimeters. The reference positions closer to the pipe wall
than 4 mm have a skewed data problem that must be resolved. This is
done by adding a zero magnitude velocity vector at the same distance
from the reference point as the data point if the data point is further
from the reference point than the pipe wall is. This alleviates the
skewed data problem at the pipe wall by reducing the biasing influence
of the data point.
There should be a minimum of 2 or 3 data points in each volume of
influence or approximately 8 data points per cubic centimeter. This
means that in the 30 cm3 volume of this reference grid influenced area,
240 data points are necessary. The small prism viewing volume is ap-
proximately 50 cm3 so 400 data points or particle positions per frame of
cinematographic film are needed in the volume to use the reference grid.
This is 800 images on each cinematographic film frame which is four or
five times the data density of this work. This is the main reason that
this program was not used. It was mentioned earlier that reworking the
program for following particle images and matching image pairs should be
capable of such an increase in data density.
216-FOLLOW -FORT -SYSPRINT
LEVEL 2.1 ( JAN 75 )
OS/360 FORTRAN H EXTENDED
DATE 76*251/15.
REQUESTED OPTIONS: NODECK,NOLOADOPT=0*LC(48)
OPTIONS IN EFFECT: NAME(MAIN) NOOPTIMIZE LINECOUNT(48) SIZE(MAX) AUTODBL(NONE)
SOURCE EBCDIC NOLIST NODECK NOOBJECT MAP NOFORMAT GOSTMT NDXREF ALC NOANSF NOTER
ISN 0002
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0020
0021
0022
0023
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
C THIS PROGRAM WILL REFERENCE THE DATA TO A GRID SYSTEM IN THE PIPE
DIMENSION RATRO(23),THETAO(16)9XO( 9)*ADDU(23.15, 9),ADDV(23,.16, 9
1),AD.DW(23.16, 9),SUMWTF(23.16. 9) ,U23.16, 9).V(23,16. 9)*W(2316
1 9)
C SET INITIAL VALUES TO ZERO
DO 40 J=1,23
DO 41 K=1 16
DO 42 L=1-9
ADDU( JK,L)=0.0
ADDV(JKL)=0 0
ADDW(JK,L)=0.O
SUMWTF(JK ,L)=0.0
42 CONTINUE
41 CONTINUE
40 CONTINUE
C THE INITIAL VALUE OF IFRAM IS THE FIRST VALUE OF FRAME
IFRAM=2
C READ IN DATA CARDS
DO 50 I=1,.400
READ(5,105) IPARTIFRAMEX,R,THETA.UX*VR*WT
105 FORMAT(I 4,2XI2. 12X.6F6.2)
RATIOR=1.-R/63.5
C THIS PART ORDERS BY FRAME
86 IF(IFRAME.EQ.IFRAM) GO TO 88
C SOLVE FOR VELOCITIES AT THE GRID POINTS
DO 90 J=1,22
DO 9.1 K=1,15
DO 92 L=1,8
IF(SIUMWTF(J.K,L).EQ.3O.0) SUMWTF(J,KL)=1.
U( J.KL)=ADDU(JKL)/SUMWTF(JKL)
V( JKL)=ADDV(,JKt.L)/SUMWTF(JKL)
W(J,KoL)=ADDW(JSK.L)/SUMWTF(J.KL)
92 CONTINUE
91 CONTINUE
90 CONTINUE
C OUTPUT BY FRAME AND X WITHIN FRAMES
DO 94 L=1,8
WRITE(6.106) IFRAMXC(L)
106 FORMAT(1H1.,1OX,37HVELOCITY COMPONENTS =OR R=3500 FRAME=,12,2X,2HX=
1 F4.1)
WRITE(6,107)
FIGURE 4.7 Computer program to reference velocity field to known locations
MAIN: OS/360 FORTRAN H EXTENDED
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
ISN 0063
ISN 0064
ISN 0065
ISN 0066
ISN 0068
ISN 0069
ISN 0070
ISN 0072
ISN 0073
ISN 0074
ISN 0076
ISN 0077
ISN 0078
ISN 0079
ISN 0080
FIGURE 4.7 conti
107 FOR AT(1X, 3HY/R,42X,5HTHETA)
WRITE(6,108) (THETAO(K),K=1,5)
108 FORMAT(12X.F5.11,2XF5.1.12X,F5.1,12XF5.1 .12XF5.1)
WRITE(6,109) (RATRO(J) (U(J.KL,) V{J ,L) W( J K,L) K=1,5) J=1 22)
109 FORMAT(1X,F4.2,IX,IH( 3F5.1.1H),1H(e3F5.1,lIH) 1H(,3F5.1,1H).IH(,3F
15.,1 H),1H(, 3F5*,1 H))
WRITE(6,110)
110 FORMAT(/)
WRITE(6,108) (THETAO(K),K=6.10)
wRITE(6,109) (RATRO(J).(U(J.KL),V(J,KL),W(JKL),K=6,10),J=1,22)
WRITE(6,110)
WRITE(6,108) (THETAO(K).K=1 .15)
WRITE(6.109)(RATR)(J).(U(J.KL),V(J.KtL)W(JKL)K=11.15),J=1.22)
94 CONTINUE
IFRAM=IFRAM+1
IF(IPART*EQ.O) GO TO 93
87 00 43 J=1,23
DO 44 K=1 16
DO 45 L=i1 9
ADD(JK, L)=0.0
ADDV(JKtL)=0.0
ADDW(J,K,L)=0.0
SUMWTF(JKL)=0.0
45 CONTINUE
44 CONTINUE
43 CONTINUE
GO TO 86
88 CONTINUE
C THIS PART WILL FIND GRID THAT DATA POINT INFLUENCES AMD ADD THE
C WEIGHTED VELOCITIES
DO 60 J=1,22
RATRO(J)=( J-1 )*002+0.02
RO=( .-RATRO(J)) *63.5
IF(R LT.(R.O-5.0).OR.R.GT.(RO+5.0 ) GO TO 60
DO 70 K=1415
THETAO K) 16 .-2.*K
C THIS IS 5.5MM AT 1/2R AND 11MM AT R(63.5)
IF(THETA.LT.ITHETAO(K)-10.).OR.THETA*GT*(THETAO(K)+10*)) GO TO 70
0O 30 L=1,8
XO(L)=L*4.+5.
IF(X.LT.(XO(L)-5.).OP.X.GT.(XO(L)+5. ) GO TO 80
H=O
DX=XO(L)-X
DRH=RO*SIN(THETAO(K)*1.570795/90.)-R*SIN(THETA*1.570795/90.)
DRV=RO*COS(THETAO(K)*1.570795/90.)-R*CDS(THETA*1.570795/90.)
RADIUS=SOR.T(DX*DX+DRV*DRV+DRH*DRH)
nued
LEVEL 2.1 ( JAN 75 )
DATE 76.251/15
MAIN OS/360 FORTRAN H EXTENDED
ISN 0081 IF(R.ADIUSGT .4.0) GO TO 80
ISN 0083 IF((RO-R).GT.(63.5-RO)) H=(4.-RADIUS)/4.
ISN 0085 WTF=(4.-RADIUS)/4.
ISN 0086 ADDU(JKL)=ADDU(J.K,L)+UX*WTF
ISN 0087 ADDV(J,K.L)=ADDV(JKL)+VR*WTF
ISN 0088 ADDW(JK, )=ADDW(JK.KL)+WT*WTF
ISN 0089 SUMWTF(J.KL)=SUMWTF(JK,L)+H+WTF
ISN 0090 80 CONTINUE
ISN 0091 70 CONTINUE
ISN 0092 60 CONTINUE
ISN 0093 50 CONTINUE
ISN 0094 93 STOP
ISN 0095 END
FIGURE 4.7 continued
FIGURE 4.7 continued 0
LEVEL 2.1 ( JAN 75 )
DATE 76.251/15
4.6 Method for Displaying Reconstructed Flow Fields
At this stage of the data reduction process, a variety.of quanti-
tative information is recorded about the motion of the fluid through the
recorded flow field. This descriptive information is contained in the
instantaneous positions and velocities of the trace particles which
travel with the flow, and their time-space averaged statistical quantities.
The statistical quantities are in form for analysis but the instantaneous
quantities must be put into a more manageable form. The instantaneous
data points or reconstructed particle path sampling points are the posi-
tions of the particles for a series of arbitrary time intervals where the
cinematographic film rate determines the time intervals. Since the par-
ticle density and sampling rate are relatively high, there is a tremen-
dous amount of data to be analyzed even for short times. The easiest
method of displaying large quantities of instantaneous flow data is by
graphical representation. Two graphical methods of displaying the flow
field in time were chosen. The first method was plotting the velocity
vectors in the fixed volume at fixed times; the second was plotting the
reconstructed particle paths. The computer graphics unit was found to
be a handy way to pictorially display large quantities of data since it
does not require any manual operations. The graphics unit used was the
fast-draw general graphics package provided by the Gould Corporation.
The package consists of a number of library subroutines that are used
to construct and manipulate the graphics display, and the hardware to
generate the plot.
The instantaneous velocity fields are displayed at specific times
using a three dimensional perspective view, five R-theta field segment
views, and five X-R field segment views. This is a quick and easy method
for the illustration of flow events, since there are no manual operations.
The three dimensional perspective view or volumetric display shows the
entire experimental volume at any given instant in time while the R-theta
and X-R segment views clear up any confusion in interpretation of velocity
vector positions and orientations. The volumetric display makes two and
three dimensional motions easy to visualize since it is not necessary
to assimilate the instantaneous flow field from several two dimensional
plane views. Since only the volumetric view is necessary once any con-
fusing points are resolved, it is easy to follow events in time. An
example of the graphical displays drawn by this computer program are
shown in Figures 4.8, 4.9, and 4.10. Methods of display that require
the assimilation of several plane views are difficult to use to follow
three dimensional activity in time.
The velocity field computer program, which uses the fast-draw package
to construct its drawings, is constructed in three sections. The first
section draws the three dimensional velocity field (Figure 4.8). The
second section draws the R-theta segment velocity fields (Figure 4.9).
The third section draws the X-R segment velocity fields (Figure 4.10).
The second and third sections will be discussed together since their
formats are the same. This discussion will be a very general flowchart
of the program with emphasis on alterations for modifying the program.
The complete program is listed in Figure 4.12.
The three dimensional velocity field display consists of the velocity
vectors, a segment of pipe wall, and the viewed flow region boundaries.
All dimensions and .positions are a multiple of the actual size. This
proportionality is the scale factor which is arbitrarily set by the user.
The graphical representation of the flow region can be rotated to any
arbitrary viewing angle by altering the rotation angles in line 7 and 8
of the program in Figure 4.12. The angles (-80 and 600) were chosen to
give the clearest view of the velocity field. Drawing the pipe section
and flow region boundaries are the next operations performed by the
computer. The manipulation of the "pen" (movea, linea, etc.) are the
same ones that would be used when drawing the figure manually. The
last operation of this section is to draw the velocity vectors. The
position and velocity data are read by the computer at this point. The
positions are converted from millimeters to inches by multiplying by
"SCALE i". This conversion must be made because the fast-draw routines
only recognize dimensions in inches. Once the position is located, the
velocity vector is drawn. Its scale is arbitrary and set by multiplying
the velocity in millimeters per second by "SCALE 2". The fast-draw
package does not have an arrowhead subroutine for three dimensional
vectors so one was derived and programmed. The arrowheads were made
proportional to the length of the vector and located in the vector-viewing
vertical plane for clarity. The length is 30% of the vector with an
included angle of 30.
The R-theta and X-R plane segment views are used to clear up any
confusion in interpretation of the velocity vectors positions and orien-
tations in the three dimensional views of the velocity field as well as
to analyze flow events in two dimensions. Figures 4.9 and 4.10 are
divided into five segment views and a sixth view which superimposes all
the segments. The R-theta and X-R segments are 8 mm and 80 thick, res-
pectively. The superimposing of all segments in one sketch is used to
check for motions on the scale of the plane dimensions. Care must be
taken when viewing the segment displays so as not to be misled by large
differences in size and/or orientation of velocity vectors which are very
close together in a display. The vectors may be up to 8 mm away in the
R-theta segments or 80 away in the X-R segments from each other in the
direction of the normal to the plane. Segments this large were used
here because of the large scale of the turbulence and relatively small
quantity of data. The thickness of the segments can be varied by
changing the values of "LOLIMX" and "LOLIMT". When the lower limit
(lolim) equals the upper limit (uplim), the segments are reduced to
planes which is necessary when using data output of the reference grid
program. The lengths of the velocity vectors in the R-theta and X-R
plots are calculated by multiplying the plane velocities by "SCALE 3"
and "SCALE 4", respectively. The arrowheads in these sections are two
dimensional counterparts of the arrowheads in the first section of this
program. The two dimensional vector subroutine of the fast-draw package
draws arrowheads of constant size so it was not used in the present
displays.
The last computer program in the computerized data reduction method
draws the three dimensional reconstructed particle path display (Figure
4.13). The input data for this programare the reconstructed particle
paths from the three dimensional position and velocity program of section
4.4. This graphical display shows the motion of the fluid as it traverses
the experimental volume in the test pipe. A sample of the graphical dis-
play is shown in Figure 4.11 and the complete program is listed in Figure 4.13.
This display is the same section view with the same ability to be rotated
as the three dimensional velocity field display. The particle paths are
constructed by connecting the particle path sampling points with straight
lines. "SCALE 1" converts the dimensions of the positions from millimeters
to inches so the scale of the figure and the pathlines are the same factor
in the "INIT" statement's argument (line 2 in Figure 4.13).
- -I -
-9
-
/ -9
'-9
--9 r\
9-c-
r
#\J
K;l 4\y
FRRME 6
.FIGURE 4.8 Three dimensional velocity field
-9
(N-
/
/ /
/.I A -
.. 4-.
421
X=18-26
I ./
I I
THETR 20 -20
I #
0 THETA
X=34-U2
10 -
20 -20
X=RLL
THETA 20
-
Y/R
0.2
L
-2
S- I \1
FIGURE 4.9 R-eaplane
S.. I -I .
THETFI A -20 :
segment velocity fields
vI
_. 1,
THETR 20 -20
-.4
S TET
-20 0
I I
0
X=26-34
0
0
-I o
o0 0o
2
I
i .
-------
X=2-10
X=10-18
-20
-20
c
THETR=20-12
.-~- *- .
THETA=12-t4
THETR=U-(-4)
25 X
THETR= (-U) (-12)
JI
50 0
25 X
THETfA= -123 -20)
JI
50 0
I
25
THETAR=LL
S----__.
-b--- --~"
I I I
0 25 X 50 0
.FIGURE 4.10 X-R plane segment velocity fields
xI I50 .0 2
X 50 .0 25
'f/F
0.2
0.4
T/R
0.2
t"
-7
~ 'Y'
'4~a ~ag~t~~
----s
~c~-~L~sq~
~u iA~
5~-~p
C~b 7"
X 5
0 o
0
------35-~p
~%~-~
~--~-
~ccr
Sc3~n
--
-I
PARTICLE PATHS
FIGURE 4.11 Three dimensional particle path display
210-FOLLOW -FORT -SYSPRINT
LEVEL 2.1 ( JAN 75 ) DS/360 =ORTRAN H EXTENDED DATE 76.251/15
REQUESTED OPTIONS: NODECK9NOLOAD.OPT=0,LC(48)
OPTIONS IN EFFECT: NAME(MAIN) NOOPTIMIZE LINECOUNT(48) SIZE(MAX) AUTODBL(NONE)
SOURCE EBCDIC NOLIST NODECK NOOBJECT MAP NOFORMAT GOSTMT NOXREF ALC NOANSF NOTE
C FAST-DRAW PROGRAM TO PLCT VELOCITY FIELD IN PIPE (R=3500)
C IN 3-D VIEW AND PLANE CUTS (R-THETA) AND (X-R)
C THIS SECTION WILL PLOT 3-D VIEW OF FLOW FIELD
ISN 0002 DIMENSION X( 90),R(i90) ,THETA(90),UX(90),VR(90) WT(90),RATIO(90)
C TO CHANGE SIZE OF 3-D PLOT CHANGE *FACT' IN NEXT _INE
ISN 0003 CALL INIT (*SIZE'O28.*20.**FACT',2.3,*ORIG', 8B.10.)
C LABEL THE PLOT
ISN 0004 CALL TXSIZ(O100)
ISN 0005 CALL TEXT("FRAME 225')
ISN 0006 CALL TXPLT (-.75.-3.4,0*,0o)
C THE INITIAL ALIGNMENT IS A=X. B=Ys C=Z
C THE 2-D0 VIEW IS THE XY PLANE
C THE XY PLANE IS LOOKING DOWN AT THE PIPE WITH PRISM ON RIGHT
ISN 0007 CALL ROTAT (4,*AB'-80.)
ISN 0008 CALL ROTAT (5,*AC,.60.)
C USING THE DIMET AND ROTAT IN COMMENT CARDS WITHOUT THE FIRST TWO
C ROTAT CARDS KILL PLOT ISOMETRIC VIEW OF PIPE FROM DOWNSTREAM SIDE
C CALL DIMET (8,1.0)
ISN 0009 CALL PERSP (9, 0.* 0.,20.)
ISN 0010 CALL XZ
C CALL ROTAT (5*'AC*,90.)
C DRAW PIPE SECTION
ISN 0011 CALL RESERVE (0,0.4)
ISN 0012 CALL DEFINE (CIR (1).'CEN',0*,0*..SANG'. 30.,'FANG',330.,*CW*,
1'RADO,2.5)
ISN 0013 CALL DEFINE (CIR (2),'CEN*0*O,0.'SANG** 330..'FANG', 30..*CCW*,
1'RA)3,2.75)
ISN 0014 CALL DEFINE (CIR (3)9,CEN',0.,0.,*SANG' .30.,'FANG*,3309.'CW*,
I RAD'. 1.125)
ISN 0015 CALL DEFINE (CIR (4),*CEN'.0.*O. SANG'.,330.**FANG', 30.,'CCW,
I'RAD, 1.125)
ISN 0016 CALL MOVEA (2.165,0.*1.25)
ISN 0017 CALL- DRAW (CIR (1)')
ISN 0018 CALL MOVEA (2.381,90--1.375)
ISN 0019 CALL MOVEA (2.381,0.,1.375)
ISN 0020 CALL LINEA (2.1655.0*.125)
ISN 0021 CALL LINEA (2.165*-2.,1.25)
ISN 0022 CALL DRAW (CIR (1))
ISN 0023 CALL LINEA (2.381.-2*,-1.375)
ISN 0024 CALL DRAW (CIR (2))
ISN 0025 CALL LINEA (2.165,-2..1.25)
FIGURE 4.12 Computer program to construct velocity fields
MAIN OS/350 =ORTRAN H EXTENDED
ISN 0026 CALL MOVEA (2.381.-2.* 1375)
ISN 0027 CALL LINEA (2.381o.0.1.375)
ISN 0028 CALL MOVEA (2.16590..-1.25)
ISN 0029 CALL LINEA (2.165,-2.,--125)
ISN 0030 CALL MOVEA (2o165,.0.,-1.25)
ISN 0031 CALL DASH (.075*,075)
ISN 0032 CALL LINEA (.9743,0..-.5625)
ISN 0033 CALL DRAW (CIR(4))
ISN 0034 CALL LINEA (2.16b50. 1.25)
ISN 0035 CALL MOVEA (.9743,0...5625)
ISN 0036 CALL LINEA (.97431-2.,*5625)
ISN 0037 CALL LINEA (2.165,-2.*1.25)
ISN 0038 CALL MOVEA (.9743,-2.,*5625)
ISN 0039 CALL DRAW (CIR(3))
ISN 0040 CALL LINEA (2.165,-2.,-1.25)
ISN 0041 CALL MOVEA (.9743--2.,-e5525)
ISN 0042 CALL LINEA (.9743.0.,-.5625)
ISN 0043 CALL NDASH
ISN 0044 CALL MOVEA (0., 0.,0.)
ISN 0045 CALL LINEA (0.*-.10,C.)
C VELOCITY FIELD PLOTTING
C X1FAST-DRAW=YD(DATA).. DX FAST-DRAW=DYD
C Y FAST-DRAW=-XD(DATA), DY FAST-DRAW=-DXD
C Z FAST-DRAW=ED(DATA), DZ FAST-DRAW=DZD
C MUST HAVE *00000' CARD AT END OF DATA
C READ DATA CARDS
ISN 0046 ITOTAL=O
ISN 0047 DO 5 I=1,9.0
ISN 0048 READ(5100) Y,ZX(I),R(I),THETA(I),UX(1),VR(I),WT(I),DYDZ
ISN 0049 100 FORMAT(8X4.10F6.2)
ISN 0050 IF(X(I).EQ.O.) GO TO 6
ISN 0052 ITOTAL=ITOTAL+1
ISN 0053 RATIO(I)=1.-R(I)/63.5
ISN 0054 X1=X(I)
ISN 0055 DX=UX(I)
C SCALE CONVERTS MM TO INCH
ISN 0056 SCALE1=0.0394
ISN 0057 XI=KI*SCALE1
ISN 0058 Y=Y*SCALE1
ISN 0059 Z=Z*SCALE1
ISN 0060 CALL MOVEA (Y,-X1,Z)
C SCALE2 CONVERTS MM TO INCH AND SCALES VELOCITY
ISN 0061 SCALE2=0.005
ISN 0062 DX=3X*SCALE2
ISN 0063 DY=DY*SCALE2
ISN 0064 DZ=DZ*SCALE2
FIGURE 4.12 continued
LEVEL 2.1 ( JAN 75 )
* DATE 76.251/15
MAIN OS/360 =ORTRAN H EXTENDED
ISN 0065
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
0066
0067
0068
0069
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
ISN 0088
ISN 0089
ISN 0090
ISN 0091
ISN 0092
ISN 0093
ISN 0094
ISN 0095
ISN 0096
ISN 0097
ISN 0098
ISN 0099
ISN 0100
ISN 0101
ISN 0102
ISN 0103
ISN 0104
ISN 0105
FIGURE 4.12 -
CALL LINE (DY,-DXDZ)
C THIS PORTION WILL CONSTRUCT ARROWHEADS ON THE VECTORS WHICH ARE IN
C THE PLANE OF VECTOR-Z AXIS, 30% OF VECTOR LENGTH, AND 30 DEGREE
C INCLUDED ANGLE. TANGENTOF 1/2 ANGLE(15) IS4* 0.26795
A=DK*DX+DZ*DZ
B=SQRT(A*DX*DX+A*A)
C=0.
IF(B.EQ.O.0) C=1.
VU1=A/(B+C)
VU2= DX*DY/(B+CJ
VU3=-DY*DZ/( B+C)
R1=-.3*DY
R2=- 3*(-DX)
R3=-.3*DZ
D=0.3*(SQRT( A+DX*DX))*0.26795
I1X=R1+VU1 *D
P1Y=R2+VU2 *D
P1Z=R3+VU3*D
P2X=R1+VU1 *(-D)
P2Y=R2+VU2 *(-D)
P2Z=R3+VU3*(-D)
CALL LINE (P1X.P1Y,P1Z)
CALL MOVE (-P1X,-PIY.-P1Z)
CALL LINE (P2X.P2Y.P2Z)
5 CONTINUE
C THIS SECTION WILL PLOT 5 CUTS AND A VIEW OF ALL PLANE VELOCITY VECTORS
C FOR BOTH THE R.-THETA PLANE AND THE X-R PLANE
C DRAW PLOT OF R-THETA CUTS
6 CALL FACT (1.0)
CALL CLEAR. (4,5.9)
CALL XY
CALL ORIG (18.6915*2)
CALL TXSIZ (0.0875)
CALL TXPLT (-1.6*-1.844,0.,+1l,0.2a*)
CALL TXPLT (-1.6,-I.266,0. +1,Y/Ra )
CALL TXPLT (-1.6,-0.588,0.,+1,'0*.4*)
CALL TXPLT (-1.6,1.156,0.*+1,*C.2a')
CALL TXPLT (-1.6.1.734.0. +1,'Y/R@a)
CALL TXPLT (-1.6,2.312,0., +1,0.48*)
CALL TXPLT (0.012.70,0.,0*X=2-10@')
CALL TXPLT (3.0,2.70,0.,0.'X=10-18@*)
CALL TXPLT (6.0,2.70.0.,0**X=18-26a@ )
CALL TXPLT (0.0.-0.30,0.,0,'X=26-34@*)
CALL TXPLT (3.0,-0.30,0.,0,'X=34-423~).
CALL TXPLT (6.0,-0.30,0.*0s'X=ALL@')
DO 10 J=1,2
continued
LEVEL 2.1 ( JAN 75 )
DATE 75.251/15
MAIN OS/360 =3RTRAN H EXTENDED
.ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
ISN
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0115
0117
0118
0119
0120
0121
0122
0123
0125
0127
0129
0131
0133
0134
0135
0136
0137
0138
0139
0140
ISN 0141
ISN 0142
ISN 0143
ISN 0145
ISN 0146
ISN 0147
ISN 0148
ISN 0149
ISN 0150
ISN 0151
ISN 0152
FIGURE 4.12 -
DO 20 I=1,3
XORIG=18.6+3.0*(I-1)
YORIG=15.2-3.0*(J-1)
CALL ORIG (XCRIG,YORIG)
CALL MOVEA(-1.4.0.1)
CALL LINEA (-1.4,0.)
CALL LINE (1.4,0.)
CALL LINEA (1.40.1)
CALL MOVEA (0.,0.1)
CALL LINEA (0..0*)
CALL TXPLT (-1.5,-0.1.0.*+l1.,-2Oa*)
CALL TXPLT (0*.-0*.1,0*0,'0')
CALL TXPLT (0.7,-0.1,0.0,*THETA@')
CALL TXPLT' (1 .4-0.1,C.-1, *20a.)
UPLIMX=8.*I+2*+(J-1)*24.
LOLIMX=UPLIMX-8.
C DRAW VECTORS IN PLOT OF R-THETA CUTS
DO 30 K=1IITOTAL
IF(THETA(K).LT.(-20.).OR.THETA(K).GT.20.) GO TO 30
IF( ATIO(K) LT.O OR.RATIO(K) GT.*0.45) GO TO 30
IF(X(K).LT.2.OR.X(K).GT.42.) GO TO 30
IF(I*J.EQ.6) GO TO 31
IF(K(K).LT.LOLIMX.OR.X(K).GT.UPLIMX) GO TO 30
31 CONTINUE
C MULTIPLY RATIO AND THETA TIMES SCALE FACTORS TO =IT S ACE
VERTD=RATIO(KJ*5.78
HORIZD=THE.TA(K)*0*.07
C MULTIPLY VR AND WT TIMES A SCALE FACTOR TO ADJUST VECTORSTO PLOT SCALE
SCALE3=0.045
VERT'V=-VR(.K) *SCALE3
HORI ZV=WT(K)*SCALE3
CALL MOVE (HORIZDVERTD)
CALL LINE (HORIZVIVERTV)
C DRAW ARROWHEADS ON VECTORS OF 30% VECTOR LENGTH AND 30 DEGREE
C INCLUDED ANGLE. TANGENT OF 1/2 ANGLE(15) IS 0.26795
C=SQRT(VERTV*VERTV+HORIZV*HORIZVI
8=0.
IF(C.EQOO..) B=1.
VUI=VERTV/(C-e)
VU2=-HGRIZV/(C+B)
R1=-0.3*HORIZV
R2=-0. 3*VERTV
D=0.3*0.26795*C
PIANG=R1+VUI*D
P 1RAT=R2+VU2*D
P2ANG=RI+VU1 *(-D)
continued
LEVEL 2.1 ( JAN 75 )
DATE 75.251/15
OS/360 =ORTRAN H EXTENDED
ISN 0153 P2RAT=R2+VU2*(-D)
ISN 0154 CALL LINE (PIANG.P1RAT)
ISN 0155 CALL MOVEI (-PIANG,-PIRAT)
ISN 0156 CALL LINEI(P2ANG.P2RAT)
ISN 0157 30 CONTINUE
ISN 0158 20 CONTINUE
ISN 0159 10 CONTINUE
C DRAW PLOT OF X-R CUTS
ISN 0160 CALL ORIG (18.6,5.2)
ISN 0161 CALL TXPLT (-1.6,-1.844,0.,+1,0.2a')
ISN 0162 CALL TXPLT (-1.6,-1.266.0*,+1,*Y/Ra@)
ISN 0163 CALL TXPLT (-1.6,-0.688.0.*+1,90.4@*)
ISN 0164 CALL TXPLT (-1.6,1.156,0.,+1,*0.2@')
ISN 0165 CALL TXPLT (-1 6,1.734,0.*1, 'Y/RZ')
ISN 0166 CALL TXPLT (-1.6,2.312,0.,+1,*0.4a')
ISN 0167 CALL TXPLT (0.0.2.70,0.,0,THETA=20-12}*)
ISN 0168 CALL TXPLT (3.0,2.70,0.,0*THETA=12-4@')
ISN 0169 CALL TXPLT (6.0,2.70,0.,0.'THETA=4-(-4)a')
ISN 0170 CALL TXPLT (0.0,-0.30,0.0,*THETA=(-4)-(-12)@')
ISN 0171 CALL TXPLT (3*.0-0.30.0.,0,'THETA=(-12)-(-20)@')
ISN 0172 CALL TXPLT (6.0,-0.30,.0..0*THETA=ALL@')
ISN 0173 DO 40 J=1,2
ISN 0174 DO 50 I=1,3
ISN 0175 XORIG=18.6+3.0*(I-1)
ISN 0176 YORIG= 5.2-3.0*(J-1)
ISN 0177 CALL ORIG (XORIGYORIG)
ISN 0173 CALL MOVEA(-1.4,0.1)
ISN 0179 CALL LINEA (-1.4,0.)
ISN 0180 CAL. LINEA (1.4,0.)
ISN 0181 CALL LINEA (1.4,0.1)
ISN 0182 CALL MOVEA (0*,0.1)
ISN 0183 CALL LINEA (0.,0.)
ISN 0184 CALL TXPLT (-1.4,-0.1C..+1'O@*)
ISN 0185 CALL TXPLT (0.,-0.1.0.,0,'25@*)
ISN 0186 CALL TXPLT (0.7,-0.1,.0,0'X@')
ISN 0187 CALL TXPLT (1*.4-0.1.0.,-1**50@)
ISN 0188 UPLIMT=20.-(I-1)*8.-(J-1)*24.
ISN 0189 LOLIMT=UPLIMT-8.
ISN 0190 DXORIG=XOR;IG-1.4
ISN 0191 CALL ORIG (DXJRIGYORIG)
C DRAW VECTORS IN PLOT OF X-R CUTS
ISN 0192 DO 60 K=1,ITOTAL
ISN 0193 IF(THETA(KA).LT.(-2C.).OR.THETA(K).GT.20.) GO TO 60
ISN 0195 IF(RATIO(K).LT.,.OR.RATIO(K).GT.0.45) GO TO 60
ISN 0197 IF{X(K).LT,0.OR.X(K).GT.45.) GO TO 60
ISN 0199 IF(I*J.EQ.6) GO TO 61
FIGURE 4.12 continued
MAIN
DATE 76.251/15.
LEVEL 2.1 (
JAN 75 )
|
