ju of Standards

in. Bldg.

TECH NOTE

Computer Solutions for

Thermal-Acoustical Oscillations

In Gas-Filled Tubes

363

Q

"*.

O

^»tAU Of '

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

I

1

UNITED STATES DEPARTMENT OF COMMERCE

C. R. Smith, Secretary

NATIONAL BUREAU OF STANDARDS • A. V. Astin, Director

NBS TECHNICAL NOTE 363

ISSUED NOVEMBER 30, 1967

COMPUTER SOLUTIONS FOR THERMAL-ACOUSTICAL

OSCILLATIONS IN GAS-FILLED TUBES

M. T. NORTON AND R. C. MUHLENHAUPT

Cryogenics Division

Institute for Materials Research

National Bureau of Standards

Boulder, Colorado

NBS Technical Notes are designed to supplement the Bureau's regular publications program. They provide a means for making available scientific data that are of transient or limited interest. Technical Notes may be

listed or referred to in the open literature.

For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C., 20402

Price: 50 cents.

TABLE OF CONTENTS

Page

Abstract . . . „ „ „ 1

I. Introduction „ 1

2„ General Theory of Amplitude Calculations . 2

3. General Theory of Frequency Calculations ........... 7

4. Description of Amplitude Calculations. 11

5. Detail of Subroutine READ DATA in Amplitude Calculations . . 16

6. Detail of Subroutine ELEMENT in Amplitude Calculations ... 19

7. Detail of Subroutine HEAT TR ................ . 20

8. Detail of Subroutine AMPLITUD 23

9. Detail of Subroutine READ DATA in Frequency Calculations. . 35

10. Detail of Subroutine ELEMENT in Frequency Calculations . . 3 5

II. Detail of Program FREQENCY 37

12. Results and Conclusions. .................... 42

13. References • ........................... 52

14. Appendix A 53

1 5. Appendix B 66

16. Appendix C 68

17. Appendix D. Program THM LOS 70

18. Appendix E. Subroutine READ DATA for Amplitude ..... 71

19. Appendix F. Subroutine ELEMENT for Amplitude. ...... 73

in
TABLE OF CONTENTS (CONTINUED)

Page

20. Appendix G. Subroutine HEAT TR 74

21. Appendix H. Subroutine AMPLITUD ............. 76

22. Appendix I. Subroutine READ DATA for Frequency 80

23. Appendix J. Subroutine ELEMENT for Frequency ...... 81

24. Appendix K. Program FREQENCY .............. 82

List of Figures

Figure 1. Work Cycle of Simplified Thermal Oscillation Cycle . . 4

Figure 2. Typical Amplitude Determination ............ 6

Figure 3. Typical Calculated Oscillation Frequency 43

Figure 4. Test Data Courtesy of John Bannister, Air Reduction

Co. ............................ 44

Figure 5. Amplitude Versus Pipe Diameter Using Constant

Friction Factors and Text Book Heat Transfer

Coefficients. 46

Figure 6. Amplitude Versus Pipe Diameter Using Modified

Calculated Friction Factors and Text Book Heat

Transfer Coefficients. .................. 47

Figure 7. Amplitude Versus Pipe Diameter Using Modified

Calculated Friction Factors and 0.9 times Text Book

Heat Transfer Coefficients. ............... 48

Figure 8. Sonic Front Patterns at Odd and Even Calculation

Times. 54

Figure 9. Element Pair Surrounded by Sonic Fronts ....... 57

Figure 10. Sonic Front in Element N ........ 62

Figure 11. Sonic Front in Element 1 64

IV
COMPUTER SOLUTIONS FOR THERMAL-ACOUSTICAL

OSCILLATIONS IN GAS-FILLED TUBES

M. T. Norton and R. C. Muhlenhaupt

A digital computer program to determine solutions for thermalacoustical oscillations in gas -filled pipes is described. Using a typical

temperature gradient and tube length from test data, the program calculates the effect of changes in heat transfer coefficient, friction factor,

and pipe diameter. Details of the program are explained. A comparison

is made between computation and one particular test datum point. Although the calculations are based on helium as the media, the program

can accommodate any fluid treated as a perfect gas.

Key Words: Cryogenics, thermal oscillations, liquid helium

1. Introduction

When the open end of a pipe is at low temperature and the closed

end is much warmer, thermal-acoustical oscillations often occur. These

oscillations are common in cryogenic vent lines, pressure sensing lines,

and liquid transfer lines. This study is an investigation of this usually

undesirable phenomenon. The result is a digital computer program to

predict frequencies and overpressure amplitudes from pipe diameter,

temperature gradient, and gas properties. One test datum point is used

to establish the relation between text book friction factors and friction

level in an oscillating pipe.

Curves are presented showing the effect of pipe diameter, friction

factor, and heat transfer coefficient on the magnitude of the oscillations.

Although the program is specifically for helium gas, minor changes

can adapt it to other gases.

The calculations for frequency and amplitude are performed by

quite different methods. The frequency depends on acoustic velocity,
and acoustic frequencies in general are essentially independent of amplitude. So it is assumed that in thermal oscillations, the frequency is independent of amplitude and,for a given gas, depends only on the temperature

gradient. For convenience, when determining frequency, an initial pulse

of pressure is assumed.

The amplitude depends on the frequency, but the method of calculating is not the same as that for frequency. In amplitude calculations

various assumptions are made in order to approximate the shape of the

pressure profile.

2. General Theory of Amplitude Calculations

There are infinite varieties of duct configurations which can produce oscillations when inserted into a cryogenic dewar. Of these the

simplest and most common case is a constant diameter tube, closed at

the warm end, and open at the cold end. This is the same as a dewar

vent line with the shut off valve closed, or a pressure measurement line.

It is well known that these are subject to thermal oscillations. The

general theory applicable to this case is applicable to all cases.

The problem is not new nor is it confined to cryogenics. Glassblowers frequently experience it. A description of the phenomenon and

a qualitative explanation is found on page 230 of Rayleigh (1945).

A qualitative explanation of thermal oscillations, fundamentally in

agreement with Rayleigh, follows. During the pressure rise portion of

an oscillation cycle there is movement of gas from the open towards the

closed end. As a gas element travels, its pressure increases and, from

adiabatic compression, its temperature rises. The wall temperature

increases toward the closed end, and the gradient is steeper than the

temperature gradient in the gas. Therefore, there is a temperature

difference tending to transfer heat from the wall to the gas. Because heat
transfer requires time, first the pressure increases and then heat is

transferred. After the pressure in the element has peaked, the process

is reversed. Then there is a series of pressure drops and heat transfers

out of the gas. When plotted on temperature-entropy coordinates this

sequence of events has positive work area, a simplified version of which

is shown on figure 1. Thus a gas element in contact with a wall, whose

temperature is greater at the closed end, and whose gradient is steeper

than the adiabatic gas temperature rise, does work on its surroundings.

If each gas element is in contact with this type of temperature gradient

then, clearly, work is done by the gas elements in the pipe. Conversely,

any gas element in contact with a wall, whose temperature gradient is

reversed or is less steep than adiabatic compression, will have work done

on it. Furthermore, if all gas elements are in contact with a gradient

such that work is done by all elements, it is clear that the motive power

exists for continued oscillations. Likewise, if all gas elements are in

contact with a gradient such that work is absorbed by all elements, it is

clear that any accidental oscillations will quickly damp out. For the

intermediate case, where the wall temperature gradient sometimes causes

work to be done by, and sometimes on, the various gas elements, it is not

clear if chance oscillations damp out or continue. A digital computer

program can be used to sum up all the infinitesimals of work done in one

cycle in all the gas elements in the pipe.

Thus, in an oscillation cycle, positive work is a necessary criterion

for the existence ol thermal oscillations. If the work summation is positive, any accidental oscillations build up to an equilibrium steady state

condition or until heat transfer from the duct wall alters the temperature

gradient, in which case the oscillations may cease. If the work summation

is negative, any accidental oscillations quickly damp out.
kj

ENTROPY

Figure 1. Work Cycle Of Simplified Thermal Oscillation Cycle
In case the work summation is positive, at equilibrium it must exactly balance losses. There are two types of losses, friction and kinetic

energy, expelled into the dewar. Part of the friction loss is simply wall

friction opposing the gas motion in either direction. Also, when the gas

enters the pipe from the dewar there is turbulence created which is equivalent to additional wall friction. In this study the entrance loss will be

handled as an additive to the wall friction.

In one pressure cycle each element will undergo heat transfer and

will do some positive or negative work. As each element slides along the

wall or into the duct, there is frictional return of work to heat in the gas

element. After summing up for all elements in one cycle, the excess of

positive work done over frictional work must equal the kinetic energy

ejected into the dewar during one cycle.

This criterion, that the net work done after subtracting friction

must equal the kinetic energy ejected, is used in an iterative calculation

to determine the pressure amplitude of helium gas in a given pipe with a given temperature gradient. The method used is to estimate the pressure

amplitude and calculate for all gas elements in the duct, the work done

and the friction and kinetic energy ejected. If the net work exceeds the

kinetic energy, the estimate is low. If the net work is less than the kinetic

energy, the estimate is high. With a few iterations it is possible to determine the pressure amplitude for which the net work equals the kinetic

energy ejected. This is the calculated pressure.

Figure 2 illustrates the method used. The particular pipe diameter

under investigation is .3175 cm (1/8 inch). In this illustration, pipe friction is calculated using friction factors that remain constant during each

run, and heat transfer is computed using unmodified heat transfer coefficients from McAdams (1954). When the overpressure estimate is low,

the ratio of ejected kinetic energy to work done is less than unity. When
»-

o

ii a t

»~ o E

.2 io m

^

h

o

CD

E

13 X

Heat

Transfe

Coefficient

Multipl

Frequency

=

31.

2

Diameter

=

3.1* %

<^\

°oIT)

I

CM

ro

E

o

(/> c

LUor

id

c/) LU

or

o_

orLU

>

o

00

o

c

•r-l eu

<D Q

0)3

OhS

<

t—

i

o

•t-< Oh H

tMu

o

cvj 00 
 CvJ COd

>W0M 01 3>1 11X3 -OllVd
the estimate is high, the ratio exceeds unity. Numerical interpolation

helps determine the final calculated overpressure.

3. General Theory of Frequency Calculations

Except for temperature gradient and heat transfer effects, the

pressure pulses of thermal oscillations are the same as those of water

hammer. As described by Rich (1951), one common case of water hammer involves pressure pulses traveling in a pipe with its open end in a large reservoir and its closed end at a valve. In that case the initial

condition is liquid flowing in a long pipe from a reservoir through a valve

to equipment of some sort or to discharge. If the valve is suddenly

closed, there is a sequence of four steps which cyclically repeat themselves until damped out by friction. These steps are as follows.

Step one starts just after the valve closure. Liquid velocity in the

immediate vicinity of the closure is transformed into pressure. This

pressure pulse travels to the open reservoir end at sonic speed. When it arrives at the reservoir, the pipe is momentarily filled with stationary

liquid whose pressure is higher than that in the reservoir. The amount

of pressure rise in step one is given by

AP = Cp |U|

where

2

AP = pressure rise ------- dynes/cm

C = sonic velocity ------- cm/sec

3

p = fluid density -------- g/cm

|U| = magnitude of fluid velocity- cm/sec.

Step two starts just when the overpressure wave of step one

reaches the reservoir. The overpressure at the open end drops to

reservoir pressure. In the process it acquires a velocity equal to its
original velocity but in the opposite direction. This reverse velocity and

reservoir pressure wave travels from the reservoir to the closure at

sonic speed. When it arrives at the closed end the pipe is momentarily

filled with liquid at reservoir pressure, whose velocity is everywhere

the reverse of the original.

Step three starts as soon as the reverse velocity wave arrives at

the closed end. Liquid velocity in the immediate vicinity of the closure

is transformed into pressure drop. This reduced pressure pulse travels

to the open reservoir end at sonic speed. When it arrives, the pipe is

momentarily filled with stationary liquid whose pressure is lower than that

in the reservoir. The amount of pressure drop is the same as the pressure rise of step one.

Step four starts just when the underpressure wave of step three

reaches the reservoir. The underpressure at the open end rises to reservoir pressure and acquires the original velocity. This pressure and

velocity travels from the reservoir to the closure at sonic speed. Momentarily at the end of step four the liquid pressure and velocity are the same

as they were at the beginning of step one.

After step four the process repeats itself until friction damps out

the oscillations..

Except for the effects of heat transfer and gas temperature gradient,

the process is similar for thermal oscillations. Because of liquid evaporation or other causes, a small pressure rise occurs in the dewar. Momentarily the tube is filled with stationary gas whose pressure is lower than

that in the dewar. This is the same situation that exists in the pipe in the

water hammer case at the end oi step three. Again ignoring the effects of

heat transfer and gas temperature gradient, the thermal oscillations

follow steps four, one, two, three of the water hammer case.
If thermal oscillations followed the foregoing sequence of events

exactly, a comparatively simple theory ot heat transfer from the warm

regions of the tube to the dewar could be derived. Even though the actual

case is more complicated, an examination of the simple case is instructive,

During step four every particle of gas in the tube is pressurized a small amount and travels a short distance into a warmer region of the

tube. Some heat is transferred after pressurization. During step one

the gas is pressurized once more and again moves into a still warmer

region. More heat is transferred. During step two the gas is depressurized and moves into a cooler region. Some heat is transferred. Once

equilibrium has been obtained the final condition at the end of step three

and the original condition at the start of step four are identical.

This chain of events is shown on temperature, entropy coordinates

on figure 1

.

A fluid which follows the path shown must do work on the surroundings. Part of the work done is that required to overcome wall and internal

gas friction to maintain the process. The balance oi the work done is expended in expelling gas from the tube to the dewar. This gas expelled into

the dewar contains the energy that increases the liquid evaporation rate.

If the behavior oi a thermal oscillation were exactly the same as

that of water hammer, the frequency could be determined by dividing the

tube into regions of known temperature, determining the average sonic

velocity in each region, and calculating the time required for a sonic

pulse to travel from one end of the tube to the other. The period of the

oscillation by this method would be the time required for the pulse to

traverse the tube four times, and the frequency would be the reciprocal

of the period. This method was applied to the data of Bannister (1965).

As is shown in table I, the method is not sufficiently accurate.
Position Length Temp. Sonic Vel

Inches cm °K cm/sec

0-8 20.3 4. 3 12. 2 x 10^

14.4 x 10

20. 38 x 10

55.8 x 10

91. 5 x 10

101. x 10

8-16 2 0. 3 6

16-24 20. 3 12

24 - 32 20. 3 90

32 - 40 2 0. 3 242

40 - 48 20. 3 295

TABLE I Using Temperatures of Figure 2

Time

sec

1.664 x 10";?

1.411 x 10"

.999 x

10"

. 364 x 10"

. 222 x 10"

.201 x

10"

-3

Time for one traverse 4. 861 x 10 sec

Period 0. 01944 sec

Frequency calculated 51.4cps

Frequency measured 34 cps

Since the calculated frequency of table I does not agree with measurement, it is evident that the method is inadequate. The reason is

simply that a temperature gradient exists in the pipe and causes pressure

wave attenuations and reflections. Therefore, the peak pressure at any

point is a summation of multitudinous reflections and does not necessarily

advance at sonic velocity, even though the first warning of a pressure

pulse does travel at sonic speed. To calculate frequency it is found convenient to divide the pipe into a large number of constant temperature

regions, assume the initiation of a small pulse in the dewar,and calculate

the resulting pressure versus time throughout the pipe from equations

derived in Appendix A.

In the derivations, all pressures are considered as absolute. However, since the derived expressions show pressure differences, it is often

possible in the digital computer program to ignore the absolute pressure

level and to consider only over or underpressure. Of course, in determining gas densities, the absolute pressure must be used.

10
These derived relationships have been programmed for calculation

on a high speed digital computer. As a first operation the program divides

the pipe into gas elements such that each has an equal sonic traverse time.

The number of these elements may be varied. The length of time required

for sound to traverse one tube length is one calculation time. Temperature

gradients are obtained from Bannister (1965). The original pressure level

is an input variable. The amplitude of the assumed pressure pulse is

2

arbitrarily assumed to be 1000 dynes/cm , although other values may be

used. For as many calculation times as desired, the program calculates

and prints gas velocities and pressures for each of the N elements. By

examination of the printout it is possible to determine the pulse frequency.

The program presented herein determines the period between each of the

first 10 pressure minimums thereby obtaining 9 periods, and 9 corresponding frequencies. The average of these 9 frequencies is then considered

to be the true frequency. Minor modifications can increase the number

of individual periods and frequencies averaged.

4„ Description of Amplitude Calculations

For amplitude calculations some assumptions were made as to the

nature of the pressure oscillations. It was noticed during frequency calculations that even a step pressure input resulted in almost sinusoidal

pressure variation at the closed end. Furthermore, all peak pressures

