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 )