tended to occur at the same time, although the amplitude varied from end

to end. From oscilloscope observations it has been noticed that the oscillations tended towards sinusoidal. Therefore, it is assumed that at any

position in the pipe the pressure oscillation is sinusoidal. The end to end

variation of the maximum overpressure is obtained from a simple cosine

curve and the fact that at the open end there is a pressure node at which

the maximum overpressure is always zero.

11
For amplitude determinations the duct is divided into a large

number of elements, all of which have the same length. To account for

gas entering and leaving the pipe an additional quantity of elements is

assumed to exist in the dewar. Thus, there are N elements in the pipe

and M elements altogether. The period of one cycle is divided into 360

degrees. One time interval is taken to be a small fraction of the total.

Usually, a 5 degree time interval is selected as being small enough for

accuracy without requiring excessive computer cost.

Because of pressure change in the element nearest the closed end,

element one, its midpoint will oscillate back and forth in the duct. Its

temperature will change due to compression and heat transfer. The point

on the wall in contact with the midpoint will change with the movement,

resulting in temperature difference and heat transfer. There also will be

density and position changes that are somewhat different from an adiabatic

case.

Initial approximate temperatures and densities are assumed.

After heat transfer has been calculated repeatedly around one complete

cycle, conditions are different than initially assumed. However, after

a sufficient number of calculation cycles, equilibrium is obtained. Then

the position of the cold end of element one is a series of points varying

in a cycle, but invariant from one cycle to the next.

This final series of positions of the cold side of element one is the

permanent series of positions of the warm side of element two. Initial

approximate temperatures and densities are assumed in element two followed by the same type of heat transfer calculations that were performed

for element one. After numerous cycles of calculations, a series of positions is obtained for the cold side of element two. These are the permanent

series of positions of the warm side of element three.

12
In general this same series of calculations is performed for each

of the M elements. In each time interval the heat transfer, DQ, and the

infinitesimal work done, DWK, are computed. The summation of these

quantities during the final cycle gives the net heat transfer into,and

the net work done by, the I'th element, Q(I) and WORK(I). Because of

the oscillation, gas is continually entering and leaving the duct. As a result there are three types of elements. One type, like element number one,

is entirely inside the duct both at the beginning and at the end of a time

interval. A second type is partly inside the duct and partly in the dewar,

either at the beginning or at the end of the time interval, or both. The

third type is out of the duct and inside the dewar at all times.

For the first type, the calculations, though involved, are relatively

straightforward. For the third type, practically no calculations are

needed. For the second type it is necessary to divide the gas element into

a section inside the pipe and a section outside the pipe, at the beginning

and at the end of the time interval. As a result there are three routes

through the amplitude program, one for each type of calculation.

Using Fortran notation the cosine variation of overpressure with

distance along the pipe is expressed as,

where

PX(I) = PM*COSF(PI/2. *(XAVG(2)-XX(1)/XL)

2 PM = Maximum overpressure at closed end - dynes/cm

PI = 3.1415926536

XAVG(Z) = Position of element center - cm

XL = Length of duct - cm

PX(I) = Maximum overpressure in element (I) at

2

position XAVG(2) - dynes/cm

13
XX(1) = Position of closed end of duct - cm.

The above relationship holds if (XAVG(2)-XX( 1)) is less than XL. Otherwise PX(I) becomes zero.

The momentary overpressure at the end of a given time interval

in the cycle is

P(2) = PX(I)*SIN2

where

2

P(2 ) = overpressure at end of interval - dynes/cm

SIN2 = the sine function of the fraction of the cycle completed

at the end of the interval.

One of the primary variables calculated during each time interval

is the infinitesimal work done. The obvious way to make this calculation

is to multiply the average pressure during the interval by the volume increase during the interval. In Fortran notation this is

DWK = PRESS*(VOL(2)-VOL(l))

where

DWK = infinitesimal work done in the interval - ergs

PRESS = average absolute pressure in the interval - 2 dynes/cm

3

VOL(2) = element volume at end of the interval - cm

3

VOL(l) = element volume at beginning of the interval - cm . An instability in the program was finally traced to this method of computing work. An equivalent expression, derived in Appendix B, was

finally used. The equivalent method produced infinitesimal work terms

that at equilibrium were the same as those of the original method with an

error less than one part in 5000, and the instability was largely eliminated.

This method of calculation, in Fortran notation, is

14
DWK = (XK-1. )/XK*DQ-VOLM*(P(2)-P(l))/XK-DFR/XK

where

XK = specific heat ratio

DQ = infinitesimal quantity of heat transferred in a time

interval - ergs

3 VOLM = average element volume during interval - cm

P(2) = element overpressure at end of time interval - dynes/cm

P(l) = element overpressure at beginning of interval - dynes/cm

DFR = infinitesimal frictional energy loss during interval -

ergs.

Kinetic energy is calculated in two ways. One method is to deter

mine the infinitesimal kinetic energy in a given element every time the

velocity is towards the dewar, that is to say when it is positive. This is

done with the equation

DKE = 0. 5*(VEL*AF*ROE*DT)*VEL*VEL

"where

DKE = infinitesimal kinetic energy - ergs

VEL = gas velocity - cm/sec

2 AF = pipe cross section area - cm

3 ROE = average density - g/cm

DT = time interval - sec.

The terms in the parentheses are the mass of gas moving past a point

during the interval, so the equation is the usual one in which kinetic

energy is one half the product of the mass and the velocity squared.

By summing up the DKE for each time interval, the dewar bound

kinetic energy in one cycle for each element, ZKE(I), is found. The

15
kinetic energy of the N'th or M'th elements, ZKE(N) or ZKE(M), should

be representative of the kinetic energy expelled into the dewar.

A more correct value of the kinetic energy expelled into the dewar

is that at the exit from the duct. During one cycle, during each time interval, some one element is partly in and partly out of the duct. By summing the infinitesimal amount of kinetic energy of these elements, a

close value of the kinetic energy expelled into the tank is obtained. The

equation used is

EKE(J) = 0. 5*(VEL*AF*ROE*DT)*VEL*VEL

where

EKE(J) = infinitesimal kinetic energy at the duct exit during

the J'th time interval - ergs.

Kinetic energy is obtained in the program by both methods.

Approximately the same summations are obtained from each method.

The values obtained from summation of EKE(J) are used because they

more closely reflect actual operation,,

The program for calculating amplitude is called THM LOS. Since

this program does nothing except call the subroutines, READ DATA,

ELEMENT, and AMPLITUD, there is no need for discussion oi it, but

there follows a detail description of the subroutines.

5. Detail of Subroutine READ DATA in Amplitude Calculations

This subroutine is used to read information into the program and

to calculate certain derived quantities. Most of the input information is

printed out to show the calculation conditions. In the following description no explanation is made of Fortran statements which merely have

to do with the mechanics of the digital computer program such as COMMON, TYPE REAL, or FORMAT statements. Similarly no explanation

16
is made of statements which themselves are merely explanations. Subroutine READ DATA is given in Appendix E.

In line 5 a number of quantities are read which have the following

meanings:

PSTART = dewar pressure - dynes/ cm

DIAM = the duct diameter. (This is read in inches and

converted to centimeters)

FRF = Fanning friction factor

FACTOR = a multiplier used if desired to alter text book

heat transfer coefficients

AVGFREQ = the frequency of the oscillations - cycles/sec

PM = the assumed maximum overpressure - dynes/cm

N = the number of elements in the duct

NJJ = the maximum number of iterations permissible

NFR = a multiplier used if desired to alter text book

friction factors

JD = the number of degrees in one time interval. ( JD

can have the values 1, 2, 3,4, 5, 6, 8, 9, 10. The

value 5 has usually been used. One cycle is 360

degrees, so JD is a fraction of a cycle).

Since the program determines the oscillations that result from a given temperature gradient, temperature versus duct positions are needed.

This can be from test data, or from an arbitrary gradient to determine its

oscillation characteristics. In line 10 the temperature, TMP, in degrees

Kelvin is read for each of 32 duct positions, S, which are read in inches.

Viscosity of the gas, VIS, in micropoise is read at 16 temperatures,

TMRVIS, in degrees Kelvin at line 15.

17
Thermal conductivity, CONGAS, in milliwatts/cm- °K is read

at 12 temper atures, TCONGAS, in degrees Kelvin at line 20.

Viscosity and conductivity are obtained from Johnson (I960) either

directly or by extrapolation.

To determine natural convective heat transfer coefficients, a

table of Nusselt number versus the Grashof-Prandtl product is needed.

The log of the Nusselt number, ZNU, is read at 16 values of the log

of the Grashof-Prandtl product, GRPR, in line 2 5. Figure 7-10 from

McAdams (1954) supplied 13 of these values. The two largest values

used, not obtainable from figure 7-10, were obtained from figure 7-7 of

the same reference. Since the Nusselt number goes to zero if the

Grashof-Prandtl product goes to zero, it was desired to approximate this

condition. To do so the log of the Nusselt number is set to -37 when

the log of the Grashof-Prandtl product is -37.

In lines 30 and 40012 the duct positions, S, and the duct diameter,

DIAM, are converted from inches to centimeters.

In lines 35 and 40011 the viscosity, VIS, is converted from micropoise to poise.

In lines 40 and 40010 the conductivity, CONGAS, is converted from

mw/cm-°Kto erg/sec -cm- °K.

In line 45 the molecular weight of helium, WT, the universal gas

constant, RU, and the value for PI are entered. From these the helium

gas constant, R, and the cross section area of the duct, AF, are calculated. In line 50 the specific heat ratio, XK, is entered, and the specific

heats, CPG and CVG, are calculated. In line 55 the temperature oi the

gas in the dewar just outside the duct, TDEW, is entered and the density,

RHODEW, is calculated. The quantity, M, is also calculated. Without

oscillations there are N gas elements in the duct. With oscillations,
elements or portions of elements enter the duct from the dewar and

enter the dewar from the duct. To keep track of all these elements an

additional quantity of elements is considered. Thus, M represents a number of elements some of which are always in the duct, some are

partly in and partly out of the duct, and some are always in the dewar.

The remainder of the subroutine READ DATA is used for printing

out the input information. Computational units are printed except in the

case of viscosity, VIS, where convenience and space cause it to be printed

in micropoise rather than poise.

The input value of the Fanning friction factor, FRF, is not printed

in READ DATA. As explained later, the input value is subject to change.

The range of values actually used is printed from subroutine AMPLITUD.

6. Detail of Subroutine ELEMENT in Amplitude Calculations

The subroutine ELEMENT, when used with amplitude calculations,

divides the duct into N elements and M elements total in the duct and

dewar. All elements are of the same length, DXX. The position of each

boundary, XX, the wall temperature at each boundary, TEMP, the initial

average temperature of each element, TMG, the initial average density oi

each element, RHO, the mass of gas in each element, GM, the initial midpoint of each element, XMID, are also calculated in ELEMENT and

printed. These are all initial or non-oscillating conditions. This subroutine, ELEMENT, is given in Appendix F.

In line 5 the warm boundary temperature of the first element,

TEMP(l), is set equal to the temperature read in at that point, TMP(l).

Similarly the warm boundary, XX(1), is set equal to S(l). The quantities

NP1 and MP1, needed later, are set equal to N+l and M+l. Line 15 sets

the cold end boundary, XX(I+1), of elements 1 through N. Lines 20 through

33 compute the temperature, TEMP(I+1), corresponding to each position,

19
XX(I+1). Lines 35 and 40 calculate the average initial element temperature, TMG(I), the average initial density, RHO(I), the initial average

position, XMID(I), and the mass of gas in each element, GM(I)„ Lines

45 through 60 compute the same quantities for elements N+l through M.

The remainder of the program prints out the values of these quantities.

7. Detail of Subroutine HEAT TR

When called by subroutine AMPLITUD, subroutine HEAT TR calculates the friction factor, FRF, the Reynolds number, REYNO, the heat

transfer coefficient by forced convection, HHT1, the Grashof-Prandtl

product, GP, the heat transfer coefficient by natural convection, HHT2,

the total heat transfer coefficient, HHT, and the heat transferred, dQ.

When it is desired not to calculate the friction factor but to use a constant

input friction factor, the subroutine omits calculation of friction factors.

Subroutine HEAT TR is given in Appendix G.

In line 1, if the variable, ITESTFRF, has been set equal to unity

in subroutine AMPLITUD, only friction factors are desired and the program jumps to line 386.

Lines 2 through 3 calculate the momentary conductivity, CONG, of a gas element from the momentary element temperature, TAVG, by

linear interpolation in tables which have been read in the program.

Lines 4 through 5 calculate the momentary viscosity, VISC, of the

gas element from the momentary element temperature, TAVG, also by

linear interpolation.

At line 6, if the viscosity is zero, an error exists and the program

is terminated.

In line 8 the Reynolds number, REYNO, is calculated using variables determined in AMPLITUD. Since a Reynolds number is always

positive, its absolute value is used.

20
In line 10 the lowest possible value of the Reynolds number is set

to 0. 001 to avoid the possibility of division by zero. Furthermore, when

the Reynolds number is less than 0. 001 forced convection is usually

negligible- Friction may be appreciable, but this assumption merely

forces the friction factor to a single high value when the Reynolds number

would be below 0. 001.

In line 22, if the conductivity of the gas is negative or zero, an

error exists and the program terminates.

In line 25 the Prandtl number is calculated. If it is zero or negative,

an error exists and line 26 initiates program termination.

In line 30, if the momentary gas temperature, TAVG, is negative

or zero, an error exists and the program terminates.

Line 31 calculates the Nusselt number, ZNUSSELT, by equation

(9-10a) of McAdams(1954). From this, line 32 determines the forced

convective heat transfer factor, HHT1. The multiplier, FACTOR, is an

input quantity used to determine the effect on oscillations of reduced or

increased heat transfer factors.

Line 33 determines the heat driving temperature difference,

THETA. If it is zero, the program skips unnecessary calculations and

goes to line 385.

Lines 35 through 42 calculate the log of the momentary GrashofPrandtl product, GPL. If the Grashof-Prandtl product, GP, is negative,

an error has occurred and the program is sent to line 410 and terminated.

If it is zero the program goes to line 380 and skips unnecessary calculations.

Lines 43 through 370 determine the login °"^ the free convective

Nusselt number, ZNUL, and line 375 determines the corresponding

Nusselt number, ZNUS.

21
In line 380 the free convective heat transfer coefficient, HHT2, is

determined. The multiplier, FACTOR, is an input quantity used to determine the effect of increased or decreased heat transfer factors on oscillations. The overall heat transfer coefficient, HHT, is obtained bysimple addition of forced and natural heat transfer coefficients.

In line 385 the heat transfer coefficient, HHT, is multiplied

by the heat transfer area, AREA, the temperature difference, THETA,

and the time increment, DT, to obtain the infinitesimal heat transferred,

DQ, for the element in the time increment.

Lines 386 through 39 5 have to do with friction factor. In line 386

the minimum value of the Reynolds number is limited to 0. 001.

In line 2386 the decision is made whether or not to calculate the

friction factor. The program has the option of using a constant input

friction factor or oi using a modified text book friction factor. If the

value of the friction factor during input is zero, then in subroutine

AMPLITUD the variable, IFRF, is set to zero. Alternately, if the input

value of the friction factor is not zero, the variable, IFRF, is set to unity.

Thus, line 2386 causes calculation of a friction factor if the original input

friction factor is zero. Otherwise, it sends the program to line 99999

and makes no friction factor calculations.

If the Reynolds number is less than 1000 in line 387, the friction

factor is obtained from the usual laminar equation, but multiplied by NFR.

The quantity, NFR, is an input multiplier used to explore the effect of

increased or reduced friction. If the Reynolds number is less than 1000

in line 388, the friction factor has already been obtained, and control

returns to AMPLITUD.

If the Reynolds number is greater than or equal to 1000, lines

389 through 39 5 calculate the friction factor by an equation (6-8e) of

22
McAdams (1954). The multiplier NFR is used as in laminar flow. The

method is iterative. A value for friction factor is assumed in lines 389

and 390. A better estimate is found in line 392. If this estimate disagrees with the previous estimate by an appreciable amount the program

returns to line 390 and iterates again. Otherwise, the calculation of

friction factor is complete.

Lines 400 through 420 are used in case of error. In the event of

an unstable calculation, TAVG, CONG, or GP can become negative. If

this occurs, there is a print-out to aid in diagnosis and the quantity, DT,

is set equal to zero, which causes termination of the program.

8. Detail of Subroutine AMPLITUD

This subroutine is the fundamental part of the program. All the

other subroutines are for the purpose of helping subroutine AMPLITUD.

As an element changes pressure, its density changes, heat is transferred

from or to the wall, work is done by the element, friction is experienced,

and at the open end, kinetic energy is lost from the pipe and ejected into

the dewar. From an assumed value of the maximum overpressure at the

closed end, AMPLITUD computes the behavior of each element in turn

from 1 to M. Subroutine AMPLITUD is given in Appendix H.

The line before line 1 determines whether or not the program will

use input or calculated friction factors. If the input value of friction

factor, FRF, is zero,the parameter, IFRF, is set equal to zero. If the

input value is different from zero, the parameter, IFRF, is set equal to

unity. Then, if the input friction factor is zero, IFRF is zero, and subroutine HEAT TR computes friction factors for each element in each time

interval. Alternately, if the input friction factor is any desired non-zero

quantity, IFRF is unity, and subroutine HEAT TR leaves the friction factor

unchanged at its input value.

23
Lines 1 through 40009 cause the printing of a representative table

of friction factors used in the program versus Reynolds number. Where

a constant friction factor is used, this same friction factor shows for each

Reynolds number. Where a calculated friction factor is used, the friction

factor calculated for the corresponding Reynolds number is printed. The

friction factor printed includes the multiplier NFR for the calculated case. When the input friction factor is used, NFR is inoperative.

Lines 5 through 25 set the values of certain constants needed during

calculations or set initial values for variables to insure proper operation

of the program. The most important of these is the time interval, DT,

calculated in line 10. The time of one interval is the time for one cycle

times the fraction of a cycle used. Since the time for one interval is the

reciprocal of the frequency, AVG FREQ, and since the cycle fraction is

JD/360, in Fortran notation, the time interval is given by

DT = JD/360. /AVG FREQ.

Line 3 is a DO LOOP ending in line 40000 which causes the calculation of all quantities needed for each of the elements 1 through M.

Lines 31 through 55, for each gas element in turn, set initial values

for variables to insure proper operation.

In the DO LOOP from line 60 to line 40002, the calculation for each

element is caused to repeat as many times as needed. The maximum number of repetitions is the quantity NJJ, an input quantity. The value of NJJ,

in general, is kept high enough to insure that conditions calculated for the

final cycle are substantially the same as for the previous one. The quantity,

IQ, is used to indicate the adequacy of this convergence. Later in the program, if cycle convergence is adequate, IQ is increased by 1. If this

happens twice, the accuracy is sufficient and the program proceeds to the

next element.

24
In line 6l, the maximum number of cycles, JJ MAX(I), used for

the gas element I is determined. Since this quantity is printed, it is

possible to see if the input quantity NJJ is adequate.

In line 65 the quantity QI is set equal to Q(I) and WORK I to WORK(I)

for later use in determining adequacy of convergence. Later the sum of

all infinitesimal heat transfer to one gas element, Q(I), and the sum of

the work terms in one gas element, WORK(I),is obtained. At that time

QI and WORK I become the heat and work summations for the previous

cycle and a comparison is then possible.

Other calculations in lines 65 and 70 set initial values to insure

proper program operation.

Line 80, a DO LOOP ending in line 40003, causes calculations for

one cycle to occur for the element already selected in line 30. It sets the

variable J equal to the portion of a cycle, expressed in degrees, completed

at the end of this time interval. The first time interval is zero to JD degrees (If JD is 5 then the interval is 0-5 degrees. ) and so on until the

cycle is completed at 3 60 degrees. Because of line 80, the input variable

JD must be divisible into 360. /lso computed in line 80 is the sine function

of the angular fraction of the cycle completed at the end of the time interval.

The quantity, JLJD, refers to the time (or angular fraction of the

cycle) at the beginning of the interval just as the quantity, J, refers to it

at the end of the interval. In line 85 it is calculated. At the beginning of

the cycle J = JD, JLJD = 360, but at any later interval when J > JD,

JLJD = J-JD.

In lines 90 through 105 numerous variables at the beginning of the

current interval are set equal to the same variable at the end of the previous interval. Thus, RO(l), the gas density at the beginning of the interval,

is set equal to RO(2), the density at the end of the previous interval. Later,

25
a value of the density at the end of the current interval is determined,

which then becomes RO(2). The same is true for the volume of the element, VOL, the average position of the element, XAVG, the overpressure

of the element, P, the absolute pressure of the element, PRES, the mass

of the gas of an element that is in the pipe, G, and the mass of an element

that is in the dewar, GID2 and GID1.

Line 110 determines if a gas element is entirely outside the pipe

and in the dewar by means of the variable, X. Variable, X, is similar

to XX. XX is computed in ELEMENT and is the non-oscillating or initial

boundary of the various elements. Thus, XX(1) was set equal to input

position S(l) and is the closed end of the pipe. XX(Z) is the position of

the boundary between elements land 2. XX(I+1) is the boundary between

elements I and 1+1. The final boundary in the pipe is XX(N+1) which is

equal to S(32). Also in line 25, to save space, the variable, XF, was

set equal to XX(N+1) and represents the boundary between the pipe and

the dewar. Since XX is invariant with time it is single subscripted. X

varies with time as well as with the element being considered and, therefore, is double subscripted. Thus X(l, J) refers to the warm side of a given element at time J, and X(2, J) refers to the cold side at time J. To

conserve storage, the first subscript is always 1 or 2. This necessitates

a change when shifting from one element to the next. The value of X for

the previous element that was the cold side becomes the value for the

warm side o± the next element. This shift is accomplished in line 3 5 as

the program progresses from one element to the next. For element one,

the initial values o± X for all time increments are set equal to XX(1). So

control shifts from line 110 to line 690 whenever the warm side position

of the element is either at the end of the pipe or protruding into the dewar.

26
If the element is partly or wholly in the pipe at the beginning or end

of the interval, the calculation continues to line 120.

In line 120, to insure proper operation of the program, two heat

transfer parameters are set to zero.

The program now calculates heat transfer from the wall to the gas.

Lines 240 through 490 do this for elements wholly in the pipe. Lines 240

through 3 00 plus lines 500 through 2689 do it for elements partly in the

pipe and partly in the dewar.

The heat transfer calculation is necessarily iterative. From the

previous calculations on other elements, the warm end position of the element is known. Then from an approximate initial volume of the element

its approximate average position is found. Heat is transferred as if this

were its real average position. A more accurate density and volume are

calculated. From this new volume a new average position is found and

the heat transfer is recalculated. The process is repeated until the heat

transfer calculation stabilizes.

In line 240, the variable, KOUNT,has 1 added to it in order to

print out later the number of heat transfer iterations needed for convergence. KOUNT was initially set to zero in line 105 so that its value always

represents the number of heat transfer iterations for a given element and

time interval during the last cycle iteration.

In line 240, the parameter, DQ1, is set equal to AQ, the heat transfer from the previous iteration, for convergence test purposes. Note:

AQ and DQ are essentially the same.

In line 241, the maximum value of KOUNT for the given element,

KOUNT MAX(I), is found. KOUNT MAX(I) is set to zero in line 65 so that

the final KOUNT MAX(I) is for the final cycle.

27
Line 245 sets the position of the cold end of the given element

equal to the warm end position plus the length of the element. This position

varies slightly with each KOUNT because the heat transfer changes the

volume of the element. That is, VOL(2) is the element volume at the end

of the time interval. Dividing by the cross section area, AF, gives the

length. But VOL(2) is not known absolutely until after the heat transfer

and work terms are evaluated. Initially, an approximation for VOL(2) is

found in line 45. Successive calculations make it more accurate. During

the final heat transfer iteration of the final cycle the cold side position of

the given element at time J, X(2, J), is firmly established.

Line 250 determines the average position of the element at the end

of the time interval, J, by averaging the cold and warm end positions.

Line 255 establishes the gas velocity during the interval. It divides

the motion of the average position by the time interval.

Line 260 determines the average gas temperature during the interval. For an initial approximation, line 35 sets the instantaneous temperatures during a cycle to the gas temperature found in ELEMENT for the

given gas element. During iteration a more accurate value is found for

the temperature at the end of the interval, T(J), and, therefore, for the

average temperature, TAVG. The temperature determined at the end of

the previous interval establishes the temperature at the beginning of the

current interval, T(JLJD).

Lines 265 and 270 calculate the maximum overpressure occurring

at the end of the time interval when the average position is at XAVG(2).

Line 275 multiplies this maximum overpressure by the sine of the angular

position to obtain the instantaneous overpressure at the end of the interval.

The corresponding absolute pressure, PRES(2), is obtained by adding

dewar pressure to the overpressure.

28
Line 280 obtains the average absolute pressure, PRESS, and the

average density, ROE, during the interval by averaging conditions at the

beginning and end of the interval.

Line 300 transfers calculations to line 500 if the element protrudes

into the dewar at the beginning or at the end of the interval.

If the element is entirely in the pipe, line 305 calculates its average

volume, VOLM, by dividing the element mass, GM(I), by its average

density, ROE.

Line 405 calculates the average position, AVGX, of the element

during the interval by taking the average of the average positions at the

beginning and end ofthe interval.

Lines 410 through 420 calculate the wall temperature, TWI, corresponding to the average position of the element during the interval.

Line 425 calls the subroutine HEAT TR. The quantity ITEST FRF

is set to zero so that heat transfer will be calculated. Had it been set to

unity, only friction factor information would have been returned.

Sometimes the program is unstable and troubles occur in subroutine HEAT TR. If this occurs, DT is set to zero in the subroutine and

line 430 causes termination of the program.

Line 435 calculates the infinitesimal friction work, DFR, done

in the interval by a formula derived in Appendix C. Since the friction is

always a positive number, the absolute value is used.

Line 440 calculates the infinitesimal net work, DWK, after subtracting friction work, by a formula derived in Appendix B.

Line 455 calculates the internal energy change of the element, QW,

from the first law of thermodynamics by subtracting the net work, DWK,

29
from the heat input, DQ, returned from HEAT TR. The quantity AQ,

used later, is set equal to DQ.

The gas temperature at the end of the time interval, T(J), is calculated in line 460. To the temperature at the beginning, T(JLJD), is

added the temperature rise computed by dividing the internal energy

increase, QW, by the mass of gas in the element, GM(I), and the specific

heat, CVG.

Lines 465 and 47 calculate the element density and volume at the

end of the interval.

During the first iteration of a heat transfer calculation in an interval, it is not possible for equilibrium to have been established, so line

475 forces the program back to line 240 for at least one more iteration.

Sometimes on undue number of heat transfer iterations occur and it

is necessary to proceed without complete convergence. If there have been

more than 10 attempts, KOUNT will exceed 10, and line 480 sends the

program to line 700 to continue on with the program. This seldom happens.

If it does, the quantity, KOUNT MAX(I), calculated in line 241 exceeds 10

and shows in the print-out.

A test for heat transfer convergence is in line 485. If the heat

transfer computed is appreciably different from the previous computation,

the program returns to line 240 and tries again. To avoid division by zero

the quantity AQ, if zero, is changed to 0. 000001.

If the program gets to line 490,it has completed heat transfer and

work calculations for this time interval of this element, and the calculations

have been for an element entirely in the pipe. It now skips calculations

designed for elements partly or wholly in the dewar and control is shifted

to line 7 00.

30
If the element is partly in the dewar and partly in the pipe, either at

the beginning or the end of the interval, control is transferred to line 500 by

line 300. For heat transfer purposes it is necessary to determine the length

of the element touching the wall. Four quantities are defined, CI, C2, Wl,

and W2. These are the positions of the cold side at the interval beginning,

the cold side position at the interval end, and the warm side at the beginning and end. Depending on when and by how much the element protrudes

into the dewar, these quantities may be the sides of the element or the

ends of the pipe. Lines 500 through 530 determine these positions.

Lines 550 and 551 calculate the infinitesimal amount of kinetic

energy ejected into the dewar during this time interval.

Line 600 calculates the heat transfer area at the start, AREA1,

at the end, AREA2, and the average area, AREA.

Line 605 calculates the average position of the element for heat

transfer purposes.

Lines 610 through 620 determine the wall temperature, TWI, at

the average element position.

Line 625 determines heat transfer into the element by calling subroutine HEAT TR. The quantity ITEST FRF is set equal to zero so that

HEAT TR returns heat transfer as well as friction factor information.

Line 63 terminates the program if a calculation instability is

found in HEAT TR.

Line 635 calculates the infinitesimal friction work done during the

interval by a formula derived in Appendix C.

Line 639 calculates the overage volume of the element in the pipe.

Since it is the volume in the pipe rather than the entire element volume that

is desired, the quantities G(l) and G(2) are used for the mass of gas instead

of GM(I).

31
Line 640 calculates the net work done in the interval by a formula

derived in Appendix B.

Line 655 calculates the internal energy change in the element.

Line 670 calculates the gas volume of the given element that is in

the pipe, VOLIPZ, the mass of this gas, G(2), and the mass of gas of this

element in the dewar, GID2, at the end of the interval.

Lines 674 and 67 5 calculate the density of the fraction of the gas

element that is in the dewar at the end of the interval, ROID2. This

quantity is used in line 67 6 to determine the volume of the dewar fraction

of the gas element, VOLID2, and the total volume of the element, VOL(2).

Note: The calculations oflines 674through 676 possibly are not completely

rigorous, but they do not affect the results obtained using kinetic energies,

EKE.

Lines 68 through 687 calculate the temperature of that portion of

the gas element that is in the pipe at the interval end. Three routines are

necessary, one for flow out of the pipe, one for flow into it, and one special

case where there is no gas in the pipe at the end of the interval. It should

be noted that T(J) is always the temperature at time J of the gas fraction

in the pipe.

Line 688 computes the density of the gas fraction in the pipe at the

interval end.

In lines 688 through 689 criteria similar to those of lines 475

through 485 are used to test heat transfer convergence.

At line 2689, the program skips calculations used for elements that

are entirely in the dewar.

32
Lines 690 through 697 are for the case of the gas element entirely

in the dewar. No computations are really needed. Those that are made

are for the purpose of maintaining consistency of printout.

At line 7 00 all calculations have been performed on the gas element

for the current time interval. In line 701, the infinitesimal work, DWK,

just calculated,is added to the work term for the element so that at the

end of the cycle the term, WORK(I), will represent the total work done by

the element in the interval. Similar steps are taken to obtain the total

heat input to the element, Q(I), and the total friction in the element, FR(I}.

In lines 703 and 7 04 the dewar bound kinetic energy of the element, ZKE(I),

is computed.

If the convergence from cycle to cycle is not yet adequate, the

quantity IQ will be zero and detail printing is not desired so line 7 05 sends

the program to line 40003 for return to line 80 to advance to the next time

interval, or if a cycle has been completed, to line 723.

If only one cycle has been computed for the element, line 723

transfers control to line 40002 which returns control to line 60 to

start a second cycle.

Lines 724 through 750 concern cycle convergence, printing, and

program termination.

There are three criteria for cycle convergence. For the elements

that are wholly in the pipe, the heat received by the element must equal

the net work done by the element. For all elements at equilibrium, the

heat received by an element in a calculation cycle must equal the heat received in the previous calculation cycle. The work done in a calculation

cycle must equal that done in the previous one. These three criteria are

calculated in the form of three ratios, QLWOQ, DELQIOQI, and DELWKOWK.

They are defined as follows. QLWOQ is the difference between the heat and

33
work, divided by the heat input. DELQIOQI is the difference between this

cycle's and the previous cycle's heat, divided by this cycle's heat transfer.

DELWKOWK is the same applied to work done. In all three criteria the

absolute value is used.

If the heat transfer and the work done are both zero, the criterion

is satisfied and line 725 sends the program to line 750.

Line 726 avoids division by zero.

Line 730 and 73 5 calculate the three criteria.

If the values of any two of the three criteria are less than 0. 0001

in lines 738 through 740, the program is sent to line 750.

The quantity NJJL2 was calculated in line 6 and is equal to two less

than the input quantity NJJ, which limits the number of cycles that can be

calculated. If the number of cycles is as great as NJJL2, it will be desired

to print detailed information before proceeding to the next element. The

program is therefore sent to line 750 by line 742.

Line 7 50 adds 1 to the variable, IQ. Then during the next cycle

various details will be printed. If cycle convergency has occured twice,

the variable IQ increases to two and detailed printing occurs again. When

IQ reaches three, cycle convergency is undoubtedly satisfactory and line

60 causes the program to proceed to the next element.

Line 40002 is the end of the cycle or JJ iterative DO LOOP. Line

40000 ends the element DO LOOP. Most of the remainder of the subroutine concerns printing. There are however a few calculations remaining,

In lines 760 through 40001 the work done in all elements is summed

as is the heat transfer and the friction. There are now two types of

kinetic energy that have been calculated. In each element the infinitesimals

have been summed to find the total dewar bound kinetic energy, ZKE(I), in

34
element I. The kinetic energies in either the N'th or M'th element, ZKE(N)

or ZKE(M), could serve as an approximation of the kinetic energy ejected

into the dewar.

The kinetic energy has also been calculated for each time interval

in line 550 and summed up in line 78 as SUMEKE. The quantities ZKE(N),

ZKE(M), and SUMEKE are all printed and found to agree closely.

The purpose of the program is to compare the sum of work done

by all the elements with the kinetic energy ejected into the dewar. The

best comparison is a ratio of the kinetic energy to the work. Lines 785

and 787 compute three such ratios. The quantity EKEOSW is SUMEKE/

SUM WORK; ZKEMOSW is ZKE(M)/SUM WORK and ZKENOSW is ZKE(N)/

SUM WORK. Because it is based on the most reliable method the quantity

EKEOSW is used.

9. Detail of Subroutine READ DATA in Frequency Calculations

The subroutine READ DATA for frequency calculations is very

similar to the one for amplitude. The difference is that much of the input

needed for amplitude is unnecessary for frequency. The input quantities

have the same meaning in each case. Subroutine READ DATA for frequency

is given in Appendix I.

10. Detail of Subroutine ELEMENT in Frequency Calculations

The subroutine ELEMENT, when used with frequency calculations,

divides the duct into N elements. The length of each element is adjusted

so that it takes sound the same length of time to cross each. Thus, the

length of each element divided by its velocity of sound is the same, and the

quantity, DXXOVC, is the same for each element. The same quantities

are computed and printed for each element as they were for amplitude

calculations. Additionally the velocity of sound for each element, C(I),

is computed and printed.

35
Line 10 sets the temperature of the warm end boundary of element

1, TEMP(l), to the input value of the w~rm end of the pipe, TMP(l). It also sets the position of this boundary, XX(1), equal to the input position

of the warm end of the pipe, S(l). In addition line 10 computes the value

of NP1 = N+l, which is needed later.

Line 20 is an estimate made of the quantity DXXOVC. This estimate is the length of the whole pipe, S(32) - S(l), divided by the number

of elements, N, and also divided by the velocity of sound in the median

element, SQRTF(XK*R*TMP(16)).

Lines 24 through 40002, for each element in turn, calculate the

sonic velocity in an iterative routine, and determine the boundary positions

of the element. Then lines 47 through 49 test to see if the cold end of the

N'th element is in its proper position. If not an adjustment is made to

the quantity, DXXOVC, «-nd the program returns to line 24.

Thus, line 25 estimates the sonic velocity of the element, C(I),

using the warm side temperature TEMP(I). Line 28 calculates the position of the cold side. Lines 29 through 33 calculate the temperature of

the cold side. In line 37, the average temperature of the element is

computed. Then the tentative sonic velocity of the element, CCC, is

computed using the average element temperature. Line 38 tests to see

if the tentative and estimated velocities are nearly the same. If they are

not, line 39 re -estimates the sonic velocity to be the tentative value and

the program returns to line 27. If they are nearly the same, line 43

establishes the sonic velocity as being the tentative value and the program

proceeds to the next element.

At line 47 the positions of all the elements have been computed.

In line 47 a test is made to see if the cold side position of the N'th element

closely coincides with the position of the dewar end of the pipe. If this is

36
1

not sufficiently close, the quantity, DXXOVC, is adjusted in line 49 and

the program returns to line 24 repetitively until it is sufficiently close

after a number of iterations. When it is sufficiently close the program

goes to line 50 and a number of quantities are computed and printed in

a manner similar to that of subroutine ELEMENT in amplitude calculations.

In addition, three other quantities are calculated. In the foregoing iteration

the velocity of sound, C(I), was determined for each element. In line 55,

the parameter, Z(I), is calculated. This is simply the product of the sonic

velocity and the density of the gas in the element, RHO(I). Also computed

is its reciprocal, A(I). These two parameters are needed in pressure and

velocity calculations which determine frequency. The remainder of

ELEMENT is for print out purposes.

11. Detail of Program FREQENCY

The program computes the frequency of thermal oscillations from

the known temperature gradient in the pipe, using equations derived in

Appendix A. In the theory, a sudden overpressure is assumed in the

dewar which travels towards the closed end according to the velocity of

sound and the derived equations. In the program of this report the equivalent assumption is made that there is a sudden underpressure in the pipe.

Frequencies calculated by the two methods are identical; the usage is one

of convenience. In either case the initial gas velocities are zero.

The time required for sound to traverse one element, DXXOVC,

is one computation time. At the start of the program, with everything

at rest, the time, J, is considered to be one. After the passage of one time

interval, DXXOVC, J becomes two. At the time [ J = 2] a pulse has traversed from the dewar across element N. At the time [ J = l] the program

inserts initial values of overpressure, P, and velocity, U. At the time

[ J = 2] the only change has occurred in element N, and the program determines the proper values. After these first two preliminary calculations,

37
the program computes, at each time J, the overpressure and velocity

in each element. Because of computer storage limitations, the maximum

value of J is 102. Then, if desired, the pressures and velocities in the

N elements are printed for the values of J from 1 to 100. The values of

pressure and velocity in each element at times [ J = l] and [ J = 2] are

now set equal to the values at [J=10l] and [ J = 102], and the calculations

proceed for another hundred time intervals. On each occasion that J is

thus shifted, the quantity KOUNT is increased by 100 and the quantity J2

is always equal to J+KOUNT, with the result that J2 is the number of the

calculation time. Because of the nature of the equations used, a different

route through the program is used for times when J is even or odd.

When first devised, the frequency was determined by plotting the

closed end pressure versus time. Then the number of cycles per second

was determined by examination of the plot. In the program presented here,

the frequency is determined in the program itself. Because of the somewhat erratic nature of the overpressure, the routine used to calculate frequency is quite complicated.

In principle the method is as follows. During each calculation time

for each element the quantity, PMIN( I, JK), is determined. If the overpressure, P(I, J), at the time is less than the previously determined PMIN(I, JK),

then the value of PMIN(I, JK) is adjusted downward to this new value.

Simultaneously the quantity, TMIN(I, JK), is set equal to the calculation

time, J2, at which the overpressure went to a new minimum. The subscript, I, refers to the element number, and the subscript, JK, to the

cycle number that is in progress. Initially the cycle number and JK are

unity and the closed end overpressure is -1000. The closed end pressure

then increases to a maximum overpressure and decreases. Eventually

it becomes negative. At this time a new low is being approached for a

38
new cycle and JK is increased by one. When a new minimum pressure,

PMIN(I, JK), occurs at a new time, TMIN(I, JK), the process is repeated.

Without some safety provision, the cycle number, JK, would increase prematurely. This results from the fact that during a pressure

trend that is generally up or down, there may well be a number of temporary reversals. This reversal occurring just as the overpressure crossed zero would erroneously increase JK. The safety provision introduced

is the quantity, MAXCOUNT. At each change of JK, MAXCOUNT is set

to zero. During each calculation, if the overpressure is slightly over

zero, one is added to MAXCOUNT. Until MAXCOUNT is greater than

N, the quantity JK cannot change. By this time the danger of inadvertently

having the overpressure cross zero is past.

TMIN of the first element is used as a criterion of the interval of

one cycle. This quantity and others are computed for all the elements.

These are,

TMAX(I, JK

PMAX(I, JK

UMAX(I, JK

TMIN(I, JK)

PMIN(I, JK)

UMIN(I, JK)

The interval

, the calculation time of the maximum overpressure

, the maximum overpressure of cycle JK

, the maximum velocity in cycle JK

the calculation time of the minimum overpressure

the minimum overpressure of cycle JK

the minimum velocity in cycle JK

ZINTERVL, of one cycle is obtained by subtracting

a given value of TMIN from the next value of TMIN. The period, PERIOD,

is this quantity multiplied by the calculation time, DXXOVC. The reciprocal of any one period is the frequency, FREQ, and the average frequency,

AVG FREQ, is the average of the first nine frequencies calculated. The

average frequency is then used in the subroutine AMPLITUD. Program

FREQENCY is given in Appendix K.

39
Lines 5 and 10 call the subroutines READ DATA and ELEMENT.

Line 20 sets the initial value of KOUNT at - 100 so that when the

calculation starts and 100 is added, the value of KOUNT during the first

100 calculations will be zero. The quantity SAMBDA is the value of the

initial assumed under pressure in the pipe. GAMBDA is a small quantity

used as a criterion to determine MAXCOUNT.

Line 25 sets the values of several quantities which are used later.

Lines 30 and 40000 sets the pressure in all elements at times [ J = 1]

and [ J = 2] to - SAMBDA and the velocities to zero.

Lines 40000 and 3 5 set the initial values of quantities which may

change later.

The first cycle starts at its minimum; therefore, line 40 sets

PMIN equal to - SAMBDA. For other cycles it is set to + SAMBDA to

insure the low value calculated later is the proper one.

Line 50 sets the velocity and pressure in the N'th element and

second calculation time equal to initial values.

In line 60 the value of KOUNT is increased to zero before the first

series of a hundred calculations, and increased by 100 after each series.

Lines 65 through 141 are the main DO LOOP that calculates the

pressures and velocities. In line 65, if the calculation time is the 3rd,

or any other odd time, the program goes to line 101. If even, it goes to

line 121.

Lines 101 through 115 are the DO LOOP that calculates odd time

overpressure, velocity, and parameters associated with PMAX and PMIN.

The overpressures and velocities are calculated in lines 105 through 113

by formulas derived in Appendix A. Lines 2114 through 2129 determine if

40
the overpressures and velocities calculated are a maximum, a minimum

or neither for the cycle and set up a new cycle when MAXCOUNT is greater

than N and the closed end overpressure, P(l, J) is negative.

Lines 121 through 5141 perform the corresponding calculations for

even values of J.

At line 150 the overpressures and velocities have just been calculated by the DO LOOP ending in line 141 for a hundred values of J. In line

150, if KOUNT equals zero or is greater than 1000, detailed printing of

pressures and velocities is desired and carried out by lines 160 through

225.

In line 226, if not enough cycles have been completed, the program

jumps to line 229 where the overpressures and velocities of the 101st and

102nd calculations times are transferred to the 1st and 2nd time, and

control returns to line 60 for computation of another hundred values of J.

In line 226, if enough cycles have been completed, the program

continues to line 6227 for final calculations.

Line 6228 finds the number of calculation times, ZINTERVL(J), for

each cycle.

Line 6229 converts ZINTERVL(J) to cycle period, PERIOD(J), by

multiplying by the calculation time, DXXOVC.

Line 6230 takes the reciprocal of PERIOD(J) to find the frequency

of each cycle, FREQ(J).

Line 6231 adds all nine frequencies together.

Line 6232 determines the final answer needed for AMPLITUD, which

is the average frequency, AVGFREQ.

The remainder of the program is for print-out.

41
12. Results and Conclusions

Typical results from the frequency program are shown in figure

3. To obtain this plot a version of the program is used which does not

compute the frequency. Instead, the closed end overpressures are printed and then plotted. The plot shows the almost sinusoidal overpressure

obtained. It also shows the tendency for the calculated overpressure to

be erratic even though the general trend is clearly recognizable. Since

the temperature profile for the calculations is from figure 4 which is from

the work of Bannister (1965), it is possible to compare measured with

calculated frequencies.

The curve of figure 3 shows three oscillation cycles over a range

of 321 calculation points which followed 3984 previous calculation points.

Over the whole calculation, the period averaged 107.2 calculation intervals of 0. 000302 seconds each. From this the calculated frequency is 30.9

cps. From figure 3, the period is about 105 calculation intervals, or

0. 0317 seconds, and the frequency is 31.5 cps. From the program of

this report which averaged the first 10 cycles, the frequency is 31. 25 cps.

It should be noted that figure 3 is from a slightly different program

than the one described in detail in this report. In figure 3 an initial over2

pressure in the dewar of 10130 dynes/cm is assumed. In the described

2

program an initial underpressure of 1000 dynes/cm is assumed throughout the pipe. This difference has no effect on the calculated frequency.

The calculated frequency is needed in the determination of pressure amplitude. In all the amplitude determinations of this report the

temperature profile of figure 4 is used, so the frequency of each calculated point is taken as 31. 25 cps. The program is designed however to

accommodate any predescribed temperature profile by modifying input data.

42
% ^ *

2

ujo/ SduAp ' 3ynSS3UdU3A0

i
3

Q>

d

U

C

is o

"} +->u

3

id

<D ti %

U

• H<

<o

*> u

<: <u (j -t-> .«5 tn

•1-1

k» rt 5> « ^

k. o * ^ l-D * ^ <-MO

Uj >^ <o 03 ^

a; *F 3

<\J <o o

Q

to

<u H

*. ' 3dfllVd3dlN31

44
The amplitude calculations depend on heat transfer coefficients

and friction factors. Neither of these quantities is precisely known for

the type of fluid movement that exists in thermal oscillations. In McAdams

(1954) is information for fully developed flow. Unfortunately the flow in

ftiis case is far from fully developed. Furthermore, it is intended that a high value of friction factor is to be used, to allow for inlet turbulence.

Accordingly, the results are presented for two types of heat transfer coefficients. With the one type, heat transfer coefficients are obtained directly from the text book. With the other, they are multiplied by 0.9.

Two methods are used for presenting friction factor. With one

method, friction factors are assumed to be constant during a run regardless of the momentary Reynold*s number. Then a series of runs are

presented, each with a given constant friction factor. With the other

method, friction factors are calculated from a text book formula for each

momentary Reynold's number, but the calculated value is increased by a constant multiplier for each run. Then a series of runs is presented,

each with a constant friction factor multiplier.

Three figures are presented, using a variety of friction factors

and two coefficient multipliers, showing the effect of pipe diameter on

the amplitude of oscillations. Figure 5 uses unmodified text book heat

transfer coefficients and a series of four constant friction factors. Figure 6 uses unmodified text book heat transfer coefficients and a series

of five friction factor multipliers. Figure 7 uses text book heat transfer

coefficients times a multiplier of 0.9 and the same series of five friction

factor multipliers.

The general shape of all the curves is the same. At the small

diameter end, the oscillations tend towards zero amplitude. As the

diameter increases, the calculated amplitude increases to a maximum,

45
Helium Gas

Heat Transfer

Coefficient Multiplier = 10

Frequency = 3l.25cps

\ \-tt

\ \ \o \ \ \

"* V \ \o

Q \ Y

3

\ \ V

\ \ \

Q \o \

\ \°

\

"^ \

vv \o

\o

) \o->

o \ o \

\& \

\ x*""* \° \ \o \

V

Figure 5

4

INSIDE

8 10

DIAMETER , mm

Amplitude Versus Pipe Diameter Using Constant Friction Factors

and Text Book Heat Transfer Coefficients

46
to 4

<uc

:>»

Helium Gas

o Heat Transfer

->\ o Coefficient Mu 1 tiplier =1 .0

/^ ( 3 \ i Frequency = 31.25 cps

1 -\ 1°

I"

3

1 a

1°

lo

1 -%

1 c

1 *'

1 —\

1 "

1 "(S

\ n V

01

\ \

(J1 \

V"4\ \

\

V

5 \ \ \

2 4 6 8

INSIDE DIAMETER , mm

Figure 6. Amplitude Versus Pipe Diameter Using Modified Calculated

Friction Factors and Text Book Heat Transfer Coefficients

47

10
7

m

g

<M £

o

c

UJ

a:

CO

</) UJ

a:

o_

a:LU

>

o

Helium Gas

Heat Transfer

Coefficient Multiplier =0.9

Frequency = 3l.25cps

\-n

\o

\ ^*" \o

\-n

to

\o

\c

A l "

\ Y^

\ Q V

\ \-^

3^ \ \

^^

2 4 6 8

INSIDE DIAMETER , mm

Figure 7. Amplitude Versus Pipe Diameter Using Modified Calculated

Friction Factors and 0. 9 times Text Book Heat Transfer

Coefficients

.

48

10
and then decrease asymptotically towards zero at large diameters. As

might be expected, the higher friction factors result in lower peak overpressures in figures 5, 6, and 7. In figures 6 and 7 the effect of heat

transfer coefficients may be seen. In all cases, reducing the heat transfer coefficient multiplier from 1„ to 0.9 reduced the peak overpressure.

As an approximation, the reduction is about 30 percent.

It is of interest to compare the calculatjd overpressures with the

data from Bannister (1965). For a 0. 63 5 cm (1/4 inch) diameter pipe

2

the measured peak to peak pressure oscillation was 220, 000 dynes/cm

2

(3. 19 psi) which corresponds to a peak overpressure of 110, 000 dynes/cm . 2 On figure 5 the constant friction factor that produces 110, 000 dynes/cm

is a value of 0. 0065. On figure 6, by interpolation, the friction factor

multiplier, which produces this overpressure, is about 1. 57, and on figure 7, by extrapolation, tne value is 1.37.

The constant friction factor method is not favored. At small diameters, the Reynold's number will be smaller and the friction factor

larger than at large diameters. As a result, the maximum overpressure

occurs at too low a diameter as is the case on figure 5.

It is not usually expected that heat transfer coefficients are known

precisely. Therefore the effect of diameter on amplitude can be followed

equally well with a heat transfer coefficient multiplier of 1.0 or 0. 9. The

value of 1. from figure 6 is therefore used. From the one datum point

the friction factor multiplier of 1. 57 has been selected. From the figure

it is found that the maximum overpressure, for the given temperature

gradient, occurs at 0. 3175 cm (1/8 inch) and its value is about 480, 000

dynes/cm (7 psi), and double that value for peak to peak variation. At

0. 63 5 cm (1/4 inch) the value coincides with test data and at large diameters tends towards zero.

49
Considering that no allowance exists in the program for entrance. turbulence, and further considering that friction factor determination

cannot be precise, it is quite reasonable to have to multiply by 1. 37 or

1. 57 to correlate calculated oscillation overpressures with test data.

It is concluded that the program in its present form calculates

overpressure peak amplitudes that are approximately correct. Therefore,

the qualitative explanation of the phenomenon, given in this report, which

generally follows Rayleigh (1945), is correct. Further, the heating effect

which evaporates helium is the result of the operation of the work cycles

which expel kinetic energy and warmed gas into the liquid.

Thus the method of oscillation suppression, by inserting a knotted

string in the pipe is in effect increasing the coefficient of friction.

Similarly the insertion of glass wool or other packing does the same thing.

Clement and Gaffney (1954) found that a vacuum insulated line does not

oscillate. This case corresponds to a heat transfer coefficient multiplier

of about zero. Since a 0. 9 multiplier greatly reduces oscillations, it is

expected that a multiplier close to the value of zero would essentially

eliminate them. Furthermore, since the oscillations are a thermal cycle, a source of heat is needed to provide the cycle work. Complete insulation

eliminates the source. In the same reference it is found that closure of

the cold end suppresses oscillations. For this case there is no kinetic

energy ejected into the dewar. Any chance oscillations that start come to

equilibrium when the sum of all the work terms equal the kinetic energy.

This equilibrium occurs when the work terms are all zero, or the oscillations are non-existent.

Although the program does simulate frequency and amplitude of

thermal oscillations, it has a defect which should be improved upon. The

assumption that the maximum overpressure varies as a cosine function

from one end to the other of the pipe is not needed. The assumption that

50
the closed end pressure is sinusoidal with time is correct and should be

retained. Then using the gas laws and Newton's laws of motion the

boundaries of the elements can be calculated one by one from the closed

to the open end. Simultaneously the instantaneous pressure variation at

all of the boundaries can be obtained. Such a calculation would be an

improvement over the present one and is recommended for future considerations. Additional work is needed, particularly in experimental

phases such that more appropriate comparisons can be made with the

mathematical model.

Although calculations have been performed and comparisons

made for one particular case, the program is expressly designed to

accommodate different fluids, tube lengths and temperature gradients

along the tube.

51
13. References

Bannister, J. (1965), Air Reduction Co. , Murray Hill, N. J., Private

Communications

.

Clement, J. R. and J. Gaffney (19 54), Thermal Oscillations in Low

Temperature Apparatus, Advances in Cryogenic Engineering,

Vol. 1, p. 302, K. D. Timmerhaus, Ed., Plenum Press, Inc.,

New York, N. Y.

Johnson, V. J. (Oct. I960), A Compendium of the Properties of Materials at Low Temperature (Phase 1), Part 1, Properties of Fluids,

WADD Technical Report 60-56, Parti.

McAdams, W. H. (1954), Heat Transmission, Third Edition, McGrawHill Book Co., Inc., New York, N. Y.

Rayleigh (1945), Theory of Sound, Vol. II, Dover Publications, New

York, N. Y.

Rich, G. R. (1951), Hydraulic Transients, McGraw-Hill Book Co. , Inc.,

New York, N. Y.

52
14. Appendix A

To analyze the behavior of pressure pulses traversing a non

homogeneous medium it was necessary to make an assumption. The gas

temperature in the tube varies continuously from one end to the other.

If an attempt is made to analyze this situation, an infinite number of

reflected pressure pulses would be obtained, which cannot be handled in a satisfactory manner. It therefore was assumed that the behavior of the

continuous temperature gradient could be closely approximated by an

equivalent gradient composed of a large number of constant temperature

elements. As the number of these elements becomes sufficiently large,

the behavior of the continuous gradient and the discontinuous gradient

should become identical. Such an assumption is quite common in

analyses of physical situations. The assumption makes no requirement

of the length of each gas element other than the tacit one that each element

length should be small. As a convenience in analysis, the element lengths

were selected such that the time required for sound to traverse each gas

element would be the same. As a further convenience in analysis an even

number of elements was used.

In the derivations, all pressures were considered as absolute.

However, since the derived expressions showed pressure differences, it was often possible in the digital computer program to ignore the absolute

pressure level and to consider only over or underpressure. Of course,

in determining gas densities, the absolute pressure must be used. The

quantity, J, is used in analysis. When J is one the overpressure has just

occurred. After the length of time it takes sound to traverse one element,

J has increased to two. After each such time increment, J is increased by

one.

At time [J= l] , with reference to figure 8, the process is just

ready to start. In each element the overpressure and gas velocity are

53
Element)

I- I Z 4 5 G 7

Element boundaries

Time- , J~- odd

S0/7/C barriers

Time , Lf-ei/cn

Time , J~ odd

Time , J= a/en

Figure 8. Sonic Front Patterns At Odd And Even Calculation Times.

54
zero. In the dewar there is a small overpressure X. A moving sonic

front exists at the tube end between the dewar overpressure and the tube.

At time [ J = 2] the dewar overpressure has penetrated across the

gas element nearest the dewar, element [I =N] . Because of the discontinuity in gas properties between element [ I =N] and [ I = N - l] , a portion of the overpressure is transmitted into element [I = N - l] and

a portion is reflected through element [ I = N] . The values of the resulting gas velocities and pressures will be derived later. At time [ J = 2]

there is no longer a sonic front at the tube end but there now is one between elements [ I = N] and [ I = N - l] .

At time [ J = 3] the transmitted wave has traversed element

[ I = N - 1] and is at the border of elements [ I = N - l] and [ I = N - 2] . At the same time the reflected wave is at the tube end between the dewar

and element [ I = N] . At time [ J = 3] there is a sonic front between the

dewar and element [ I = N] , and between elements [ I = N - l] and [ N -2] .

At each succeeding time there will be a new pattern of sonic fronts.

After the transmitted pulses have reached the closed end it may be seen

that the pattern will be continuously repetitive. At each time [ J = odd number] there will be sonic fronts at the two tube ends and at the dewar end

of all even gas elements. (This was the reason for using an even number

of elements. ) At each time [ J = even number] , there will be no sonic

fronts at the tube ends, but there will be some at the dewar end of all odd

gas elements.

For a tube with N = 8 gas elements this situation is illustrated on

figure 8 .From figure 8 it may be seen that at any time [ J = odd] all gas

elements are grouped in pairs, with each pair surrounded by sonic fronts.

55
For any time [ J = even] most of the gas elements are grouped in such

pairs, but elements, [ I = l] and [ I = N] , are arranged independently.

In the following derivations the gas velocities and pressures are

derived in terms of gas velocities and pressures at previous times. The

first derivation is for element pairs surrounded by sonic fronts. Separate

derivations then follow for the cases, [ Time J = even, Element 1=1 and

I = N].

On figure 9 at time [ J -2] there are two regions, [ I] and [I - 1],

surrounded by sonic fronts. Later at time [ J] they are again in the same

condition. At time [ J - 2] two pressure pulses have just started towards

the center. At time [ J - 1] they have met at a density discontinuity, and

at time [ J] they are widely separated again. At the intermediate times

[t ] and [t ] the pressures, velocities and other gas functions 9 are

known on each side of the sonic fronts in terms of functions cp at times

[J -2, J - 1, and J] .

At time [t ] on the left side, from continuity

dM=(C

I-l. t- U

I-l, J-2)PI-1, J-2

Adt (1)

dM=(C

I-l,.

1

- U

.-I, J-1 )PI-1, J-l

Adt - (2 >

From momentum it follows that

(p

i-i, j.r p

i-i, j- 2

> = dM < u

i.i, j.r ui-i, j-2>/<Adt)

• < 3 '

Combining (1) and (3) results in

(P - P ) = V 1-1, J-l 1-1, J-2 ;

(C

I-1, t

x

~ U

I-1, J-2> P I-1, J-2

(U

I-1, J-l" UI-1, J-2>'

56

(4)
/<kgion _ZV

<Pt-,,T-2

Region X

0i,J-*

Time = J-2

T/me = tf

Time - t, +d£

7?wb = J-/

Time - t<

Time = t? fdt

0T-, 0.'*,* 77/ve = J~

Figure 9. Element Pair Surrounded by Sonic Fronts.

57
Combining (1) and (2) leads to

(C

I-l.t - U

I-l,J-2> P I-l,J-2 " (CI-l >tl

- U

I-1, J-l> "i-l. J-l • < 5 >

By algebraic manipulation of (5)

(C

I-1, f

U

I-1, J-2>

=

(6)

(u

i-i, j-r u

i-i, j-z) p i-i, j-i /(p i-i, j-r p i-i, j-2> •

Solving (4) for (C - U ) and multiplying by (6) 1~ 1 j £ .. x~ 1 1 J — £* produces

.2 = P

I-1, J-l ' PI-1, J-2 P I-1, J-l ^ 4 l"

I

" 1

' J " 2

"I-l. J-l " "l-l. J-2 'i-l, J-2

At time t on the right side, by the same operations

(P

I,J-2 " P

I, J-l^ = < C

I, tl

+UI, J-2> P

I, J-2 < U

I, J-l" UI, J-2>

(4)'

(P - P ) P ^1 *' J " 2

"

(P

I, J-2" PI, J-l) P

I, J-2

U

At time t on the left side

2

(P

I-1, J-l"

P

I-1, j)

=

(C

I-1, t 2

+U

I-l, J-1> P I-1, J-1< U

I-1, J" U

I-1, J-l)

(4)»

58
(C + U T1 ) 2

= Pl- 1 ' J -' " "•'' J ^i . (7)" H. t

2

l-l, J-i p j_ x -p j f ^

At time t on the right side

< p

i, j" Pi, j-i> < c

i, « 2

- u

i, j-i> \ j-i < u

i, j" ui, j-i» < 4>'"

P - P p (C -U \ 2 = _J±_i 1> J " 1

*> J (7)ih

^z *• J " x PI, J"\ J-l P

I, J-l '

'

It may be noted that where the velocity of the fluid is small compared with sonic velocity and where the density in a given element changes

only slightly with the passage of a pulse, the velocity of sound relative to

the wall is the usual expression for sonic velocity.

Since each sonic front during any time period is in a medium of

constant properties, each front is at a constant velocity. From this fact

and from continuity it follows that pressure and gas velocity between two

adjacent sonic fronts is the same on each side of the density discontinuity.

P

I-1, J-2

= P

I, J-2 * 8)

°I-1. J-2 " UI, J-2 . (9)

P = P

1-1, J I, J (10)

u

i-i, j- uW <">

For small pressure pulses the equations may be simplified. The

density in a given region will not change appreciably, and the gas velocities will be small compared with sonic velocities. After this

59
simplification and using (9) through (12), (4) becomes

P -P =Cp(U -U ) (12)

1-1, J-l I, J-2 1-1 M I-1 |U

I-1, J-l I, J-2' lx '

P

I, J-2 " P

I, J-l " C

I PI (UI, J-l - \ J-2> < 12>'

P

I-1, J-l " P

I, J

" C

!-l

P I-1

(U

I, J- UI-1, J-l' < 12) "

P

I, J

" P

I, J-l " CI h

(U

I, J

" U

I, J-l» • (12'

Adding (12) to (12)' produces the intermediate equation

P -P =C p(U -U J

1-1, J-l I, J-l 1-1,

H I-1 l I-l.J-1 I, J-2'

+T C p (U - U ) . I M

I

v I, J-l I, J-2;

Adding (12)" to (12)'" gives

P -P =Cp(U-U ) 1-1, J-l I, J-l 1-1 H I-1 V I, J 1-1, J-l 1

+ c

i

p

i (u

i, j

" u

i, j-i' •

Upon adding the intermediate equations

< u

i, j

- u

i, j-2' 2 \\ iXXl'

1 <13)

Upon subtracting the intermediate equations

1 \

/ 1 \

u

i. J-i ^c

< u

i, j + u

i, j-2> =

2

1 iV7 i x —

13>' V/ ic

.-i '.V

60
Similarly

(p

i, j

- p

i, j-1 = z rt vtt. \

J '\ • (14»

P

I J-l

(C

I-1 "l-l* + PI-1 J-l

(C

I "i*

,pw+ p

t j-z>

2 v^CmIi - (14''

Equations (13) and (14) determine the gas velocity and pressure

at any time [ J = odd] for any gas element [i = even] in terms of previously determined velocities and pressures. In each element [ (I - 1) = odd]

,

the velocities and pressures are identical to those in the elements [I =

even] because of (8), (9), (10), and (11).

Equations (13) and (14) also determine the gas velocity and pressure at any time [ J = even] for all gas elements [I = odd] except for the

first element [ I = 1] . In each element [ (I - 1) = even] , the velocities

and pressures are identical to those in the elements [I = odd] .

For all times [ J = even] the pressure in the element [ I = N] is

dewar pressure. The gas velocities are obtained with the aid of figure 10.

Using the principles of momentum and continuity as before

(P - P ) TT TT

V N, J-l N, J; N. J" N, J-I-P^^^C^^+U^ j^) '

For small pulses the density does not change appreciably from one time

to another. Also the velocity of the gas is small compared with sonic

velocity. So for small pulses for all times [ J = even] the gas velocity in

the element [ I = N] may be found from,

61
Element N

0>WjJ--/ dcwar Time t/V * odcJ

9*1,7-1 N,J

°N,t,

T/rne t,

9m,j-/ I 9/i/r

a%t

T/'me t, feft

0*,j T/'me J~= eifc/7

Figure 10. Sonic Front in Element N.

62
U -U _ P

N, J-l ' PN, J U

N, J

U

N, J-l" P NCN

• < l6 >

For all times [ J = even] , from continuity, the velocity in the

element [I = l] is zero. The gas pressures are obtained with the aid of

figure 11. Using the momentum and continuity principles again and

assuming small pressure pulses yields

P

l, J

" P

l, J-2 " " 2 Cl Pl U

l, J-l

=

" 2 C

l P

l U

2, J-l • < 17 >

Equations have now been derived which show gas velocity and

pressure in each of N gas elements at any calculation time J provided the

velocities and pressures are known for the previous two calculation times,

For the first calculation time [ J = 1] , all velocities and pressures are

zero. For the second calculation time [ J = 2] all velocities and pressures are zero except for the element [ I = N] . In that element the pressure is ^ above original and the velocity is obtained from equation (16).

Nomenclature for Appendix A

Symbols

2

P Pressure dynes/cm

C Sonic velocity relative to tube cm/sec

3

p Density g/ crn U Gas velocity (positive toward dewar) cm/sec

2

X Initial overpressure dynes /cm

M Mas s g

cp Sonic function, pressure, density,

velocity as determined by context.

2 A Area cm

63
Element /

0,,J~-2 T?'me> J~2 - euzn

Time t/

0*>" Time tt +dt

0./.j--/ h*' T/me J-/ =catt

'2. J-/

T/'me t>

2< J~f Time tz+dt

<Pi. j T/'me J~= ever?

Figure 11. Sonic Front in Element 1,

64
Nomenclature for Appendix A (continued)

Subscripts

I Gas element, I = 1 at closed end, N at open end

N Gas element closest to dewar

J Calculation time. At J = 1 process is just ready to start.

Time between J and J 4* 1 is time required for sound to

traverse one gas element.

Where double subscripts are used, the first refers to a

gas element, the second to the calculation time interval.

65
1 5. Appendix B

The infinitesimal net work done by a gas element in a short time

interval is,

dW = P dV - dF

where

dW = work done - - ergs

P = pressure - - dynes/cm

dV = volume change - - 3

cm

dF = friction work - - ergs • From the gas laws we have,

PV = MRT

dV

V

dT dP

T " P

PdV = MR dT-VdP

where

R = gas constant

T = gas temperature

M = gas mass

ergs/g-°K

°K

g •

The temperature rise may be expressed as,

dQ - dW dT = M C

v

where

dQ = heat received - - - ergs

C = specific heat - - - ergs/g-°K.

Substituting,

PdV = ^?- (dQ - dW) - VdP

v

66
PdV = -^- dQ - -^- (PdV -dF) - VdP

v v —R R

dQ + —dF - VdP

v v PdV = —

V

(v_i) (Y-1) VdP PdV = ^—J dQ + ^—^ dF - -^=-

y y y

dW = PdV - dF = i^-^-J dQ -

where,

/ VIl) dQ _VdP_dF

Y / Y Y

Y = specific heat ratio.

In Fortran notation this becomes,

DWK = C9*DQ-VOLM/XK*(P(2)-P(l))-DFR/XK

where,

DWK = net work - - - ergs

XK = specific heat ratio

C9 = (XK-1)/XK

DQ = heat received - - - ergs

3 VOLM = average volume - - - cm

2

P(2) = overpressure at end of interval - - - dynes/cm

2

P (1) = overpressure at start of interval - - - dynes/cm

DFR = friction work ergs .

67
16. Appendix C

When a fluid is moving in a pipe there will be a frictional force

between any gas element and the wall. This is expressed in equation

(6-2c) of McAdams (1954) as

T = £ v—

y

where

t = force per unit wall area - - - dynes /cm

f = Fanning friction factor

3

p = fluid density - - - g/cm

V = fluid velocity - - - cm/sec

From this, the frictional work overcome by an element with a

given area is,

- PV \ FdX^f ii~ ) A V dt

where

where

F = wall force - - - dynes

dX = distance moved - - - cm

2 A = wall contact area - - - cm

V = fluid velocity - - - cm/sec

dt = time increment - - - sec

In Fortran notation

DFR = ABSF(FRF/2. *AREA*ROE*VEL**3*DT)

DFR = friction work - - - ergs

2 AREA = element wall contact area - - - cm

3 ROE = average density - - - g/cm

VEL = fluid velocity - - - cm/sec

DT = time increment - - - sec

68
The absolute value is used because friction work is always a positive

number to be subtracted from the work done.

69
17. Appendix D. Program THM LOS

PROGRAM THMLOS

CALL READ DATA

CALL ELEMENT

CALL AMPLITUD

END THMLOS

70
18. Appendix E. Subroutine READ DATA for Amplitude

SUBROUTINE READ DATA

C FRICTION FACTOR (FRF) IS READ» BUT IT MAY BE CHANGED IN (HEAT TR)

C IF THE INPUT FRICTION FACTOR (FRF) IS ZER0»THE QUANTITY (IFRF) IS SET EQUAL TO

C ZERO IN SUBROUTINE (AMPLITUO) , AND THE SUBROUTINE (HEAT TR) COMPUTES THE HEAT

C TRANSFER COEFFICIENT ACCORDING TO A FORMULA IN MC ADAMS» MODIFIED By THE

C MULTIPLIER (NFR)

.

C

C IF THE INPUT FRICTION FACTOR (FRF) IS NON ZERO, THE QUANTITY (IFRF) IS SET EQUAL

C TO UNITY IN SUBROUTINE (AMPLITUD) AND THE INPUT VALUE OF THE FRICTION FACTOR

C IS USED REGAROLESS OF REYNOLDS NUMBER (REYNO)

.

COMMON/l/pSTART»DI AM, FRF, FACTOR, AVGFREQ,PM,N,NJJ,PI,AF,XK,CPG»M,

1 TMP(32)fS(32)»TMRVIS(l6) »VIS (16) ,TCONGAS ( 12) » CONGAS ( 12) ,ZNU<16),

1 GRPR<16) »CVG,RHODEW,TDEW,R,jD

C0MM0N/4/NFR

TYPE REAL NFR

5 READ800ltPSTART,DI AM, FRF, FACTOR, AVGFREQ,PM,N,NJJ. NFR, JD

8001 FORMAT(6F10.0»2I5»F5.0»I5)

10 READ 8000f (TMP(I) 1*1,32) $ READ 8000 » <S < I )

»

1=1 ,32)

15 READ 8000, (TMR VIS ( IN) IN=1 , 16) SREAD 8000, (VIS (IN) , IN=1 , 16)

20 READ8000HT CON GAS ( IN) , IN«1 , 12) $READ8000» (CON GAS ( IN) , IN*1 , 12)

25 rEAD 8000* (Z NU ( IN) , IN»1 , 16) SREAD 8000, (Gr Pr ( IN) , IN = 1 . 16)

8000 FORMATU6F5.0)

CONVERSION OF INPUT DATA TO UNITS USED IN CALCULATIONS

30 DO 40012 K*1.32 $S (K) =S (K) *2.54

40012 CONTINUE $DIAM=DIAM*2.54

35 DO 40011 Kxl,16 $vIS <K) *yIS <K) /I 000 000.

40011 CONTINUE

40 DO 40010 K=l,12 SCONGAS (K) .CONGAS (K) •10 000.

40010 CONTINUE

C VISCOCITY OF GAS»VlS»IS REAO IN MICROPORE AND CONVERTED TO POISE

C C0NDUCTIVITY»C0N GAS, IS READ IN MILLIWATTS/ (CM DEG K)AND CONVERTED TO

C ERG/(S CM DEG K) , TMR VIS, DEG K.GOES WITH VIS

C T CON GAS, DEG K, GOES WITH CON GAS, FOR THE GRASHOF

C PRANDTL PRODUCT.GR PR,VERSUS THE NUSSELT NUMBER, Z NU, THE LOG 10 OF THE

C QUANTITIES ARE READ IN WITH NO CONVERSION. THE LENGTH S IS READ IN

C INCHES AND CONVERTED TO CENTIMETERS. THE SAME IS TRUE FOR DIAM

C 45 WT=4,003SRU=8.317E07$R3RU/WT$PIs3.l4l5926536$AF=PI/4.«DIAM*«2

50 XK=5./3.$ CPG=R*XK/(XK-1.)$CVG«R/(XK-1.)

55 T DEWs4.3 $M=N*N/2 SRHODEw»PSTART/R/TDEW

COMPUTATIONAL VALUES ARE PRINTED EXCEPT VISCOCITY IS PRINTED IN MICROPOISE

60 PRINT 9310

9310 FORMAT (//)

65 PRINT 9000.PSTART, DIAM, FACTOR, AVGFREQ,PM,N,NJJ, NFR, JD

70 PRINT9001,(TMR VIS (K) ,K-1 , 16) SPRINT 9003» (VIS (K) ,K»1 , 16)

75 PRINT 9002»(T CON GAS (K) ,K»1 ,12) SPRINT 9004, (CON GAS (K) »K*1 , 12)

80 PRINT 9005, (Z NU (K) f Kal , 1 6 ) $ PRINT 9006» ((3 R p R <K)«K»ltl6)

85 PRINT 9009, (TMP(K) ,K=1, 32) SPRINT 9010, (S (K) ,K»1 ,32)

9000 FORMAT (5X6HPSTART,6X4HDIAM,4X6HFACT0R»3X7HAVGFREQ,

1 8X2HPM,9X1HN,7X3HNJJ.7X3HNFR»8X2HJD/

1 1XF10.2«F10.4,F10.6,F10.4,F10.0,2I10,F10*2,U0)

9001 FORMAT (//1X7HTMR VIS,2F7.2, 14F8.2)

9002 F0RMATQX9HT CON GAS»F5.2»F7.2»14F8,2>

9003 FORMAT (1X9HVIS (E06) . 1 (6PF5.2) , 1 (6pF7.2) » 14 (6pF8.2)

)

9004 FORMAT (1X6HCONGAS,F8.0»F7.0»14F8.0)

9005 F0RMAT(1X6HZ NU ,F7.3» 15F8.3)

71
9006 F0RMATUX6HGR PR ,F7.3» 15F8.3)

9009 FORMAT (1X6HTMP ,F7.3» 15F8.3)

9010 F0RMATUX6HS ,F7.3» 15F8.3)

99999 END READ DATA

C

72
19. Appendix F. Subroutine ELEMENT for Amplitude

SUBROUTINE ELEMENT

CONSTANT LENGTH OF EACH ELEMENT TYPE

CALCULATION OF ELEMENT LENGTHS AND TEMPERATURES

C0MM0N/l/PSTART»0IAM»FRFtFACTOR»AVGFREQ,PM,N.NJJ,PI,AF,XK f CPGtM,

1 TMP(32)»S(32)»TMRVIS(16) . VIS <16) tTCONGAS < 12) tC0NGAS(12) ,ZNU(16)

.

1 GRPR<16>iCVG»RH0DEWtTDEW»R,jD

C0MM0N/2/TEMP(97) ,XX(97> »RH0(96) »XMID(96> »GM(96> tTM6(96>

TYPE REAL NFR

5 TEMP(1)=TMP(1) S XX(1)«S<1) SNP1«N+1 $MP1=M*1

10 0XXe(S(32)-S(l) )/N

15 DO 40000 1*1, N $XX(I*1)«XX(I)*D XX

20 DO 40001 K*2,32 $IF <S (K) .GE.XX (1*1 ) ) Go TO 33

40001 CONTINUE $K*32

33 TEMP(I*l>sTMP<K)*<TMP<K)-TMP<K-l)>/<S<K)-S<K-l))»<XX<I*l)-s<K>>

35 TMG(I)«(TEMP(I)*TEMP(I*1))»0»5 $RHO ( I ) aPSTART/R/TMG ( I

)

40 X MID(I)«(XX(I)*XX(I*D)*0.5 $GM < I ) *RHo ( I

) #AF»DXX

40000 CONTINUE

45 DO 40002 I*NP1 »M$TEMP

(

I*J ) «TMG ( I

)

*T DEW

50 XX(X*1)"XX(I)*D XX $XMID{I)«(XX(I)*XX(I*1))»0.5 $RHO(I)=RHO DEW

55 GM(I)sRHO(I)*AF«D XX

40002 CONTINUE

60 PRINT 9310

9310 FORMAT (/)

65 PRINT 8038* (TEMP(K),Kel.NPl)

8038 FORMATOX 4HTEMR, 10F8.2»7F7.2)

70 PRINT 8039, (XX <K) »K=1 »NP1 ) 8039 FORMATOX 2HXX, 10F8.2»7F7.2)

75 PRINT 9206,(X MID (K) »K*1 ,N)

9206 F0RMAT(4X 5HX MID, 10F8.2»6F7.2)

80 PRINT 9205* (TMG(K),K»ltN)

9205 FORMAT (6X 3HTMG, 10F8.2»6F7,2)

85 PRINT 9000»(GM(K) ,K«1»N)

9000 FORMAT (7X2HGM,8E12. 4)

90 PRINT 9207i(RHO(K),K=l»N)

9207 FORMAT (6X3HRH0,8E12. 4)

9208 F0rMAt(8X1HC,8E12.4)

95 PRINT 9310 $PRINT8038» (TEMP (K) , KrNPl , MP1

)

100 PRINT 8039»(XX(K) ,K«NPl,MPl) $ PRINT 9206» <X MID (K) ,K=NP1 ,M)

105 PRINT 9205* (TMG(K) ,K«NP1»M) $ PRINt9000» <GM (K) »K=NP1 ,M)

110 PRINT 9207t(RHO(K)fK"NPl.M)

99999 END ELEMENT

C

73
20. Appendix G. Subroutine HEAT TR

SUBROUTINE HEAT TRUTESTFRF)

COMMON/ 1/PSTART»DI AM »FKF, FACTOR, AVGFREQ,PM,N,NJJ, PI , AF,XK,CPG»M,

1 TMp(32)fS(32)»TMRVlS(16) t VIS (16) »TC0NGAS(12) » CONGAS < 12) ,ZNU < 16)

,

1 GRPR<16) »CVG,RH00EW»TDEW,R,JD

C0MM0N/2/TEMP(97) ,XX(97) »RH0(96) »XMID<96) *GM(96) ,TMG(96)

C0MM0N/3/G(2) ,P<2) ,T( 360) ,PRES (2) ,X (2,360) ,XAVG(2) , EKE (360).

1 VOL( 2),R0( 2) ,FR(96) ,ZKE(96) tQ(96) ,W0RK(96) ,IFRF,

1 VEL,R0E,THETA,AREA,DQ,DT,TWI,TAVG,REYN0,HHT1,HHT2,HHT

C0MM0N/4/NFR

TYPE REAL NFR

1 IFdTESTFRF.EQ.DGO TO 386

2 DO 40014 K=2,12 $IF(T CON GAS(K).GE.T AVG ) GO TO 3

40014 CONTINUE $K=12

3 CON G=C0N GAS(K)*(CON GAS(K)-CON GAS(K-1 ))/(T CON GAS(K>-

1 T CON GAS(K-1 ))<MT AVG -J CON GaS(K))

4 DO 40013 K=2,16 $IF<TMR VlS(K).GE.T *VG ) GO TO 5

40013 CONTINUE $K=16

5 VISC=VIS(K)+(VIS(K)-VIS(K-1 ))/(TMR VlS(K)-TMR VIS(K-1 ))*

1 (T AVG -TMR VIS(K) ) 6 IF(VISC.EQ.0.)GO TO 400

8 REYNO=ABSF(DIAM»VEL*ROE/VISC)

10 IF(REYNO.LT.0.00l)REYNO=0.001

20 CONTINUE

22 IF(CONG)4Q0,400»25

25 PRANDTL=CPG *VISC/C0N G

26 IF(PRANDTL) 400*400*30

30 IF(TAVG)400, 400,31

31 ZNUSSELTso,023*REYNO#*0.8»PRANDTL**0.4

32 HHyl= ZNussELj*CON G/DIAM#FACtOr

33 THETA=TWI-T AVG SIFUHETA) 35,385,35

35 GP=CP6*ROE**2*980.66/T AVG*DIAM**3/VISC/C0NG*ABSF (THETA)

37 IF(GP) 410,40,40

40 IF(GP.EQ.0.)ZNUS=0.

41 IF(GP.EQ.O.)GO TO 380

42 GP L=LoGF(Gp)»0. 43429

43 DO 40016 K=2,16 $IF(GR PR(K).GE.GP L)GO TO 370

4OOI6 CONTINUE $K=16

370 ZNUL=ZNU(K)*(ZNU(K)-ZNU(K-1 ))/

1 <GR PR(K)-GR PR(K-1 ))*(GP L-GR PR(K))

375 ZNUS=10.«»ZNUL

380 HHT2=ZNUS«C0N G/DI AM*FACTOR $HHT =HHT1 + HHT2

385 DQ=HHT»AREA*THETA*DT

386 CONTINUE $IF( REYNO .LT. 0.001 ) REYNO»0.001

2386 IFdFRF.EQ.DGO TO 99999

387 IF(REYNO.LT.1000)FRF=16. /REYNO »NFR

388 IF(REYNO.LT.1000)GO TO 99999

389 FRF=0.01«NFR

390 FRF1=FRF

392 FRF=(l./(3.2* .43429*LOGF(REYNO*SQRTF(FRFD >*1.2> ) **2*NFR

393 IF(ABSF( (FRF1-FRF)/FRF) .GT.0.000DGO To 390

395 GO TO 99999

400 PRINT 9000.CON G, VISC,PRANDTL,TWI ,T AVq

9000 FORMAT (//I6X5HCON G» I6X4HVISC, 13X7HPRANDTL, 17X3HTWI , I5X5HT AVG/

1 1X5E20.5)

410 PRINT 9001»CON G, VISCGP, TWI ,T AVG

9001 FORMAT (//16X5HC0N G, 16X4HVISC, 18X2HGP , 17X3HTWI , 15X5HT AVG/

74
1 1X5E20.5)

420 DT=0.

99999 END HEAT TR

C

75
21. Appendix H. Subroutine AMPLITUD

COMPUT

12

3

40009

9004

9005

4

5

40004

6

10

15

25

30

31

32

35

40008

40

45

50

55

60

61

65

70

80

85

90

95

100

105

110

120

240

241

245

250

255

260

265

270

275

SUBROUTINE AMPLITUD

ATION ASSUMES PRESSURE VARIATION IS A COSINE FUNCTION

DIMENSION PX(96) KOUNTMAX (96) »JJMAX(96)

COMMON/1 /PSTART.DI AM* FRF FACTOR AVGFREQ f PM,N,Nj J, PI »AF«XK.CPG,M,

TMp(32)»S(32) 1TMRVISU6) .VIS(16) tTC0NSAS(12) »CONGAS(12) ,ZNU(16) , GRPRU6) »CVG»RHODEW»TDEW»R,jD

C0MM0N/2/TEMP(97) ,XX(97) »RH0<96) »XMI0(96> »GM<96) »TMG(96)

C0MM0N/3/G(2) »P(2)*T( 360) ,PRES(2> »X (2,360) t XAVG (2) *EKE (360)

.

VOL( 2),R0( 2)»FR(96) ,ZKE(96)»Q(96) ,W0RK(96)»IFRFf

VEL, ROE* THETA, AREA »DQ,UT,TWI»TAVG,REYNO,HHTl»HHT2»HHT

TYPE REAL NFR

IF(FRF.NE.0.)IFRF=1 $IF (FRF.EQ.0. ) IFRF=0

PRINT 9004

DO 40009 K*l»7 $REYNO=10.*»K $ITESTFRF*1 SCALL HEAT TR(ITESTFRF)

PRINT 9005»REYNO,FRF

CONTINUE

FORMAT (//10X,5X5HREYNo»7X3HFRF)

FORMAT(lOXf2El0.3)

PRINT 9006

PRINT 9006

DO 40004 J*l,360 $X (2» J) «XX d ) $EKE(J)aO.O

CONTINUE

NJJL2*NJJ-2

DT»JD/360./AVGFREQ $NL1«N-1 $ NP1=N*1 $ SUM W0RK=SUM Q=0.

SUMFRsGlD2aSUMEKEaO. $ROlD2eRHODEW

MP1*M*1 SXFaXX(N*f) $C8sPl«DlAM $C9=» (XK-1 . ) /XK $XL=XF-XX(1>

DO 40000 Ial.M

JJMAX(I)»0

G(2)*GM(I)

DO 40008 K-1O60 $ X (1 »K) -X (2»K) $T(K)*TMG(I)

CONTINUE

IQ=0

V0L(2)»GM(I)/RH0(I)

X(2*360)»X(1*360)*VOL(2)/AF

R0(2)=RHo(I) $XAVG(2)=(X(l,360)+X(2,360))*0.5$Q(I)=0.

DO 40002 JJsl,NJJ $IF(IQ.GT.2)G0 TO 40002

IF(JJ.GT.JJMAXd) ) JJMAX(I)«JJ

QI»Q(I) $KOUNTMAX(I)=0 JWORKl«WORK ( I

)

FR(I)= ZKE(I)=Q(I)=WORK(I)«0,

HO 40003 J=JD,360,JD $R2«J»Pl/180«$SlN2xSINF <R2)

IF(J.EQ.JD) JLJD=360 $IF ( J.GT. JD) JLJD=J-JD

RO( 1>*R0( 2) *VOL( l)*VoL( 2)

XAVG( l)aXAVG( 2) $R0lDl»R0ID2

PRES( 1)=PRES( 2) $G( 1 ) sG ( 2) $ P < 1

> "P < 2)

KOUNTsO $GID1=GID2

IF(XUtJ) ,GE.XF.AND.X(1»JLJD) .GE.XF)GO TO 690

DQ=AQ so. CONTINUE $KOUNT=KOUNT*l $Dq1»Aq

IF (KOUNT.GT.KOUNTMAX ( I ) ) KOUNTMAX ( I ) aKOUNT

X( 2, J )sX(ifJ )*VOL< 2J/AF

XAVG( 2)«(X( 2,J )*X(1,J ))*0.5

vEL=(XAVG( 2)-XAvG< 1 ) ) /Dj

TAVG=(T(JLJD)*T(J) )*0.5

IF (XAVG (2) ,LT.XF)PX(I)=PM*C0SF(PI/2.»(XAVG(2)-XX(1))/XL)

IF(XAVG( 2).GE.XF)PX(I)«0.

P( 2)sPX(I)#SIN2 $PRES( 2)-P< 2)*PSTART

76
280 PRESS=(PRES( 1>*PRES( 2) ) »0.5 $ROE= (RO ( 2)*RO( 1>)»0.5

300 IF<x(2»JLJD).GT.xF.0R.X(2»J).6T.xF)G0 TO 500

305 V0LM=6M(I)/R0E

400 AREAsC8*(X(2»J)-X(l»J)*X< 2 ,jLJD)-X(l,JLJD) )«0.5

404 G(2)aGM(I)

405 AVGXs(XAVG( 1)* XAVG( 2)) *0.5

410 DO 40005 K«»2»NP1 $IF (XX (K) .GE.AVGX) Go TO 420

40005 CONTINUE $K«NP1

420 TWI=TEMP(K)* (TEMP (K) -TEMP (K-l) ) / (XX <K) -XX (K-l )

) » (AVGX-XX <K) ) 425 ITESTFRFbO SCALL HEAT TRdTEST FRF)

430 IF(DT.EQ.0.)GO TO 99999

435 FR=ABSF(FR F/2. »AREA*R0E*VEL*»3»DT)

440 DWK «C9»DQ-V0LM/XK*(P( 2)-P( 1 )

) -OFR/XK

455 QWaDQ-DWK $AQS DQ

460 T(J)=T(JLJD)*QW/GM(I)/CVG

465 RO( 2)"PRES( 2)/R/T< J)

470 VOL( 2)«GM(I)/RO( 2)

475 IF(K0UNT.LT.2)G0 To 240

480 IF(KOUNT.GT.10.)GO TO 700

485 IF(AQ.EQ.0.)AQ*0.00000l $ IF (ABSF ( (AQ-DQ! ) /AQ) .GT.O.OOOl ) GO T0240

490 GO TO 700

500 IF(X( 2»J ) .GE.XF)C2=XF $IF(X( 2»J ) .LT.XF) C2«X < 2. J ) 510 IF(X(1 tJ ).GE.XF)W2*XF $IF(X(1 ij ) .LT.XF) W2=X ( 1 .J ) 520 lF(X(2tJLJ0).GE.XF)Cl=XF $IF (X (2f JLJD) .LT.XF) Cl=X <2» JLJD)

530 IF(X(1,JLJD).GE.XF)W1=XF $IF (X (1 , JLJD) •LT» X F) W1=X ( 1 , JLJD)

550 IF(VEL.GT.0.0)EKE(J)»0.5*(VEL*AF»ROE»DT)*VEL*VEL

551 IF(VEL.LE.0.0)EKE(J>*0.0

600 AREAls(Cl-Wl)#C8 $ AREA2« (C2-W2) #C8 $AREA= (AREAl*AREA2) #0.5

605 AVGX=0.25»(W2*W1+C2*C1)

610 DO 40006 KR-1.NL1 $K=N -KR $IF (XX (K) .LE. AVGX) GO TO 620

40006 CONTINUE

620 TWI=TEMP(K)*(TEMP(K)-TEMP(K*1))/(XX(K)-XX(K*1))*(AVGX-XX(K))

625 ITEST FRF*0 SCALL HEAT TR(ITESTFRF)

630 IF(DT.EQ.0.0)GO TO 99999

635 FR=ABSF(FR F/2. *AREA*R0E*VEL*«3*DT>

639 V0LM=(G(1)*G(2) )»0.5/ROE

640 DWK aC9*DQ-V0LM/XK«(P( 2)-P( l))-OFR/XK

655 QW=DQ-DWK $AQ*Dq

670 V0LIP2=AREA2»DIAM/4. $G (2) «RO (2) *V0LIP2 $GID2»GM ( I

) -G (2)

674 IF(GID2.GT.GIDl.AND.GID2.GT.0t)

1 R0ID2=(R0ID1«GID1*R0(1)*(GID2-GID1))/GID2

675 IF(GID2.LE.GIDl)R0ID2=R0IDl

676 V0LID2SGID2/R0ID2 SVOL <2) aVOLIP2*VOLlD2

680 IF(G( 2).LE.0.)T( J)»T DEW $ IF(G( 2).LE.0.)GO TO 687

681 IF(G(2)-G(1>)685.685»683

683 T<J)=(T(JLJD)«G(1)*TDEW»(G(2)-G(1) ) ) /G (2) *QW/CVG/G (2)

684 GO TO 687

685 T(J)=T(JLJD)*QW/CVG/G(2)

687 CONTINUE

688 RO( 2)*PRES( 2)/R/T( J) $IF (K0UNT.LT.2 ) GO TO 240

2688 IF(KOUNT.GT. 10)60 TO 700

689 IF(Aq.Eq.O.)Aq«0. 000001 $ IF (ABsF

(

(Aq-Dq1 ) /AQ) .GT. 0.0001 > GO T0240

2689 GO TO 700

690 DQ=DQl=AQsDWK=QW=D FRaQ, $REYNO«.001 $HHTlaHHT2aHHTaO.

691 ITEST FRFal SCALLHEAT TRdTEST FRF)

692 G( 2)a0. $T( J)aT DEW $r0( 2)«R0ID2 $VOL (2) aQM ( I) /RO ( 2)

77
693 X< 2»J )=X(1,J )*VOL< 2)/AF $XAV6< 2)=<X<2 ,J )*X<1,J ))*0.5

694 VEL=(XAV6( 2)-XAVG( 1))/DT

695 PRES( 2)=PSTART

696 TWI=0. $K0UNT=K0UNT*1

697 IF(KOUNT.GT.KOUNTMAXd) ) KOUNTMAX ( I ) =KOUNT

700 CONTINUE

701 WORK(I)=WoRK<I)*DWK $ U(I)bQ(I)*DQ $FR < I

) =FR (I ) *D FR

703 IF(VEL.GT.o.)DKE=0.5*(VEL»AF«ROE*DT)*VEL«VEL

704 IF(VEL.LE.0.)OKE=0. $ZKt < I

) sZKE (I ) *DKE

705 IF(IQ.LT.1)G0 To 40003

711 IF(J.EQ.90)PRINT 9001

9001 F0RMAT(///3XlHI,3XlHJ,iX2HJJ»lX5HK0UNT, 5X6H P<2)»4X7H R0(2)»

1 5X6H T(J) »5X3HTWI,8X3HVEL.9X2H0Q»8X3HDWK,8X3HDFR»3X8H X(2,J),

1 1X8H V0L(2)»5X6H G (2> /63X ,6X5HREYNO» 7X4HHHT1

.

1 7X4HHHT2.8X3HHHT.8X3HFRF)

712 IF<J.EQ.90.Or.J.EQ.180.Or.J,EQ.270.Or.J.EQ.360>

1 PRINT 9000tI.J»JJfKOUNT»P ( 2)»RO( 2)»T( J) ,TWI VEL»DQ,DWK.

1 DFR,X(2»J)»V0L(2) ,G(2)

.

REYNO.HHT1.

1 HHT2tHHT,FRF

9000 F0RMAT(lXI3»I4.l3,l6,Fll.l,E11.4,Fii,5,F8.3,4EH.3iFll.5»F9.4,

1 E11.4/63Xt6Ell.4)

720 IF(J.EQ.360)PRINT 9003.U < I ) , WORK ( I ) » ZKE ( I ) »FR(I)

1 ,QLWOQ»DELQIOQI»OELWKOWK

9003 F0RMAT(/12X4HQ(I) 8X7HWORKU). 9X6HZKE < I ) » 10X5HFR < I

)

1 ,10X5HQLWOQ»7X8HOELQIOQI»7X8HDELWKOWK

1 /1X4F15.3.3F15.5/)

40003 CONTINUE

723 IF(JJ.LT.2)G0 TO 40002

724 QI_WOQ=DELQIOQI=DELWKOWK = 0.0

725 IF(Q(I) .Eq.O.O.AND.WoHK(I) ,Eq.0.0)Go To 750

726 IF(Q(I) .EQ.0.G)Q< I) =0.000001 S IF (WORK (I) .EQ. 0.0) WORK < I ) =0.000001

730 QLWOQ=ABSF( (Q ( I

> -WORK ( I ) )/Q(I) ) $DELQlOQl=ABSF ( (Ql-Q ( I ) ) /Q ( I )

)

735 DELWKOWKsABSF< (WOrKI-WOrK ( I

)

)/WOrK(I) ) 738 IF(QLWOQ.LT. 0.0001. AND. DELQIOQI.LT. 0.0001)60 TO 750

739 IF(QLWOQ.LT.0.000l.AND.DELWKOWK,LT.0,000l)G0 TO 750

740 IF<DELQlOQI.LT. 0.0001. AND. OELWK0WK.LT.O.OOODG0 To 750

742 IF(JJ.GE.MJJL2)G0 TO 750

745 GO TO 40002

750 IQ=IQ*1

40002 CONTINUE

40000 CONTINUE

756 PRINT 9006

9006 FORMaT(IHI)

757 PRINT 9501

9501 F0RMAT(1X,4X1HI.1X8HK0UNTMAX»1X5HJJMAX,

1 5X5HPX(I) ,3X7HWORK(I) » 3X7HSUMW0RK.6X4H0 ( I )

.

1 6X4HSUMQ. 5 X5HFR(I) 5X5HSUMFR*4X6HZKE < I )

)

760 DO 40001 1=1, M $SUMq=SUMq*q < I ) $SUM WORK=SUMWORK*WORK ( I

)

761 IF (SUMWORK.EQ.0.)SUMWORK = 0. 000001

765 SUMFR=SUMFR+FR(I)

775 PRINT 9500»I»KOUNTMAX(I) JJMAX(I) »PX(I)

1 W RK(I),SUMW RK,Q(I) »SUMq,FR(I) »SUMFR, ZKE ( I

)

9500 FORMAT(1X,I5,I9.I6,8F10.0)

40001 CONTINUE

780 DO 40007 J=l,360 $SUMEKE=SUMEKE*EKE < J)

9310 FORMAT(//)

78
40007 CONTINUE

785 EKEOSW=SUMEKF7SUMW ORK $ZKENOSw=ZKE (N) /SUMwORK

787 ZKEMOSW=ZKE(M) /SUMWORK

PRINT 9008

9008 FORMAT (///1X»9X6HSUMEKE»9X6HEKE0SW,8X7HZKEN0SW,8X7HZKEM0SW)

PRINT 9007»SUMEKF.»EKEOSW»ZKENOSW»ZKEMOSW

9007 fORMAT(1X,4F15.4)

PRINT 9006

Do40010 J=JD«360»JD

PRINT 9009 f J,EKE(J)

9009 FORMAT (1X,6HEKE (J) , I3»EU,4)

40010 CONTINUE

99999 END AMPLITUD

C

79
22. Appendix I. Subroutine READ DATA for Frequency

SUBROUTINE READ DATA

COMBINATION OF STATEMENTS FOR FREQUENCY TYPE READ DATA

COMMON/l/N»TMP(32) ,S(3?) ,RH0DEW» TDEW.M.XK.R. TEMP (97) .XX(97) *PIi

1 C(96) «TMG(96) »XMI0(9fc) .RH0<96) .GM (9ft) . A (9ft) t Z (96) .AFtPSTARTt

1 DXXOVC

READ 8001tPSTART»DlAM,M

8001 FORMAT(2F10.C»1 110)

RFAD 8000* (TMP(I) tl = i»32) $READ 8000» <S ( I ) » 1 = 1 .32)

8000 FORMAT(16F5.0)

CONVERSION OF INPUT DATA TO UNITS USED IN CALCUl ATIONS

COMPUTATION USES CENTIMETERS RUT INPUT TS IN INCHES

DO 40012 K=1.32 $S(K)=S(K)*2.54

40012 CONTTNUE $DI AMbDI AM*2.54

WT=4.003SRU=8.317E07$R=RU/WT$PI=3.1415926S36$AF=PI/4,«DIAm**2

XK=5./3. $TDEW=4.3 $m=n*N/2 $RHODEWsPSTART/R/TDEW

COMPUTATIONAL VALUES ARE PRINTED

PRINT 9000»PSTART.DIAM,N

9000 FORMAT (1X13HPSTART»DI AM. N»F8.0»F8.4» 18)

PRINT QQo9« <TMP(K) »K=i ,3?) SPRINT 90l0 * <S (K) .K = l 32)

9009 F0RMAT(1X5HTMP 16F8.2)

9010 F0RMAT(1X5H S 16F8.2)

99999 END REAO DATA

C

80
23. Appendix J. Subroutine ELEMENT for Frequency

subroutine: element

common/1 /nttmp (32) »s(32) .rhodew* tdew»m, xk »r» temp (97) .xx (97) » pi

1 C(96) »TMG(96) »XMID(9ft) »RH0(96) »GM <9ft ) » A (9ft ) » Z (96) .AF.PSTART.

1 DXXOvC

CALCULATION OF ELEMENT LENGTHS AND TEMPERATURES

10 TEMP(1)=TMP(1) $ XX(1)=S(1) $NP1=N*1

?0 XX OV C»(S<ip)-S(l))/N/SQRTF(XK*R*TMp(16))

24 00 40002 I=1»N

25 C(I)»SORTF(XK*p*TFMP(I))

?7 CONTINUE

28 XX(I*1)=XX(I)*D XX OV c*c<n

29 00 40000 K=2.32 5 IF (S <K) .GE.XX ( I *1 ) > GO TO 33

40000 CONTINUE * K=3 ? 33 TFMP(I*1)=TMP(k)*<TMP<k>-TMP(k- 1 > > / (S (K) -S (k- 1 > >

<XX <I*1 ) -S («>>

17 TMG<T)=(TEMP<I)*TEMP(I+1) )/2. $ CCCsSQRTF (XK*R*TmG ( I ) ) 38 IF(ABqF(CCC-C<I) )/C(I)-. 001)43*43*39

39 C(I)aCCC $ GO TO 27

43 C(I)=CCC

4000? CONTINUE

47 IF(AB«!F(XX(N*1)-S(32) ) /S <32) -.0001 ) 50 »5o .49

49 O XX OV C=D XX OV C*(1.*\.*<(S(3?)-XX(N*1))/XX(N+1) )) SGOTO 24

50 DO 40017 I = 1«N $X MIO(I) = (XX(I)*XX(I + l) )/2,

SS RHO(T)=P START/R/TMG(I) SZ < I ) =C < I

) *RH0 ( I ) $A ( I

) =1 ,/Z < I

)

ftO GM(I)=RHO(I)*AF * (XX ( I + 1

) -XX ( I ) ) 40017 CONTINUE

ss print 9310

9310 FORMAT*/)

<?0 PRINT 9038. <TEMP(K) iK=liNPD

8038 format (ix 4HTEmp«ioP8.2»7F7,2)

9* PRINT 8039« (XX(K) »K*] »NPi) ! 8019 FORMATOX 2HXX. 10F8.2.7F7.2)

lnft PRINT R2o6*< x MID(K) »K=1.N)

92n6 FORMAK4X 5HX MID* 10F8. 2, 6F7.2)

loS PRINT 92 5» (TMG(K) iK=i»N)

92oS FORMAT (6X 3HTMO.

i

F8.2,6F7.2)

110 PRINT R000» <GM(K) .K=] *N)

9 000 F"ORMAT(7X2HGM»8El2.4)

120 PRINT 92 7» (RHO(K) iK=l»N) $ PRINT 92o8» <C (K) • Kel .N)

9207 FnRMAT(6X3HRH0»8El2.4)

9208 F0PMAT(8X1HC»8F12.4)

EVD FLFMENT

81
24. Appendix K. Program FREQENCY

PPOGpAM FREQENCY

C PROGRAM OF JUI Y 1967 WITH OMISSIONS AND PRINTOUT DELFTIOnS OF OCT 6 1967

COMMON/ 1/N.TMP (32) tS(32) ,RHODEW»TDEW»M, XK.R. TFMP (97) »XX(97) »PI»

1 C(96) «TMG(96) »XMI0(9ft) ,RH0(96) .Gm (96) A (96) .Z (96) .AF.PSTART.

1 DXXOvC

COMMON/2/P(96,]02) »U<96»102) .PMAX(96*12) .UMAX (96. 1 2) .UMlN (96. 12)

.

1 TMIN(96»12) .TMAX(96H2) dmIN(96«1?)

l.ZlNTFRVLd?) »PFRI0D(12) »FREQ(12)

5 call read data

10 call element

15 Print 9310

?0 kounfs-100 $ sambda^ooq. $gambda=0 . 1*sambda

?5 n.j=1^2 snj less 1=nj-1 $ nj less 2=nj-2$nl.1=n-1 $n plus 1=n*1

30 DO 4o000J=l*2 $DO 40000 I=1*N $ P d » j) =-SAMBDA $ U(I.J)=0.

40000 CONTINUE 5 DO 40001 JKsl.10 $ DO 40001 I=1»N

^«; p max<i»jk)=u max(i»jk)=ii min ( i » jk ) =0 .$t min(i.jk)*t max(i.jk)=1.

40 ifuk.fq.1jp min(i.jk)=-samboa $if ( jk .gt, 1 ) p min ( i . jk) =sambda

40001 continue $ jk=i $ max coiinT=o

calculation starts

50 u(n.?)=(-1)*sambda*a(n) % p<n.2)=0.

60 kount=koumt*l00 s sum frfq=0.

65 do 141 j=3»nj *j2=j+k0unt $if ( < j/2) »2-j) 101 » 121 » 101

computation of odd time pressure and velocity

101 DO 115 »I=2»N,2

105 P(L D=P(I«J-2)+2. <1'(U(I-i.J-l)-U(I»J-l5 )/(A(I-i)*A(IJ)

107 P(T-i *J)=P(I»J)

111 U<i»J)=-U(I»J-2)*2.<MlJ<I-T • J-l) *A (I) U(I» J-1 )*A ( i-l) )

1 /(A(T-1 )*A(I) ) 113 Ud-i .J)=U(I.J)

21J4 IF(P(I.J) .GT.P MAX(I»JK))T MAXdtJKjsJ?

2115 iF(Pd.J) .GT.P MAX(I,JK))P MAX ( I , JK) =P ( I . J)

2116 lF(Ud,J) .GT.U MAXd.JKMU MAX ( I , JK) =U ( I , j)

2117 iF(Pd.J) .I.T.P MIM(I,JK))T MIN(I.JK)=J2

2iTfl iF(Pd.J) .LT.P MIN(I»JK))P MIN(I.JK)=P(I.J)

2119 IF(Ud.J) .LT.U MlMd.JKJlU MTN d . JK) sU ( I . J)

2120 iF(Pd-l.J) .GT.P MAX( I-]tJK))T MAX ( I~l . JK ) =J2

2121 iF(Pd-l.J) .GT.P MAXd -1»JK))P MAX ( 1-1 JK) =P ( 1-1 . J)

2122 lF(UfI-l .J). GT.U MAX ( 1-1 » JK ) ) U MAX ( I"l . JK ) =U ( 1-1 . J)

2123 iF(Pd-l.J) .LT.P MINd-l.JK))T MIN ( 1-1 . JK) =J2

21?4 IF(P(I-1.J) .LT.P MlNd-l.JK) )P MIN <I-1 . JK) =P (1-1 , J)

2125 IF(Ud-l.J) .LT.U MIN H-l • JK) JU MIN ( t-1 . JK) =U ( 1-1 , J)

21P6 IF(I.GT.?)G0 TO 115

2127 IF(P(I.J) .GT.GAmBOA)mAX COUNT=MAx COUNT+1

212B IF(MAX COUNT. Gf.N. AND. P(1»J) .LT.O. )JK=JK*1

2129 IF(M*X COUNT. GT.N. AND. P(l »J) .LT.O. ) MAX CDUNT=0

115 CONTINUE

117 GO TO 141

COMPUTATION OF EVEN TIME PRESSURE AND VELOCITY

121 Ud.|)=0. SP<ltJ)=Pd»J-2)-2.*U(2»J-i)/A(l)

3122 IF(P(1,J) .GT.P MAX(1,JK))T MAX(1,JK)=J2

31P3 IF(Pd.J) .GT.P MAX(1,JK))P MAX ( 1 , JK) =P d » J)

3124 IF(U(1 »J) .GT.U MAX(1,JK))U MAX ( 1 , JK) =U (1 > J)

31?5 iF(PlliJ) .LT.P MIN(ltJK))T MIN(1,JK)=J2

31?6 IF(P(1.J) .LT.P MIM(1,JK))P MTN(1.JK)=P(1.J)

3l?7 TF(Ud,J) .LT.U MIN«1,JK))U MTN ( 1 JK ) =U ( j , j)

31?B IF (P(l . J) .GT.GAMBOA)MAX COUNT=MAX COUNT*)

82
3129 IF(MaX COUNT. GT.N, AND. P<1 »J) .LT.O. )JK=JK+1

3130 IF(MAX COUNT. GT. M. AND. P(nJ) .LT.O. ) MAX COUNT=0

1?5 00 134 1=3. N|_] »2

1?7 P(T»J)=P(I»J-2)*2.»( U(I-ltJ-l)-U(I,J-l) )/(A(I-l)*A(I) ) 129 P{I-l.J>=pn»J)

131 U(T.J)=-U(I»J-2)*2.*(U(I-i tJ-i)*A{I)*U(I,j-i)*A<I-i>)

! /(A(T-i )*A(D ) 132 U(I-HJ)=Ud'J)

4133 IF(P(I.J) .GT.P MAX<I,JK))T MAX(I.JK)=J2

4P4 IF(P(T«J) .GT.P MAX(I.JK))P MAX ( I , JK) =P ( I , J)

4135 IF(UU.J) .GT.U MAXlltJKMU MAX ( I , JK ) oU ( I , j)

4136 IF(P(I»J).LT.P MZN(X«JK))T MIN(I,JK)=J2

4137 IF(P«T*J) .LT-P MlN(ItJK))P MTN < I . JK ) =P ( I , J)

41T8 IF<U(I,J) .LT.U MIN(I»JK))U MIN ( I » JK) =U ( I . J)

4139 IF(P(I-1.J) .GT.P MAX( I-],JK))T MAX ( 1-1 JK) =J2

4140 IF(P(I-1.J) .GT.P MAX(I -1 . JK) ) P MAX ( 1-1 , JK) =p ( 1-1 , J)

4141 lF(U(i-l »J) .GT.U MAX<I-1 tJK) )U MAX ( 1-1 JK ) =U ( 1-1 »J)

4142 IF(P(I-1»J) .LT.P MIN(I-I.JK) )T MIN < 1-1 . JK) =J2

4143 IF(P(I-1'J).LT.P MIN(I-1,JK) )P MIN ( T-l . JK) =P ( 1-1 1 J)

4144 lF(IJ(I-i»J> .LT.U MIN(I-I.JK) )U MTN ( T-l , JK ) =U < 1-1 » J)

134 CONTINUE

135 P(N. |)rO. $U(N»J)=U(N»J-1) *(P<N.J-1)-P(N.J) )*A(N)

5136 IFfP(N.j) .GT.P MA)UN,JK>)T MAX(N,JK)=J2

5137 iF(PfN.J) .GT.P MAX(N,JKMP MAX <N, JK) =P (N, J)

5l3fl IF(U(N,J) .GT.U MAX(N,JK))U MAX (N, JK) =U (N, J)

5n<) lF(P(NtJ) .LT.P MIN(N,JK))T MIN(N,JK)=J2

5l4f) IF(P<N,J) .LT.P MIN(N.JK))P MIN (N, JK ) =P <N, J)

5141 IF(U(N.J) .LT.U MIN(N,JK))U MIN <N, JK ) =U (N» J)

141 CONTINUE

150 IF(K0UNT.GT. . AnO.KOUNT.LT. 1000 > GO TO 226

16n 00 2?5 J=1,NJ LESS 2 $ J2=J+K0UNT

220 IF(N.LF.16)PRIMT 9220 J2. <P ( I « J)

.

1=1 »N)

221 IF(n.lE.16)PRInT 9221

»

(U(ItJ) .1=1. N)

222 IF(N.GT.16)PRINT 9220 » J2, (P <I J)

,

1=1 » 16)

2?4 TF(N.GT.t6)PRTNT 9222» (P ( 1 • J)

»

I=17»N)

2?3 lF(N.GT.l6)PFriNT 9221 » <U ( I . J) . \ = \ » 16)

5223 T.F(N.GT.l6)PRlMT 9221 <U < I • J)

,

I=\7 ,n)

225 CONTINUE

226 IFUK.LT. 11)60 TO 229

6227 DO 4n003 J=l»9

6228 ZINTFRVL(J)=T MlN(l»J*l)-T MIN(l.J)

62?9 PFRIOO(J)=ZINTFR\/1 <J)*l) XX OV C 6230 FRFQ(j)=i./PERIOD(J)

6231 SUM fPFQ = SUM FREO-t-FRFO(J)

40003 CONTTNUE

6232 AVG FRFQ=SUM FrEQ/9.

6233 OMEGA=2.*Pr*AVG EPEG

6234 PRINT 93to

P310 FORMAT (//)

6235 DO 4«i005 J=l»l0

6236 PRINT 9300»J »(T MIN ( I , J) 1=]., 16 ) , <P MIN ( I .J) 1=1 » 16 ),

1 (IJ MTN(I, J) .1 = 1,16 ) , <T MAXU,J) »I = 1»16 ).

I <P MAX(I, J) ,1=1*16 ),<U MAX(T,J) ,I=1»16 ) 9300 FORMATM07H0MAXTMUM AND MTNIMUM PRESSURES AND VELOCITIES AND CORRF

1SPONDTNG CALCULATION TIMES IN ELEMENTS 1 THROUGH 16.

11X12HTHTS IS THE I2.9H TH CYCLE

83
11X5HT MINIGFS.O./IXFHP MIM16F8.0 » /1X5HU MIM16F8.0./

1 1X5HT MAX16F8.o«/lX5MP MAX16F8 .0 • /l X.5HIJ MAXlfFfi.O)

400nS CHNTTNUF

6?39 PRINT 930?. U»J= 1»9) <ZlNTEPVL<J) »J=1 *9) (PFRIOD(J) »J = 1»9) »

1 tPRFn(.J) »J=1»9) ,AV<5 FRFQ.OMFGA *D XX OV C 930? FORmAT(131h6ZInTERVl IS TMF nUmRER OF CA|CUiATIon UnITS In OnF CYC

1LF. PERIOD IS THE TIME LENGTH OF ONF CYCLE AMD THE RECIPROCAL OF

1THF PFRTOn /126H I*; THE FREQUENCY. THE TOP LINE GIVES THE CYCLE Nil 1MRFR FOLLOWING THE INITIAL PULSE. AVG FQFO IS THE AVERAGE OF ALL F 1RF0UEMCIES /

UXftHCYCl.E 9112,/lXRHZlwTEPVL 9Fl2. /1X8HPFRI0D 9F12.5*/

1 lX8HrRFQ 9F12.5»/9H0AVG FRFQ lFip.e; *SX5H0MEf!A F12.5*

1 3X9H0 XX OV C»E12.3)

6240 IF(JK.GE.11)G0 TO 99Q99

??9 00 2M 1 = 1.

N

7?T0 P(I»1 )aP(I.NJ LESS 1 ) SU < I 1 1 ) =U ( 1 1 NJ LESS 1)

7231 P(I»?)=P(I»NJ) $U(I«?)=U(I.NJ)

231 CONTINUE

232 GO Tn 60

92?n F0RMaT(1X.I4»1HP .16F8.0)

R2?1 FORMAT(5X.iHU ,i6F8.j)

92?? FORMAT rSX»lHP .16F8.C)

999Q9 EmO FrEOEnCY

C

g4 GPO 842 -367
NBS TECHNICAL PUBLICATIONS

p

o:

r,

PERIODICALS

JOURNAL OF RESEARCH reports National

Bureau of Standards research and development in physics, mathematics, chemistry, and engineering.

Comprehensive scientific papers give complete details of the work, including laboratory data, experimental

rocedures, and theoretical and mathematical analyses. Illustrated with photographs, drawings, and

charts.

Published in three sections, available separately:

Physics and Chemistry

apers of interest primarily to scientists working in these fields. This section covers a broad range of hysical and chemical research, with major emphasis

n standards of physical measurement, fundamental

constants, and properties of matter. Issued six times

a year. Annual subscription: Domestic, 35.00; for- eign, 36.00*.

NONPERIODICALS

Applied Mathematics Series.

tables, manuals, and studies. Mathematical

• Mathematical Sciences

Studies and compilations designed mainly for the

mathematician and theoretical physicist. Topics in mathematical statistics, theory of experiment design,

numerical analysis, theoretical physics and chemistry, logical design and programming of computers

and computer systems. Short numerical tables. Issued quarterly. Annual subscription: Domestic,

32.25; foreign, 32.75*.

• Engineering and Instrumentation

Reporting results of interest chiefly to the engineer

and the applied scientist. This section includes many

of the new developments in instrumentation resulting

from the Bureau's work in physical measurement,

data processing, and development of test methods.

It will also cover some of the work in acoustics,

applied mechanics, building research, and cryogenic

engineering. Issued quarterly. Annual subscription:

)omestic, 32.75; foreign, 33.50*. TECHNICAL NEWS BULLETIN

The best single source of information concerning

the Bureau's research, developmental, cooperative

and publication activities, this monthly publication

is designed for the industry-oriented individual whose

daily work involves intimate contact with science

and technologyfor engineers, chemists, physicists, research managers, product-development managers, and

company executives. Annual subscription: Domestic,

31.50; foreign, 32.25*.

•Difference in price is due to extra cost of foreign mailing.

Building Science Series. Research results, test methods, and performance criteria of building materials, components, systems, and structures. Handbooks. Recommended codes of engineering

and industrial practice (including safety codes) developed in cooperation with interested industries,

professional organizations, and regulatory bodies.

Special Publications. Proceedings of NBS conferences, bibliographies, annual reports, wall charts,

pamphlets, etc. Monographs. Major contributions to the technical literature on various subjects related to the Bureau's scientific and technical activities.

National Standard Reference Data Series.

NSRDS provides quantitative data on the physical

and chemical properties of materials, compiled from

the world's literature and critically evaluated.

Product Standards. Provide requirements for

sizes, types, quality and methods for testing various

industrial products. These standards are developed

cooperatively with interested Government and in- dustry groups and provide the basis for common

understanding of product characteristics for both

buyers and sellers. Their use is voluntary.

Technical Notes. This series consists of communications and reports (covering both other agency

and NBS-sponsored work) of limited or transitory

interest. CLEARINGHOUSE

The Clearinghouse for Federal Scientific and

Technical Information, operated by NBS, supplies

unclassified information related to Governmentgenerated science and technology in defense, space, atomic energy, and other national programs. For

further information on Clearinghouse services, write:

Clearinghouse

U.S. Department of Commerce

Springfield, Virginia 22151

Order NBS publications from:

Superintendent of Documents

Government Printing Office Washington, D.C. 20402
U.S. DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20230

POSTAGE AND FEES PAID

U.S. DEPARTMENT OF COMMERCI

OFFICIAL BUSINESS