
Citation 
 Permanent Link:
 http://digital.auraria.edu/AA00005830/00001
Material Information
 Title:
 Computer simulation of shock tubes
 Creator:
 Burris, John Randall
 Publication Date:
 1986
 Language:
 English
 Physical Description:
 vi, 65 leaves : illustrations ; 29 cm
Subjects
 Subjects / Keywords:
 Shock tubes ( lcsh )
Computer simulation ( lcsh ) Computer simulation ( fast ) Shock tubes ( fast )
 Genre:
 bibliography ( marcgt )
theses ( marcgt ) nonfiction ( marcgt )
Notes
 Bibliography:
 Includes bibliographical references (leaf 42).
 General Note:
 Submitted in partial fulfillment of the requirements for the degree of Master of Science, Department of Mechanical Engineering, 1986.
 Statement of Responsibility:
 by John Randall Burris.
Record Information
 Source Institution:
 University of Florida
 Rights Management:
 All applicable rights reserved by the source institution and holding location.
 Resource Identifier:
 17800848 ( OCLC )
ocm17800848
 Classification:
 LD1190.E55 1986m .B87 ( lcc )

Downloads 
This item has the following downloads:

Full Text 
COMPUTER SIMULATION OF SHOCK TUBES
by
John Randall Burris
B. S ., University of Colorado, 1983
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
of the requirements for the degree of
Master of Science
Department of Mechanical Engineering
1986
This thesis for the Master of Science degree by
John Randall Burris
has been approved for the
Department of
Mechanical Engineering
by
Date
i(kc
iii
Burris, John Randall (M.S., Mechanical Engineering)
Computer Simulation of Shock Tubes
Thesis directed by Professor William Clohessy
The use of computational analysis to model the
flow in straight shock tubes has been made in the past,
but variable diameter shock tubes, and nonideal shock
tubes have not been considered in any detail. This study
modifies an existing onedimensional, Eulerian, finite
difference code for a straight shock tube, to match data
obtained from an expanding diametershock tube and a
nonideal shock tube. The finite difference scheme em
ployed to solve the shock tube equation is called the
Godunov method, which integrates the Euler equations for
each cell using Green's theorm. A description of the non
ideal precursor wave is also presented. Results obtained
from the numerical computation agree well with the exper
imental data.
ACKNOWLEDGMENTS
I would like to take this opportunity to express my
appreciation to Roy Heyman, for his valuable advice and
guidance which made this project possible. And to the
Martin Marietta Denver Aerospace company for it's finan
cial support of this effort.
Special thanks are also due to my committee chair
man Dr. William Clohessy for his guidance and encourage
ment, and to Dr. Ken Ortega who has given me the inspira
tion to continue learning. I am indebted to Grady Romine
for his understanding and support through this effort.
Finally I wish to thank Sue, Beau, and Adrienne
for their endless love and patience, and for making
everything in life so worthwhile.
V
CONTENTS
CHAPTER
I. INTRODUCTION............................. 1
II. PHYSICAL MODEL........................... 4
Computer Code......................... 9
Expanding Nozzle Shock Tube.......... 15
Precursor Shock'Tube................. 18
III. GODUNOV METHOD.......................... 24
IV. RESULTS................................. 34
REFERENCES................................. 42
APPENDIX
A. Computer Listing and Subroutines.... 43
B. Expanding Nozzle Code.......... 54
C. Precursor Code.......................... 60
vi
FIGURES
Figure
2.1. Initial conditions in a pressure
driven shock tube...................... 4
2.2. Flow in a shock tube after the
diaphram is broken..................... 5
2.3. Cell boundary......................... 12
2.4. Thunder tube facility................. 16
2.5. Computer model of expanding cell... 17
2.6. Precursor wave development......... 19
2.7. Detailed precursor wave............... 21
2.8. Model for precursor shock tube..... 22
2.9. Configuration for precursor
shock tube............................ 23
3.1. One dimensional model for
computation........................ 2 6
3.2. Diagram of solution................... 26
4.1. Experimental data from the
Thunder tube.......................... 35
4.2. Experimental data from the
nonideal shock tube.................. 36
4.3. Computer generated data for
the expanding shock tube.............. 37
4.4. Computer generated data for
the nonideal shock tube.............. 38
4.5. Figure 4.3 imposed upon figure 4.1. 40
4.6. Figure 4.4 imposed upon figure 4.2. 41
CHAPTER I
INTRODUCTION
Once thought to be doomed by the computer, it is
becoming evident that shock tubes are just as important as
ever before. Someday computers will do almost all
aerodynamic predicting. However/ computational fluid
mechanics involves solution of some version of the
NavierStokes equation (in our case the Euler equation).
This still requires input of some physical parameters
which must be obtained from experiment. The passage of a
high strength shock wave over an object produces loading
on the object. Surface pressures and loads caused by the
blast wave are highly transient and generate severe
pressure gradients. Successful predictions of these loads
is necessary if the designer of any blast resistant
structure is to account for them. Such problems are
encountered in the design of missle bases, space launch
facilities, nuclear generating plants, and explosive or
propellent manufacturing facilities. In the past, the
study of shockinduced flow over bodies has been primarily
done by experiment. The experiments involve large shock
2
tubes or blast tubes and are extremely expensive to
perform. To reduce test costs a computer program has been
developed to simulate the shockinduced flow. The analysis
can be employed to determine the explosive configuration,
the optimum gas pressures, densities, and molecular
weights, and the shock tube configuration to achieve a
desired test condition. Although this idea of calculating
blast waves of large amplitude in a shock tube is not new;
variable area ducts and nonideal wave predictions have
not previously been considered in any detail. This study
then develops computer codes to model these two new
applications for shock tube testing as well as general
shock tube applications.
The first application deals with the expansion in
a shock tube. In order to accomadate a large scale model a
shock tube located at Sandia Corporation in Albequerque
was enlarged from a 6 foot diameter to a 19 foot diameter.
The driver section, containing the shock producing
explosives, has remained at a 6 foot diameter to reduce
costs and explosive handling problems, while the test
section has a 19 foot diameter. This configuration then
introduces a 40 foot expansion section from the 6 foot
driver section to the 19 foot test section. The second new
application is in trying to model the precursor phenomena.
A precursor or non ideal blast wave is formed when a
3
!
nuclear weapon is exploded at some height over a thermally
absorbing surface. The thermal radiation forms a hotair
layer above the earth surface through which the blast
front propagates at a higher velocity. The precersor wave
or nonideal wave apparatus is a straight tube but it is
unique in that various gases are placed in the tube in
order to develop a simulated nonideal wave. Such
facilities have been constructed at' SRI International and
at Kirkland Air Force Base. Test data for both of these
tests showed large descrepencies with the original
computer predictions. The computer code required extensive
modification to model, these new configurations and yield
accurate predictions for future test data.
CHAPTER II
PHYSICAL MODEL
Before we examine our method of solution let us
develop the model for the physical situation of a one
dimensional shock tube. Our blast tube facility can be
idealized by the classical pressure driven shock tube
shown in figure 2.1.a below.
Driver section Driven section
Distance
(b)
Figure 2.1 Initial conditions in a pressure driven shock tube
A diaphram seperates the region of highpressure "driver"
gas on the left, (region 4) from the region of lowpressure
5
gas on the right (region 1). The pressure distribution is
shown in figure 2.1.b. The gases in regions 1 and 4 will
generally have different temperatures and molecular
weights. The pressure discontinuity that exists after
rupture of the diaphram propogates into the test gas as a
strong shock and into the driver gas as a centered
rarefaction, in other words, when the diaphram in figure
2.1.a is broken a shock wave propogates into section 1 and
an expansion wave propogates into section 4. This picture
is sketched in figure 2.2.
Expansion wave
propogating to
the left
*
Contact surface
moving at velocity
of the gas behind
the shock wave
.
Normal shock wave
propogating to the
right with velocity W
Pressure
Figure 2.2 Flow in a shock tube after the diaphram is broken
As the normal shock wave moves to the right with velocity
6
W, it compresses the gas considerably thus resulting in
high pressures and temperatures of the gas behind it
(region 2) and induces a mass motion with velocity Up. The
contact surface is defined as the interface between the
driver and the driven gasses, and moves with velocity Un.
r*'
Note that the pressure and velocity are preserved across
the contact surface: ^3=^2 anc* u3=u2=up' i3ut t^ie entropy
changes discontinuously. From our figure it is seen that
the expansion wave propogates to the left, smoothly and
continuously decreasing the pressure in region 4 to the
lower value P3 behind the expansion wave. The flowfield in
the tube after the diaphram is broken (figure 2.2) is
completely determined by the conditions in regions 1 and 4
before the diaphram is broken. For our blast tube there is
no diaphram and the driver section contains explosives
whose sudden release of energy creates the high pressure
of region 4 in figure 2.1.b. Besides this difference the
mechanics of the flow field produced in the blast tube is
similar to that produced in the pressure driven shock
tube.
Returning now to figure 2.1 we define the ratio
P4/Pj_ as the diaphram pressure ratio. Along with the
initial conditions of the driver and driven, gas, P4/P1
determines the strength of the shock wave developed when
the diaphram is removed (or as in our case a large
7
pressure gradient is produced). As we know U3=U2=Up and
^2=I?3 across the contact surface. And it can be shown from
the RankineHugoniot equations that the particle velocity
induced behind the shock can be given as
a i
Up = u2 =
7i
Ei
Pi
2li
7l+l
+
7l ~ 1
7i + l
1/2
Also, since we are dealing with a simple centered
expansion wave we know from the Riemann invarients
u +
2 a
71
constant through the expansion wave.
We can evaluate this constant by applying the equation to
region 4
2a4 2a4
u4 +r = 0 H = constant
71 71
from which we obtain
a 7 1 / u
d4 2 \ d4
because
a = VtRT
8
and because the flow is isentropic
P_
Pa
1
7 1 / u
2 \a4
7/71
then
P_
Pa
1
71 / u_\
2 W
27/71
If we now apply this equation between the head and tail of
the expansion wave we have
P
Pa
741 / 3
2 \a4
274/741
Solving for U3 gives
113
2&A
741
741/274
Pa)
Since U2=U3 we can equate these equations giving
Qi {P2
7i VPi
1
271
7i + l
Ex +
Pi +
7irl
7i + l
1/2
2a4
74 1
(P3 \ 741/274
Which can be rearranged to yield
9
Pa
Pi
(74l)(ai/fl4)(i>2/A 1)
274/741
v'27i[27i + (7J+l)(ft/Pll)lJ
This is the shock tube1 equation which is used as a
fundamental part of the Godunov method, as will be shown.
This equation now gives the incident shock strength ^2^1
as an implicit function of the diaphram pressure ratio
P4/P1. Because a = y/'fRT = y/7(R/M) T the speed of
sound in a light gas is faster than in a heavy gas. Since
P2/P1 will be stronger if a^/a^ is made smaller, to
maximize the incident shock strength for a given P4/P1 the
driver gas should be a low molecular weight gas at high
temperature/ and the driven gas should be a high molecular
weight gas at low temperature. We are now ready to examine
the method of solution.
Computer Code
The computer program used in this study is called BTUBE,
which utilizes subroutines CELLPR, FIGAM, and AIR. BTUBE
is written in Fortran IV and uses the finite difference
technique called the Godunov Method (chapter III) The
basis for the Godunov method is to integrate the Euler
Equations for each cell using Green's theorm.This will
10
relate the volume changes in a cell of mass, momentum, and
energy, to the flux of these quantities across the cell
boundary. The shock tube equations are then used to
evaluate the fluxes and determine the properties in each
cell. This process is performed for each cell for each
time step. Further insight is provided if we examine the
path of the logic in the computer program BTUBE. Starting
in the main program the dimensions and initial conditions
of the cells are set up. The program has the option of
initializing cells containing different gases with
varying properties. Another option is for the tube
boundary conditions of either an opened or closed end. If
the end is opened the gas is allowed to escape. For the
closed end any incident wave is reflected back into the
tube. Driver pressures are initialized by the weight of
explosive in each cell. A time step is then calculated
based on the length of the cell and the sonic speed of the
gas. BTUBE then makes a call to subroutine CELLPR. This
subroutine first determines which condition has occured
in the cell; it determines if Wj_ and W2 are shock waves or
expansion waves. Given the pressure, density, velocity,
and values of gama across a cell boundary CELLPR then
finds the mass, momentum, and energy fluxes across the
cell boundary by solving the associated centered wave
problem using the Godunov method. These boundary
11
conditions are then used to solve the continuity,
momentum, and energy equations for each cell by a finite
difference scheme. Starting with the mass, momentum, and
energy we have
dp 9
Mass Tt + 6i(pu) = 0
d d , 2.
Momentum \PU) "I" \P "i Pu / = 0
9 9
Energy ^ (pE) + (pE + p)) = 0
E = e + u2 / 2
Equations of State
e E/ p(j 1)
Now let At=t2ti and make At sufficiently small that
the quantities p,u,P,E over S, the height of the cell, may
be considered as constant. Consider the integral of these
equations over At and S. Recall that Greens theorm for a
plane gives
dM
dY
dx dy =
+ M dx
Where M and N are continuous singled valued functions over
S bounded by the curve C and C is positive in the counter
clockwise direction. Starting with equation 2.2.1 we have
12
JLl%didzdyC*JIs{P)"~{l,)'
{pS)t2 (pS)h
and taking S as constant over
JIs It ^ dt dxdy^ JJ (pu) dx dy
Then applying Greens theorm to obtain
IL
A< (pu) dx dy = At
s &x
pu dy
Figure 2.3 shows the boundary C being made up of line
segments 123456781.
6
5
7
4
Y1 S Y2
3
l
2
Figure 2.3
13
The boundary conditions are:
Then
along 12 and 56
along 23, 45, 67, and 81
along 34, PU= ( p U)
along 78 P U= ( P U)
dy=0.0
U=0.0
and dy is positive
and dy is negative
/ pu dy = (pu)t2 y2 (pu)ti Vl
v C
And finally
dt dx dy % (pS)t2 (pS)ti + At ((pu)2 y2 {pu\ Vl) = 0
or in another form
(pS)ta = (ps)tl + ((pu)i Vn {pu)2 y2)
Equations 2.2.2 and 2.2.3 can be integrated in similar
fashion to obtain
(.puS)t2 = {puS)ti + At ((p + pu2) 1 y\ (p + pu2)2y2)
(,pES)t2 = {pES)ti + At {u (pE + p)j 2/1 u (pE + p)2 y2)
14
These are the finite difference equations which are then
used by the main program. The conditions in each cell are
now calculated for that time. The time is then increased
by. the time step and the whole process is started over
again. In summary CELLPR determines which shock condition
has occured, solves the shock tube equation, and returns
the fluxes to the main program.
Subroutine FIGAM iterates with subroutine AIR to
obtain the ratio of specific heats (gama) for the driver
cells. Subroutine AIR provides a "semiphysical fit" to
the equation of state which is generated from tabulated
values. The explosive in the driver cells is modeled by
inputing the energy of combustion and mass of the
explosive into FIGAM which which then iterates with AIR to
find the thermodynamic properties of the hot air. Thus,
the amount of energy which is put into the gas by the
explosive, is used to determine gama which is used for the
equation of state. We can then use these values to
determine the other properties in the driver cells.
Subroutines FIGAM and AIR are used only for the driver
cells, because the driver gets so hot. Gama is considered
constant for the "driven" cells of our model. A complete
computer listing of the program and its subroutines is
given in appendix A.
15
Expanding nozzle shock tube
Now that we have described the general model of a
shock tube we are ready to look at two applications which
modify this model. First consider the Sandia Labs
Thunderpipe facility in New Mexico which is illustrated
in figure 2.4. Both ends of the tube are open to the
atmosphere. The tests are powered by strands of Iremite
explosive placed inside the 6 ft. diameter "driver"
section which is approximatly 156 ft. long. As a result
of the energy release associated with the ignition of the
explosive, a blast wave is created. The blast wave first
travels thru the 40 ft. long expansion nozzle which
increases the diameter from 6 ft. to 19 ft. The shock
then moves toward the test section located 156 ft.
downstream. Although the tube is 19 ft. in diameter after
expansion it contains a flat concrete roadbed which
decreases the cross sectional area. An equivalent
diameter, due to this road bed, was found to be 17.35 ft..
In terms of the physical model the main deviation from the
simple shock tube is the smooth expansion nozzle, which is
handled in a straight forward manner. The computer
modeling of this geometry was accomplished by increasing
the radius in steps. Through the expansion section each
cells height is incrementally increased, eventually
I1
CTl
Figure 2.4 Thunder tube facility
17
stopping at the height of the driven section radius. So
the computer code was rewritten to handle two situations,
the straight sections and the expanding sections. The
uniqueness of this analysis was in the expanding section.
This approach explores the application of a one
dimensional method to a two dimensional problem. From the
finite difference equation for the continuity equation we
have,
MFS (J) = MFS(J) + DT (MFR(J 1) MFR(J)) YI
This is for the straight section of tube since the
difference term is multiplied by a constant height of YI
to obtain the mass flow. For the expanding section the
mass flux of the previous cell (MFR(Jl)) is multiplied by
the height YI, and the mass flux of the present cell
(MFR(J)) is multiplied by YI+DELY. This scheme is pictured
in figure 2.5.
DELY/2
YI+DELY
T DELY/2
_________i
Figure 2.5 Expanding cell geometry
18
The momentum and energy equations were handled in exactly
the same manner. Although this approach seems rather
simple, the results agreed well with the experimental data
as shown later.
Precursor shock tube
The second application we shall consider is the
precursor phenomena. A precursor is the term used for an
auxiliary blast wave caused by a nuclear explosion, which
propogates ahead of the main shock wave due to the
absorption of heat energy by a nonideal surface. If the
intense thermal radiation from a nuclear detonation
impinges on a thermally nonideal (heat absorbing) surface,
a hot thermal layer is formed near the surface. This layer
often includes smoke, dust, and other particulate matter,
and is formed ahead of the blast wave since the thermal
radiation travels faster than the airblast. The enhanced
sound speed in the Thermal layer near the ground surface
causes that portion of the incident blast wave which
enters the Thermal layer to be accelerated. Thus, the
shock wave propogating within the Thermal layer rapidly
outruns the incident blast wave resulting in the formation
of an aboveambient pressure region in the Thermal layer
that still has unshocked, ambient pressure air lying above
it. The thermal layer and it's interaction with the blast
19
Early Development Late Development
Figure 2.6 Precursor wave development
wave is shown in figure 2.6. Within the precursor region
considerable modification of the usual blast wave
characteristics may occur. From figure 2.6 we see that
after the precursor forms the main shock front usually no
longer extends to the ground. Between the bottom edge of
the main shock wave and the ground is a gap, through which
the energy that feeds the precursor may flow. Ahead of the
main shock front, the blast energy in the precursor is
free to follow the rapidly moving shock front in the
thermal layer and to propagate upward into the undisturbed
air ahead of the main shock front. This diverging flow
pattern within the precursor tends to weaken it, while the
20
energy which is continually fed into the precursor from
the main blast wave tends to strengthen the precursor
shock front. Interest in the precursor wave comes about
because the development of a precursor usually results in
lower peak overpressures, increased rise times to peak
overpressure, and increased dynamic pressure. The
increased dynamic pressure can be attributed to higher
particle velocities and/or greater air densities due to
dust and other debris picked up by the blast wave as it
passes over the ground. Another interesting phenomena
which is created by the precursor (but not considered in
this analysis) is the creation of vorticity. If we write
the voticity equation for inviscid, compressible flow we
have
Since the pressure gradient associated with the blast wave
is nearly orthogonal to the density gradient in the
Thermal layer underlying the ambient air, in general
where
fi = VxF
VP x Vp/ 0
21
Thus, the formation of a vortex centered along the top of
the Thermal layer with a strength which is proportional to
both the shock strength and the Thermal layer density
gradient is possible. This process is shown in figure 2.7.
below, and is worthy of further study.
Figure 2.7 Detailed precursor wave
Two test models have been developed to simulate the
dynamic and static pressure profiles of the precursor
wave. The first employs a bag of helium which is placed on
the floor of the shock tube to simulate the thermal layer,
as shown in figure 2.8. The helium then acts as the hot
thermal layer. This approach can accomodate a 1/30 scale
test, in a 20 ft. diameter shock tube. But since the
results from response tests are more accurate, and
instrumentation is easier on a larger model, it is
22
Shock wave 
Helium
Figure 2.8 Precursor shock tube
desirable to increase the scale of our test. The second
approach, which is used in this study, permits a 1/6
scale configuration in a 20 ft. diameter shock tube. This
model loads the cells in the shock tube with varying gases
having different molecular weights, starting with the
heaviest and decreasing down the tube. This formation
produces the desired one dimensional simulation of the
precursor wave. A satisfactory configuration was found to
be as shown on the next page in figure 2.9. The first
model develops the actual configuration of the precursor
wave but requires an enormous diameter to accomodate
larger scale testing. The second model does not physically
develop the precursor wave, but it does simulate the
dynamic and static pressure profiles experienced by the
precursor. We now turn our attention to the finite
difference method which was employed.
23
Distributed
explosive
density
SF
6
CO
2
Driver cells
Driven cells
Figure 2.9 Precursor shock tube
Test
section
He
CHAPTER III
GODUNOV METHOD
The problem considered here is the modification of
a onedimensional finite difference program in order to
simulate two applications of a shock tube. The finite
difference method used is an Eulerian differencing scheme
called the Godunov method. Godunovs method uses the
differential equations inconservative form; as a result,
the integration of these equations over the volume of a
cell results in equating the changes in mass, momentum,
and energy in the cell to the fluxes across the boundary
of these quantities during the time step. Pressure will be
considered as part of the momentum flux. At the start of a
time step the flow properties are taken to be uniform in
the cells on either side of a cell wall. As a result, the
fluxes are given by the solution of the centered wave
problem. Basically each cell can be visualized as a mini
shock tube, and our problem narrows down to solving the
shock tube equation for each cell. Our centered wave then
has an imaginary diaphram at the cell boundary with
conditions in the cell on that side of the boundary, and
25
conditions on the other side of the diaphram equal to the
conditions in the cell on the other side of the boundary,
as shown in figure 3.1. At the instant that the diaphram
is removed, two centered waves, W^ and W2, and a contact
surface leave the origin as depicted in figure 3.2. The
waves may be both shock waves, both expansion waves, or
one a shock wave and one an expansion wave. The Eulerian
boundary of each cell will lie in one of the four regions
shown in figure 3.2 for the entire timestep. Also, the
values of the flow variables in these regions are taken as
the values of these flow variables on the boundaries of
the cell. The main idea of the Godunov method is that of
obtaining the fluxes on the boundary of a cell by solving
the centered wave problem for that cell (as shown in
figure 3.2) Flow properties in the adjacent cells
(i1/2, j+1/2) and (i+1/2, j + 1/2) are known, and the flow
properties on the boundary of the cell correspond to the
properties on the vertical line X=0.0 Now let us examine
the equations of motion for our cell.
We shall start with wave Wj_. Assuming it is a
shock wave and applying conservation of mass and momentum
we obtain
3.1
{PV)I = (PV)n
Figure 3.1 One dimensional model for computation
Contact
Figure 3.2 Diagram of solution
27
(;p + pv2)j= (p + pV2)n 3.2
Note that the subscripts denote regions in figure 3.2 V
is the velocity of the fluid relative to the wave, P is
the pressure, and P is the density. For simplicity we set
the mass flux crossing wave Wj_ to bs so that equations 3.1
and 3.2 become
{pV)j = (PV)H = b3
Pi Pii + ba (V> Vn) = 0
If is the velocity of wave Wj_ and qn is' the gas
velocity normal to the cell boundary in figure 3.2 we
would have
Vi Vi i = (?n.j VWIJ)
= Qnjl 3.5
Substitution of equation 3.5 into equation 3.4 yields
Pi Pii + ba (qni ~ qmi) = 0 3 6
Through manipulation bs can be expressed in terms of Pjj
3.3
3.4
28
and known values in region I
ba =
PiPii
l/pil 1/pi.
,1/2
3.7
we can now write
J_ 1_ _1_
Pii Pi Pi
Pi
L Pii
1
and using the relationship
J__ i = i r(7+im+(7i)^/i
Pii Pi pi 1(7 ~ 1) Pi + (7 1) Pn
 1
we have
1_
Pi
2(PiPii)
(7 1)Pr + (7 + 1) Pii m
Substituting this back into equation 3,7 we have
bs =
Vei
3.7a
Equations 3.6 and 3.7a are the final equations when
considering wave W^ as a shock wave.
Now let us suppose that wave is an expansion
wave; then the equations of unsteady isentropic flow will
apply across the wave It would be convenient to have
29
0
these equations in the same form as equations 3.6 and 3.7a
because this allows us to express the equality of pressure
and velocity across the contact surface of our shocktube
model. Since is now a simple expansion wave propogating
upstream we know that the Riemann invariant J+ is constant
across W]_. So we have
?7l/ +
2 Cj
71
^ 2 Cn
Hr
71
3.8
Where C is the speed of sound
substitute
C2 =
T (Â£)
v^7 Rf
we have
If we
but
?nj Qnj!
V7
71
3.9
holds for an isentropic expansion. Substituting this into
equation 3.9 yields
_ 2v^
9nj 9njj 1
71
3.10
30
We can now transform this equaiton into the same form as
equation 3.5 by multiplying and dividing the righthand
side of the equation by pIIpI and rearranging to give
9.nji
71
~ (Pu Pi)
(PiPufl2(PilPii)U2''Pi
Pi, P,
3.11
Now if we define
be
7 1 / 7Pi 1 Pii/Pi
27 V 1if>i L1 (Pn/Pir1/2y
and substituting this relation into equation 3.11 gives
Pi Pil = be{qni qnn) 3.7b
Which is the same form as equaiton 3.5. Thus equation 3.5
holds across wave W^ with b=bs (equation 3.7a) for a shock
wave and with b=be (equation 3.7b) for an expansion wave.
For wave W2, from figure 3.2, we obtain a similar
set if equations only with a change of sign. We can
explain this change in sign in the following way. If wave
W2 is a shock wave our equation becomes
Piv Phi + a3 (Vjy Vi 11) = 0
3.12
31
Where as is used for the mass flux crossing wave W2. If
Vw2 is the velocity of wave W2,
Viv VlII = {Vy>2 ~~ Qniv) ~ ("^2 QnIll) = Qnm ~ 9n/v 3.13
Substitution of this value into equation 3.12 gives
Piv Pill + a(qnUJ qniv ) = 0
3.14
Note that the sign has changed from equation 3.5 because
wave W2 propogates downstream where wave W1 was propogat
ing upstream.
As seen from equation 3.14 as is analogous to bs
which is given from equation 3.7a as
CLg
(7 + 1) Pm + (7 ~ 1) Piv
,1/2
2/piv
3.15a
Now if wave W2 is as expansion wave, the sign change
occurs because the wave is a simple expansion wave
propogating downstream, and the J_ Riemann invariant is
constant across W2. Therefore the equation corresponding
to equation 3.8 is
2 Civ
^IV ,
71
2(7///
qniu
71
3.16
32
Also ae can be derived in the same way as be was to obtain
7 1 / 7Piv
ae~ 27 \jlfpjv
Similar to equation 3.14 always holds across wave W2
and equation 3.15a holds if W2 is a shock wave while
equation 3.15b holds if W2 is an expansion wave. If we now
let Pij = Pm = Pcs and Qnnj=
equations 3.5 and 3.14 we obtain
n bPiv + aPj + ab(qni qniv) 3.17a
cj 1
CL f 0
_ pl piv + aqnjv + bqni
UCs r
a + 3.17b
1 Phi/Piv
1 (PmlPjv)7'i/Ji
3.15b
with a and b being the appropriate choice of as or ae and
bs or be. These equations along with the applicable values
for a and b are the basic equations for the exact method
of Godunov. This set of equations requires an iterative
solution. For example; assume a value for Pcs; if pCs>pI
wave Wj_ is a shock wave and equation 3.7a gives the value
of b. If Pcs
3.7b is used to determine b. Similarly, if Pcs^iv wave w2
is a shock wave and equation and equation 3.15a gives the
33
value of a. And if Pcs
and equaiton 3.15b gives the value of a. Now we recompute
Pcs from equation 3.17a and continue this iteration until
the change in Pcs between iterations is very small. Ucs is
then given by equation 3.17b.
quantities needed and these come from either the shock
relations or the isentropic relation. So for Wj_ being a
shock wave we can. use the relation
and if Wj_ is an expansion wave we can use the isentropic
relation
The densities pu and Phi are now the final
_ (7 + 1)P + (7 ~ 1 )Pi
PU (7l)P + (7 + l)P/'
3.18a
3.18b
Similarly if W2 is a shock wave we have
(.7 + 1) Pcs + (7 1) Piv
p,n = (yi)P<. + h + i)Pivp,v
3.19a
and for W2 being as expansion wave
3.19b
CHAPTER IV
RESULTS
The success of any computer program is determined
by how well it matches with the experimental data. Figures
4.1 and 4.2 show the experimental data obtained from the
test section of the expanding shock tube test and the
nonideal shock wave test respectively. The computer
results for these two applications are shown in figures
4.3 and 4.4. Comparison of the computer predictions to the
test data show good agreement for both cases. Note that
for both cases the computer predicts an earlier arrival
time than is found in the test. A timing scheme for the
explosive was developed to try and simulate the burning of
the explosive and possible delay the arrival time. The
driver section was "turned on" cell by cell starting at
the last driver cell and marching backwards, at the burn
velocity of the explosive, until the end of the tube is
reached. This method did delay the arrival of the shock
some but not enough to compensate for the whole time gap.
We believe this time descrepency could possible be
explained by test procedures such as when the time starts
PRESSURE IN PSIG
72. OO
CO
(_n
Figure 4.1 Experimental data from the Thunder tube
CLfi(UJ(/)OODQ;ui I
Figure 4.2 Experimental d^ta from the nonideal
shock tube
PRESSURE (PSI)
TIME (SEC)
Figure 4.3 Computer generated data for the Thunder tube
30.0
25.0
20.0
5.0
10.0 H
30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 110.0 120.0 130.0
Time (msec)
Figure 4.4 Computer generated data from the nonideal
shock tube
39
with regards to when the explosive is ignited. Figures 4.5
and 4.6 show the computer plots superimposed upon the test
data where the computer data has been shifted in time to
account for the time gap. These figures show how closely
the computer data matches the test data.
So the computer programs utilizing the Godunov
finite difference method has been shown to be very
versatile in modeling the two applications studied in this
analysis. Further work is necessary to obtain even more
accurate predictions and to explain the timing
descrepency. However, the present model can be used with
confidence in predicting the pressure loads experienced in
testing with the expanding shock tube and the nonideal
shock tube.
PRESSURE IN PSIG
72. OO
OU OO
50. OO
UQ. OO
iio. oo
32 OO
2U OO
18. OO
0.00
0 oo
e. oo
o
Figure 4.5 Figure 4.3 imposed upon figure 4.1
_ m jo
TIME MSEC
Figure 4.6 Figure 4.4 imposed upon figure 4.2
REFERENCES
Anderson, J.D., Jr., Modern______Compressible flow. With
Historical perspective. McgrawHill, New York, 1982
Courant, R. and K.O. Friedrichs, Supersonic Flow and Shock
Waves, Interscience, New York, 1948
Gaydon, A. G. and I.R. Hurle, The_____Shock Tube in High
Temperature Chemical Physics, Reinhold, New York, 1963
Glasstones, S. and P. Dolan (ed.), The Effects of Nuclear
Weapons. U.S. Department of Defense, 1977
Liepmann, H.W. and A. Roshko, Elements of Gasdynamics.
Wiley, New York, 1957
McNamara, W. Flame_Computer. Code__for the Axisymmetric
Interaction of a Blast Wave with a Shock Layer on a Blunt
Body, Jour, of Spacecraft, Vol.4, No.6, June 1967
Shapiro, A.H., The_____Dynamics and Thermodynamics of
Compressible Fluid Flow, Vol. I. Ronald, New York, 1953
Zucrow, M.J. and J.D. Hoffman, Gas Dynamics. Whiley, New
York, 1977
APPENDIX A.
oo oouuuuuouoooouo o oo oo uoooouooo oo oo
44
PROGRAM BTUBE(INPUT,OUTPUT,STAT,DYN,STIMP.DYNIMP,
+ TAPE5=INPUT,TAPE6=OUTPUT,TAPE7=STAT,TAPE8=DYN,
+ TAPE9=STIMP,TAPE10=DYNIMP)
PROGRAM EMPLOYS GODUNOV FINITE DIFFERENCE
METHOD (SEE W. MCNAMARAJ.SPACECRAFT AND
ROCKETSJUNE 1967PP.790795) TO COMPUTE
THE PRESSURES INSIDE A DUCT.
********************
at*******************
THIS VERSION EMPLOYS
VARIABLE GAMMA USING
DOANNICKLE SUBROUTINE
********************
it*******************
REAL MFR(600).MOR(600),MFS(600),MOS(600),MFRS,MORS,MTEMP,
+ IP(600),IQ(600).IIP,IQ1
DIMENSION EFR(600).P(600),PI(600),RHO(600),U(600),EFS(600).
+ S(600),GAM(600),GAMZ(600),Q(600),UF(600),W(200),
+ P1I(600),Q1I(600),11(10)
COMMON/UNIBK/PA1,RH1,AS1.PAH
COMMON/UNBBK/JAHK,MFRS.MORS,EFRS,GAMC
COMMON/AIR1/LGAM,TEMP
DATA GC.TR.RG,RGO,ZPI/32.2,459.67,1717.6027,
+ 49753.7978.3.141592654/
READ(7,700) DX,Y,S1,S2,H
READ(7,710) JD.JA1,JA2.JA3,JN
READ(7,720) GAMA1.GAMA2,XMA1,XMA2,GAMT,EPE
READ(7,730) PAl.TA
READ(7.740) TLIM,TAD,KOUT.JNCK1,JNCK2
READ(7,750,END=5) W
NAMELIST /BANG/DX,Y,SI.S2,H, JD,JA1.JA2,JA3.JN,GAMA1,GAMA2,
+ XMA1.XMA2,GAMT,EPE.PAl.TA,TLIM,TAD,KOUT,JNCK1,JNCK2.W.II
READ(5,BANG)
WRITE(6,BANG)
5 PAHPA1*144.
LGAM0
JTERM=JN1
JLP=(JN2)/10
TARTA+TR
AS1=S($RT(GAMT*TAR*RG)
RH1=PAH/(RG*TAR)
EPA=PAH/((GAMT1.)*RH1)
RGA1=RG0/XMA1
45
c
c
c
c
c
c
RHA1PAH/(RGA1TAR)
RGA2RG0/XMA2
RHA2PAH/(RGA2*TAR)
UF0AS1
S(1)S1
S(JN)S2
P(1)PAH
P(JN)PAH
RHO(l)RHl
RHO(JN)RH1
PKD0.0
PI ( JN ) =0.0
U(1)=0.0
U(JN)0.0
Q(JN)0.0
GAM( 1)=GAMT
GAM( JN)GAMT
GAMZ(1)=GAMT
GAMZ ( JN ) <=GAMT
MFS(1)RH1*S(1)
MFS(JN)=RH1*S(JN)
M0S(l)0.0
MOS(JN)=0.0
EFSC1)=P(1) *S(1)/(GAMC1)l. )
EFS(JN)=P(JN)*S(JN)/(GAM(JN)1.)
MFR(1)0.0
MFR(JN)=0.0
M0R(1)=0.0
MOR(JN)0.0
EFR(l)0.0
EFR(JN)0.0
IP(l)0.0
IQ(l)0.0
P1KD0.0
IP(JN)=0.0
IQ( JN)=0.0
PlI(JN)0.0
Q1I(JN)=0.0
AEXPI*H*H/4.
WRITE(6,500)
WRITE(6,520) DX.Y.S1.S2.H.JD,JA1.JA2,JA3,JN.
+ GAMA1.GAMA2.XMA1,XMA2.GAMT.PA1.TA.EPE,TLIM.
+ TAD.KOUT,JNCK1.JNCK2.W(2)
WRITE(6,600)
DO 10 K=2,JTERM
10 S(K)DX*Y
DO 100 J=2,JTERM
IF(J.GT.JD) GOTO 60
RHOE=W(J)/(GC*AE*DX)
RHO(J)=RH0E+RH1
EAVG=(EPA*RH1+EPE*RH0E)/RH0(J)
ETEMP=EAVG 929.02267
RTEMPRHO(J)*0.51579185
GM=0.4
46
CALL AIR(ETEMP,RTEMP,GM)
GAMCJ)=l.+GM
P(J)GM*EAVG*RHO(J)
GAMZ(J)=GAM(J)
PI(J)(P(J)PAH)/144.
GOTO 90
60 IP(J.LT.JAl) GOTO 80
IF(J.GE.JA2) GOTO 75
PCJ)=PAH
PICJ)=0.0
RHO( J)=*RHA1
GAM(J)=GAMA1
GAMZ(J)=GAMA1
GOTO 90
75 IF(J.GE.JA3) GOTO 80
PCJ)=PAH
PICJ)=0.0
RHOCJ)=RHA2
GAMC J)=GAMA2
GAMZC J)=GAMA2
GOTO 90
80 PCJ)PAH
PICJ)=0.0
RH0CJ)RH1
GAMCJ)GAHT
GAMZ C J)=GAMT
90 MFSCJ)RHOCJ)*SCJ)
MOSCJ)=0.0
EFSCJ)=PCJ)*SCJ)/CGAMCJ)1.)
UCJ)=0.0
QCJ)0.0
MFRCJ)=0.0
MOR(J)=0.0
EFRCJ)=0.0
ipCJ)=o.o
IQCJ)=0.0
P1ICJ)=0.0
Q1ICJ)=0.0
UFCJ)=ABSCUCJ))+SQRTCPCJ)*GAMCJ)/RHOCJ))
UF0=AMAX1CUFcJ).UFO)
100 CONTINUE
c
DTDX/CTAD*UFO)
C
C
T=0.0
TPR=0.0
KOUNT=0
130 K0UNT=K0UNT+1
WRITEC7.1000) CTPR,CPICIICI)),I=1.10))
WRITEC8,1000) CTPR. CQCIICD) 1=1.10))
WRITEC9,1OO0) CTPR.(IPCllCD).1=1,10))
WRITEClO,1000) CTPR,CIQCIICI)).1=1,10))
IFCTPR.LT.0.0001) GOTO 132
IF(KOUNT.LT.KOUT) GOTO 150
132 WRITEC6,580) TPR
KOUNT=0
KL=1'
DO 145 Kl.JLP
KK=10*K
WRITEC6,590) CPICI).I=KL,KK)
134 WRITEC6.591) CUCD,I=KL,KK)
c
c
c
c
136 WRITEC6,592)
138 WRITEC6,593)
140 WRITEC6,594)
142 WRITEC6,595)
143 WRITEC6,596)
144 WRITE(6,600)
145 KL=KK+1
(RHO(I) I=KL,KK)
(GAM(I),I=KL,KK)
CQ(I),I=KL,KK)
(IP(I),I=KL,KK)
CH?(I),I=KL,KK)
WRITE(6,610)
150 DO 160 J=1,JTERM
IFCJNCK1.LT.1) GOTO 154
IF(J.NE.l) GOTO 154
P(J)=P(J+1)
RHO(J)RHO(J+l)
U(J)=U(J+1)
GAM(J)=GAM(J+1)
GOTO 156
154 IFCJNCK2.LT.1) GOTO 156
IFCJ.NE.JN1) GOTO 156
P(J+1)=P(J)
RHOCJ+l)=RHOCJ)
UCJ+DU(J)
GAM(J+1)=GAK(J)
156 CALL CELLPR(P(J),RHO(J),U( J),GAM(J),P(J+l),RHO(J+l),
+ U(J+1),GAM(J+1))
GAMZCJ)=GAMC
MFR(J)MFRS
MOR(J)=MORS
160 EFR(J)=EFRS
UF0=AS1
DO 300 J=2,JTERM
GATEMP=GAM(J)
GAMD=GAMCJ)*MFS(J)+DT*(GAMZCJ1)*MFR(J1)GAMZ(J)*
+ MFR(J))*Y
MFS(J)=MFS(J)+DT(MFR(J1)MFR(J))'Y
MTEMP=MFS(J)
STEMP=S(J)
IFCSTEMP.LT.0.0001) STEMP0.0001
RHOC J)=MTEMP/STEMP
MOS(J)=MOS(J)+DT*(MOR(Jl)MOR(J))*Y
IF(MTEMP.LT.l.E9) MTEMPl.E9
U(J)=MOSCJ)/MTEMP
EFSCJ)=EFSCJ')+DT*CEFRCJ1)EFRCJ))Y
PCJ) = C CEFSCJ)MOSC J)*UCJ)/2.)/STEMP)* CGAMCJ)l.)
0CJ)=RHOCJ)*U(J)*UCJ)/288.
UFCJ)=ABSCUCJ))+SQRTCPCJ)*GAMCJ)/RHOCJ))
GAMC J) =GAMD/MFS C J )
IFCGAMCJ).LE.1.) GAMC J)=GATEMP
UF0=AMAZ1CUFCJ),UFO)
PICJ)=CPCJ)PAH)/144.
IQ1 = C(?1ICJ')+Q( J) ) *DT/2.
I(?CJ)=I
Q1ICJ)=QCJ)
IlP=CPHCJ)+PlCJ))*DT/2.
IPCJ)=IPCJ)+HP
P1ICJ)=PICJ)
IFCJ.GE.JA1) GOTO 300
CALL FIGAMCPCJ),RHOCJ).GAMCJ))
300 CONTINUE
oo
48
C
DT=DX/CTAD*UFO)
C
TT+DT
IF(T.GT.TLIM) GOTO 400
TPR=T*1000.
GOTO 130
500 FORMATC 1H1.23X.24H""* PROGRAM BTUBE*"")
C
520 FORMATC10X,4HDX= ,F10.4,5X,3HY= ,F10.4,4X,4HS1= ,E12.2/
+ 10X,4HS2= ,E12.2,3X,3HH= ,F10.4,4X,4HJD= ,14/
+ 9X,5HJA1= ,14,9X,5HJA2= ,14,9X.5HJA3= ,14/
+ 10X.4HJN= ,14,7X,7HGAMA1= ,F8.4,3X.7HGAMA2= ,F8.4/
+ 8X,6HXMA1= ,F8.4,4X,6HXMA2= ,F8.4,4X,6HGAMT= ,F8.4/
+ 9X,5HPA1= .F8.4,6X,4HTA= .F10.4.3X,5HEPE= .E14.5/
+ 8X,6HTLIM= ,F8.4.5X,5HTAD= ,F8.4,4X,6HKOUT= ,14/
+ 7X,7HJNCK1= ,14,7X,7HJNCK2= ,I4,11X,3HW= ,F10.4//)
C
580 FORMATC/6X,3HT= ,F8.2,4HMSEC)
C
590 FORMATCIX,2HP=,IX,10(IX,F9.2))
C '
591 FORMATCIX,2HU=.IX,IOC IX,F9.2
C
592 FORMATCIX,2HR=,IX,10C1X.F9.6))
C
593 FORMATCIX,2HG=,IX,IOC IX,F9.6))
C
594 FORMATCIX,2HQ=,IX,10C1X.F9.2))
C
595 FORMATCIX,3HIP=,IOCIX,F9.2)) '
C
596 FORMATCIX,3HIQ=,10C1X.F9.2))
C
600 FORMATC//)
C
610 FORMATC////)
C
700 FORMATC2F12.4,2E12.2,F12.4)
C
710 FORMATC514)
C
720 FORMATC5F10.4.E14.5)
C
730 FORMATC2F10.4)
C
740 FORMATC2F10.4,314)
C
750 FORMATC7F10.4)
C
1000 FORMAT CllClX.ElO.4))
C
400 WRITEC6.600)
END
C
SUBROUTINE CELLPRC PA,RKOA,UA,GAA,PB,RHOB,UB,GAB)
C
REAL MA,MB,MFRS,MORS
C
COMMON/UNIBK/PA1,RH1,AS1,PAH
49
COMMON/UNBBK/JAHK,MFRS,MORS,EFRS,GAMC
C
FHL(PTE)=PHHPTE/((1.+SWT*BEQ*(UHUL)/AHRBEP*
+ (PTE1.)*SQRT(1./(BLH*PRR*(ALP*PTE+1.))))**(1./BEP))
C
C
PCK1ABS(PAPAH)
PCK2=ABS(PBPAH)
EA=PA/((GAA1.)*RHOA)
EB=PB/((GAB1.)*RHOB)
UAAABS(UA)
UAB=ABS(UB)
IF(PCK1.LT..01.AND.PCK2.LT..01.AND.UAA.LT..01.AND.
+ UAB.LT..01) GOTO 100
PAB=PA/PB
IF(PAB.GE.l.O) GOTO 10
PHPB
RHOHRHOB
OH=UB
GAMHGAB
SWT1.0
PL=PA
RHOLRHOA
ULUA
GAMLGAA 
PHHl./PAB
GOTO 20
10 PHPA
RHOHRHOA
UHUA
GAMH=GAA
SWT1.0
PLPB
RHOLRHOB
ULUB
GAMLGAB
PHHPAB
20 PRR=PH*RHOL/(PL*RHOH)
ALR=SQRT(GAML* PL/RHOL)
AHR=SQRT(GAMH* PH/RHOH)
ALP(GAML+1.)/(GAML1.)
BEP=(GAMH1.)/(2.*GAMH)
BEQ=GAMH*BEP
BLH=(GAML1.)/(2.*GAMH)
IF(PRR.LT.200.) GOTO 21
PF2PHH* (1. + SWT*BE
GOTO 26
21 IF(PRR.GT..001) GOTO 22
PF21.0001
GOTO 26
22 IP=0
PF1PHH
PF2PHH/2.
23 FP1=FHL(PF1)
FP2FHL(PF2)
IF(ABS(FP2).LE..02) GOTO 26
IF(ABS(FP2FPl).LT.l.E9) GOTO 26
IPIP+1
IF(IP.LT.20) GOTO 24
WRITE(6,800) FP1.FP2
800 FORMAT(10K,20HNO CONVERGENCEFP1 ,F12.3,2X,5HFP2= ,F12.3)
GOTO 26
50
24 PF3=PF2FP2*(PF2PF1)/(FP2FP1)
PF1=PF2
PF2PF3
GOTO 23
26 PSFPF2
IFCPSF.LT.l.E9) PSF=l.E9
PII=PSF*PL
IFCSWT.GE.0.0) GOTO 30
MA=SQRT(((GAML+1.)*PII+CGAML1.)*PA)RHOA/2.)
PI2=PII
IF( PII/PB. GT. 0.999) PI2=0.999*PB
MB=BEP*SQRT(GAMH*PB*RHOB)*C1.PI2/PB)/(1.(PI2/PB)
+ **BEP)
GOTO 45
30 PI2=PII
IFCPII/PA.GT.0.999) PI2=0.999*PA
MA=BEP SQRT(GAMH PA *RHOA)*(1.PI2/PA)/(1.(PI2/PA)
+ **BEP)
MB=SQRT(C(GAML+1.)*PII+(GAML1.)*PB)*RH0B/2.)
45 URA=UA
URB=UB
VA=UAMA/RHOA
VB=UB+MB/RHOB
PCS=CMA*PB+MB*PA+MA*MB*(URAURB))/(MA+MB)
IFCPCS.LT.l.E9) PCS=l.E9
UCS=CPAPB+MA*URA+MB*URB)/CMA+MB)
RH0A2RH0A* C CGAA+1.)PCS+CGAA1.)* PA)/C CGAA1.)*PCS+
+ CGAA+1.)*PA)
RHOB3RHOB *(CGAB+1.)* PCS+CGAB1.)* PB)/CCGAB1.)* PCS+
+ CGAB+1.)*PB)
IFCVB.GT.0.O) GOTO 60
PC=PB
RHOCRHOB
UC=UB
GAMC=GAB
JAHK.1
GOTO 90
60 IFCUCS.GT.O.O) GOTO 70
PC=PCS
RH0C=RH0B3
UC=DCS
IFCEB.LE.1.E+9.0R.GAB.GT.1.4) GOTO 65
GAMRGAB
CALL FIGAMCPC.RHOC.GAMR)
GAMC=GAMR
GOTO 67
65 GAMC=GAB
67 JAHK=1
GOTO 90
70 IFCVA.GE.0.0) GOTO 80
PC=PCS
RHOC=RHOA2
UC^UCS
IFCEA.LE. 1.E+9.0R.GAA.GT. 1.4) GOTO 75
GAMR=GAA
CALL FIGAMCPC.RHOC.GAMR)
GAMCGAMR
GOTO 77
75 GAMC=GAA
77 JAHK=2
GOTO 90
80 PC=PA
nooo noo o o o nnooooon ooo
RHOC=RHOA
UC=UA
GAMC=GAA
JAHK=2
90 MFRS=RHOC*UC
MORS=PC+MFRS*UC
EFRS=(GAMC*PC/(GAMC1.)+MFRS*UC/2.)*UC
RETURN
100 MFRS0.0
MORS=0.0
EFRS0.0
JAHK=0
GAMC=GAA
RETURN
END
SUBROUTINE FIGAMC PG,RHOG,GAMG)
THIS SUBROUTINE ASSUMES THAT PRESSURE CALL IS IN
UNITSOF LBF/FT2 AND DENSITY CALL IS IN UNITS OF
LBFSEC2/FT4. IT THEN ITERATES TO FIND THE VALUES
OF GAMMA AND TEMPERATURE APPROPRIATE TO THOSE
CONDITIONS.
PM=PG*478.8
RHOl=RH0G*0.51579185
RH0M=RH01
GM=GAMG1.0
GMOLD=GM
DO 10 1=1,40
GM=(GM+GMOLD)/2.
GMOLD=GM
EE=PM/(GM*RHOM)
IFCEE.LE.2.E+13) GOTO 5
EE=2.E+13
RHOM=PM/(GM*EE)
5 CALL AIR C EE,RHOM,GM)
IF(ABS((GMGMOLD)/GM).LT.l.E5) GOTO 20
RH0MRH01
IFCI.LT.40) GOTO 10
WRITEC6.600) GM.GMOLD
GOTO 20
10 CONTINUE
20 GAMG=GM+1.0
RETURN
600 FORMAT(1OX,22HNO CONVERGENCEGMOLD= .F9.6.2X,
+ 4HGM= ,F9.6/)
END
SUBROUTINE AIR(EEE,RRR,GMONE)
THIS SUBROUTINE ADAPTED FROM L.R.DOAN AND G.H.NICKEL
"A SUBROUTINE FOR THE EQUATION OF STATE OF AIR"
I
52
C AFWL RTD(WLR)TM632,MAY 1963WITH UPDATES BY
C C.NEEDHAM.ET.AL.IN THE DNA 1ST STANDARD
C
C EEE IS ENERGY(ERGS/GM)
C RRR IS DENSITY(GM/CM3)
C GMONE IS GAMMA1.0
C
COMMON/AIR1/LGAM,TEMP
C
C
E=EEE*1.E10
RHO=RRR/1.293E3
El=(8.5E)/.975
IF(ABS(E1).LT.5.0) GOTO 20
C
IF(El.GT.O.O) GOTO 10
C
F0=0.0
FON=EXP(E/6.63)
WS=0.0
GOTO 30
C
10 F0=EXP(E/4.46)
F0N=0.0
WS=1.0
GOTO 30
C
20 DEI0.975*RHO**0.05
EE1=8.5+0.155042*ALOG(RHO)
E1=(EE1E)/DE1
VS=1./(EXP(E1)+1.)
FOEXPCE/4.46)*WS
FON=EXP(E/6.63)*(1.WS)
C
30 BETA=0.0
IF(E.LT.1.0) GOTO 40
C
BETACO.0208*WS+0.0139*(1.WS))*ALOG(E)
E2=(E40.)/3.0
IF(ABS(E2).LT.5.) GOTO 60
C
IF(E2.GT.0.0) GOTO 50
C
40 FN=0.0
WS=0.0
GOTO 70
C
50 FN=EXP(E/25.5)
WS=1.0
GOTO 70
C
60 DE2=4.*RHO**0.085
EE2=45.*RHO**0.0157
E2=(EEE2)/DE2
WS=1./(EXP(E2)+1.)
FN=EXP(E/25.5)*WS
C
70 E3=(E160.)/6.0
BETABETA *(1.WS)+0.045 WS
FE=0.0
IFCE3.GT.C5.0)) FE=1./(EXP(E3)+1.)
GM0NE=(0.161+0.255*FO+0.28*FON+0.137*FN+0.05 *FE)*
53
+ RHO*BETA
C
IF(LGAM.LT.l) RETURN
C
RHOLN=ALOG(RHO)
RHOLSRHOLN RHOLN
GMGMONE
P=GM*EEE*RRR
IF(E.GT.120.) GOTO 100
F=0.9722.7143E3RHOLN+6.582549E5 *RHOLS
G*0.026456.418873E4*RH0LN2.338785E5 *RHOLS
H9.21E5+5.971549E6*RHOIiN+3.923123E7*RHOLS
CON1=3480.*GM*E
.CON2=F+G+H*E*E
TEMP=C0N1/C0N2
GOTO 200
C
100 ALR=AL0G10(RRR)
Cl=l.69081E5+2.99265E7*ALR
C2=0.769813+3.8618E3*ALR
TEMP=C1*EEE**02+102.6275735
C
200 BETT=(E3)/0.66
IFCBETT.GT.10.) RETURN
SWIT=1./(EXP(BETT)+1.)
TEMP=C0N1/(CON2 *(1.SWIT)+SWIT)
IF(TEMP.GT.O.O) RETURN
ALRALOGIO(RRR)
Cl1.69081E5+2.99265E7*ALR
C20.769813+3.8618E3*ALR
TEMPC1*EEE *C2+102.6275735
RETURN
END
APPENDIX B.
OOOOOOO OOOOOOOQOOOOOOQOOOO
55
PROGRAM TBTUBE
PROGRAM EMPLOYS GODUNOV FINITE DIFFERENCE
METHOD (SEE W. MCNAMARAJ.SPACECRAFT AND
ROCKETSJUNE 1967PP.790795) TO COMPUTE
THE PRESSURES INSIDE A DUCT.
********************
********************
THIS VERSION EMPLOYS
VARIABLE GAMMA USING
DOANNICKLE SUBROUTINE
********************
********************
REAL MFRC600).MOR(600),MFS(600),MOS(600),MFRS,MORS,MTEMP,
R IP(600).IQ(600),IIP,IQ1
DIMENSION EFR(600),P(600),PI(600),RHO(600),U(600),EFS(600),
D S(600),GAM(600),GAMZ(600),Q(600),UF(600),W(140).
D P1I(600),Q1I(600),11(10)
COMMON/UNIBK/PA1,RH1,AS1,PAH
COMMON/UNBBK/JAHK,MFRS,MORS,EFRS.GAMC
COMMON/AIR1/LGAM,TEMP
DATA GC,TR,RG,RGO,XPI/32.2,459.67,1717.6027.
D 49753.79741,3.141592654/
READ(7,700) DX.S1.S2.H
READ(7,710) JD,JA1,JA2,JA3,JN
READ(7,720) GAMA1,GAMA2,XMA1,XMA2,GAMT,EPE
READ(7,730) PAl.TA
READ(7,740) TLIM,TAD,KOUT,JNCK1,JNCK2
READ(7,750) Y1.YF,JEXPS.JEXPI
READ(7,760,END5) W
NAMELIST /BANG/DX,SI,S2,H,JN,JD,JA1,JA2,JA3,GAMA1,GAMA2,
N XMA1,XMA2,GAMT,EPE,PAl.TA,TLIM,TAD,KOUT,JNCK1.JNCK2,
N Yl.YF,JEXPS,JEXPI,W,II
READ(5,BANG)
WRITE(6,BANG)
5 PAHPA1*144.
LGAM0
JTERMJN1
JLP(JN2)/10
TAR*TA+TR
AS1S9RT(GAMT*TAR*RG)
RH1PAH/(RG*TAR)
EPAPAH/((GAMT1.)*RH1)
UFOAS1
s(i)si
S(JN)S2
P(1)PAH
P(JN)PAH
RHO(l)RHl
RHO(JN)RHl
PI(l)0.0
PI(JN)0.0
U(l)0.0
Q(l)0.0
U(JN)0.0
Q(JE)0.0
GAMCDGAMT
GAM(JN)GAMT
GAM2(1)GAMT
GAMZ(JN)GAMT
MFS(1)RH1*S(1)
MFS(JN)RH1*S(JN)
MOS(1)0.0
MOS(JN)0.0
EFS(1)P(1)*S(1)/(GAM(1)1.)
EFS(JN)P(JN)*S(JN)/(GAM(JN)1.)
MFR(l)0.0
MFR(tTN)0.0
M0R(l)0.0
MOR(JN)0.0
EFR(l)0.0
EFR(JN)0.0
IP(l)0.0
IQCD0.0
PiKDo.o
QlKD0.0
IP(JE)0.0
IQ(JH)0.0
P1I(JE)0.0
QlKJTO0.0
AEXPI*H*H/4.
*RITE(6,500)
WRITE(6,510) Y1
WRITE(6,511) YF
TOITE(6,512) JEEPS
WRITE(6,513) JEXPI
WRITE(6,520) DX,SI,S2,H,JD,JA1,JA2,JA3,JN,
+ GAMA1,GAMA2,XMA1,XMA2,GAMT,PA1,TA,EPE,TLIM,
+ TAD,KOUT,JNCK1,JUCK2,W(2)
WRITE(6,600)
JEXPEJEXPS+JEXPI
YIY1
JEXPS1JEXPS1
DO 10 K=2,JEXPS1
10 S(K)DX*YI
YYI
DELY(YFYI)/JEXPI
JEXPE1JEXPE1
DO 15 KJEXPS.JEXPEl
YY+DELY
S(K)DX*Y
15 CONTINDE
DO 20 KJEXPE,JTERM
20 S(K)DX*Y
DO 100 J2,JTERM
IF(J.GT.JD) GOTO 80
RH0EW(J) / (GC*AE*DX)
RH0(J)RH0E+RH1
EAVG(EPA*RH1+EPE*RH0E)/RH0(J)
ETEHPEAVG*929.02267
RTEMPRH0(J)*0.51579185
GM0.4
CALL AIR(ETEMP,RTEMP,GM)
GAM(J)1.+GM
P(J)GM*EAVG*RHO(J)
GAMZ(J)GAM(J)
PI(J)(P(J)PAH)/144.
GOTO 90
80 P(J)PAH
PI(J)0.0
RH0(J)RH1
GAM(J)GAMT
GAMZ(J)GAMT
90 MFS(J)RHO(J)*S(J)
MOS(J)0.0
EFS(J)P.(J)*S(J)/(GAM(J)1.)
D(J)0.0
Q(J)0.0
MFR(J)0.0
MOR(J)0.0
EFR(J)0.0
IP(J)0.0
IQ(J)0.0
PlI(J)0.0
QlI(J)0.0
DF(J)ABS(U(J))+SQRT(P(J)*GAM(J)/RHO(J))
DFOAMAXl(DF(J),DFO)
100 CONTINDE
DTDX/(TAD*DFO)
T0.0
TPR0.0
KODNTO
130 KODNTKOUNT+1
YIY1
WRITE(7,1000) (TPR,(PI(II(I)),11,10))
WRITE(8,1000) (TPR,(Q(II(I)),11,10))
WRITE(9,1000) (TPR,(IP(II(I)),11,10))
WRITE(10,1000) (TPR,(IQ(II(I)),11,10))
IF(TPR.LT.0.0001) GOTO 132
IF(KOUNT.LT.KOUT) GOTO 150
132 WRITE(6,580) TPR
KODNTO
KL1
DO 145 Kl.JLP
(PI(I),I=KL,KK)
(U(I),I=KL,KK)
(RHO(I),I=KL,EE)
CGAM(I),I=EL,EK)'
(Q(I),I=EL,EE)
(IP(I).I=EL,EE)
(10(1),I=EL,EE)
(S(I),I=EL,EE)
EK=10*E
WRITE(6,590)
134 WRITE(6,591)
136 WRITE(6,592)
138 WRITE(6,593)
140 WRITE(6,594)
142 WRITE(6,595)
143 WRITE(6,596)
WRITE(6,597)
144 WRITE(6,600)
145 ELEE+1
WRITE(6,610)
150 DO 160 Jl.JTERM
IF(JNCEl.LT.l) GOTO 154
IF(J.NE.l) GOTO 154
P(J)P(J+1)
RH0(J)RH0(J+1)
U(J)U(J+1)
GAM(J)GAM(J+1)
GOTO 156
154 IF(JNCE2.LT.1) GOTO 156
IF(J.NE.JNl) GOTO 156
P(J+1)P(J)
RH0(J+1)RH0(J)
D(J+1)D(J)
GAM(J+1)GAM(J)
156 CALL CELLPR(P(J),RHO(J),U(J),GAM(J),P(J+l),RH0(J+1),
+ U(J+1).GAM(J+1))
GAMZ(J)=GAMC
MFR(J)MFRS
MOR(J)MORS
160 EFR(J)EFRS
UF0AS1
JSTART=2
JENDJEEPS1 ,
210 DO 200 JJSTART,JEND
GATEMP=GAM(J)
GAMDGAM(J)*MFS(J)+DT*(GAMZ(J1)*MFR(J'1)GAMZ(J)*
+ MFR(J))*YI
MFS(J)MFS(J)+DT*(MFR(J1)MFR(J))*YI
MTEMP=MFS(J)
STEMPS(J)
IF(STEMP.LT.0.0001) STEMP=0.0001
RHO(J)MTEMP/STEMP
MOS(J)MOS(J)+DT*(MOR(Jl)MOR(J))*YI
IF(MTEMP.LT.l.E9) MTEMPl.E9 '
U(J)MOS(J)/MTEMP
EFS(J)EFS(J)+DT*(EFR(J1)EFR(J))*YI
P(J)((EFS(J)MOS(J)*U(J)/2.)/STEMP)*(GAM(J)l.)
Q(J)RHO(J)*U(J)*U(J)/288.
UF(J)ABS(U(J))+SQRT(P(J)*GAM(J)/RHO(J))
GAM(J)=GAMD/MFS(J)
IF(GAM(J).LE.l.) GAM(J)=GATEMP
UF0=AMAX1(UF(J),UFO)
PI(J)(P(J)PAH)/144.
IQl(QlI(J)+
IQ(J)=*IQ(J)+IQ1
QII(J)QCJ)
I1P(P1I(J)+PI(J))*DT/2.
IP(J)IP(J)+I1P
P1I(J)=PI(J)
IF(J.GT.JD) GO TO 200
CALL FIGAM(P(J),RHO(J),GAM(J))
200 CONTINUE
IFCJ.GE.JTERM) GO TO 220
DO 300 JJEXPS,JEXPE
GATEMP=GAM(J)
GAMDGAM(J)*MFSC
+ MFR(J)*(YI+DELY))
MFS(J)=MFS(J)+DT*(MFR(J1)*YIMFR(J)*(YI+DELY))
MTEMPMFSCJ)
STEMPS(J)
IF(STEMP.LT.0.0001) STEMP0.0001
RHO(J)MTEMP/STEMP
MOS(J)MOS(J)+DT*(MOR(Jl)*YIMOR(J)*(YI+DELY))
IF(MTEMP.LT.l.E9) MTEMPl.E9
U(J)MOS(J)/MTEMP
EFS(J)EFS(J)+DT*(EFR(Jl)*YIEFR(J)*(YI+DELY))
P(J)((EFS(J)MOS(J)*U(J)/2.)/STEMP)*(GAM(J)1.)
Q(J)RHO(J)*U(J)*U(J)/288.
UF(J)ABS(U(J))+S
GAM(J)GAMD/MFS(J)
IF(GAM(J).LE.1.) GAMCJ)GATEMP
UFOAMAZlCUF(J),UFO)
PI(J)(P(J)PAH)/144.
IQ1?1I( J)+
IQ(J)=IQ(J)+K31
Q1I(J)
I1P(P1I(J)+PI(J))*DT/2.
IP(
PII(J)PKJ)
YIYI+DELY
IF(J.GT.JD) GO TO 300
CALL FIGAM(P(J),RHO(J),GAM(J))
300 CONTINUE
JSTARTJEXPE+1
JENDJTERM
GO TO 210
220 CONTINUE
DTDX/(TAD*UFO)
T=T+DT
IF(T.GT.TLIM) GOTO 400
TPRT*1000.
GOTO 130
500 FORMATC1H1,23X,24H* * PROGRAM BTUBE*** *)
510 FORMATC10X.18HINITIAL DIAMETER .F8.2/)
511 FORMATC10X.16HFINAL DIAMETER ,F8.2/)
512 FORMATC10X,26HEXPANSION STARTS AT CELL ,14/)
513 FORMATC10X,30HNUMBER OF STEPS IN EXPANSION ,14/)
APPENDIX C.
61
c
c
C PROGRAM BTUBE
C
C
C PROGRAM EMPLOYS GODUNOV FINITE DIFFERENCE
C METHOD (SEE W. MCNAMARAJ.SPACECRAFT AND
C ROCKETSJUNE 1967PP.790795) TO COMPUTE
C THE PRESSURES INSIDE A DUCT.
C
c ********************
Q ********************
C THIS VERSION EMPLOYS
C VARIABLE GAMMA USING
C DOANNICKLE SUBROUTINE
Q ********************
Q ft*******************
REAL MFR(600),MOR(600),MFS(600),MOS(600),MFRS,MORS,MTEMP,
R IP(600),10(600),IIP,101
DIMENSION EFRC600),P(600),PI(600),RHO(600),U(600),EFS(600),
D S(600),GAM(600),GAMZ(600),0(600),UF(600),W(300),
D P1I(600).011(600),11(10)
COMMON/UNIBK/PA1,RH1,AS1,PAH
COMMON/UNBBK/JAHK,MFRS,MORS,EFRS,GAMC
COMMON/AIR1/LGAM,TEMP
DATA GC,TR,RG,RGO,XPI/32.2,459.67.1717.6027,
D 49753.7978,3.141592654/
NAMELIST /BANG/DX,Y,SI,S2,H,JD,JA1,JA2,JA3,JA4,JN,GAMA1,GAMA2,
N GAMA3,3MA1,XMA2,XMA3,GAMT.EPE,PA1,TA,TLIM,VEXP,TAD,
N KOUT,JNCK1,JNCK2,W,II
READ(5,BANG)
WRITE(6,BANG)
5 PAH=PA1* 144.
LGAM=0
JTERMJN1
JLP=(JN2)/10
TARTA+TR
AS1S0RT(GAMT*TAR*RG)
RH1=PAH/(RG*TAR)
EPAPAH/((GAMT1.)*RH1)
RGA1RGO/XMA1
RHA1PAH/(RGA1*TAR)
RGA2RG0/XMA2
RHA2=PAH/(RGA2*TAR)
RGA3RGO/XMA3
RHA3=PAH/(RGA3*TAR)
UF0=AS1
S(1)=S1
62
SCJN)S2
PCDPAH
PCJN)PAH
RHOCDRHl
RHO(JN)RHl
PICDO.O
PICJN)0.0
UCD0.0
QCD0.0
U(JN)0.0
Q(JN)=0.0
GAMCl)GAMT
GAMCJlOGAHT
GAMZ(1)=GAMT
GAKZ C JN )=GAMT
MFSCl)RHl*SCl)
MFSCJN)RH1*SCJN)
MOSCl)0.0
MOSC JlO0.0
EFSCl)PCl)*SCl)/CGAMCl)l.)
EFS(JN)P(JN)*SCJN)/(GAM(JK)1.)
MFRCl)=0.0
MFRCJN)0.0
MORCDO.O
MORCJN)0.0
EFR(l)0.0
EFRCJH)0.0
IP(1)0.0
IQCD0.0
P1ICD0.0
QiiCDo.o
IPCJN)=0.0
IQ(JH)0.0
P1ICJN)=0.0
Q1I(JN)0.0
AEXPI*H*H/4.
WRITEC6.500)
WRITEC6,520) DX,Y,S1,S2,H,JD,JA1,JA2,JA3.JN,
+ GAMA1,GAMA2,XMA1,XMA2,GAMT,PA1.TA,EPE,TLIM,
+ TAD.KODT,JNCK1,JNCK2,W(2)1
WRITE(6,600)
DO 10 K=2,JTERM
10 S(K)=DX*Y
DO 100 J=2,JTERM
60 IF(J.LT.JAl) GOTO 80
IF(J.GE.JA2) GOTO 75
PCJ)PAH
PI(J)0.0
RH0CJ)RHA1
GAM(J)GAMA1
GAMZ(J)=GAMA1
GOTO 90
75 IFCJ.GE.JA3) GOTO 77
PCJ)PAH
PI(J)=0.0
RH0(J)=RHA2
GAM( J)GAMA2
GAMZ(J)GAMA2
GOTO 90
77 IF(J.GE.JA4) GOTO 80
P(J)PAH
PI(J)0.0
RHO( J)=*RHA3
GAM(J)GAMA3
GAMZ(J)=GAMA3
GOTO 90
80 P(J)PAH
PI(J)0.0
RH0(J)=RH1
GAM(J)GAMT
GAMZ(J)=GAMT
90 MFS(J)RHO(J)*S(J)
MOS(J)0.0
EFS(J)P(J)*S(J)/(GAM(J)1.)
U(J)0.0
Q(J)0.0
MFR(J)0.0
MOR(J)0.0
EFR(J)0.0
IP(J)0.0
IQ(J)0.0
P1I(J)=0.0
Q1I(J)=0.0
UF(J)ABSCU(J))+SQRT(P(J)*GAM(J)/RHO(J))
UF0=AMAX1(UF(J),UFO)
100 CONTINUE
DT=DX/(TAD*UFO)
TPR0.0
JSJD
ZREM=0.0
K0UNT=0
130 KOUNTKOUNT+1
XST=C VEXP*DT)+XREM
IF(JS.LE.2)GO TO 20
30 IF(XST.GE.DX)THEN
IF(JS.LE.2)GO TO 20
RHOEWCJS)/(GC*AE*DX)
RHO(JS)RHOE+RH1
EAVG=(EPA RH1+EPE RHOE)/RHO(JS)
ETEMPEAVG*929.02267
RTEMPRHOCJS)*0.51579185
GM=0.4
CALL AIR(ETEMP,RTEMP,GM)
GAM(JS)=1.+GM
P(JS)=GM*EAVG*RHO(JS)
GAMZ(JS)=GAM(JS)
PI(JS)=(P(JS)PAH)/144.
MFS(JS)=RHO(JS)*S(JS)
EFS(JS)=P(JS)*S(JS)/(GAM(JS)1.)
UF(JS)=ABS(U(JS))+SQRT(P(JS)*GAM(JS)/RHO(JS))
JSJSl
XST=XSTDX
IF(XST.GT.0.0)THEN
XREM=XST
ELSE
XREM=00
END IF
ELSE
XREMDXXST
END IF
IF(XREM.GE.DX) GO TO 30
20 CONTINUE
WRITE(7,1000) (TPR,(PICII(I)),1=1,10))
WRITE(8,1000) (TPR,(Q(II(I)),1=1,10))
WRITE(9,1000) (TPR,(IP(II(I)),I=1,10))
WRITE(10,1000) (TPR,(IQ(II(I)),1=1,10))
IF(TPR.LT.0.0001) GOTO 132
IF(KOUNT.LT.KOUT) GOTO 150
132 WRITE(6,580) TPR
KOUNT=0
KL=1
DO 145 Kl.JLP
KK=10*K
WRITE(6,590)
134 WRITE(6,591)
136 WRITE(6,592)
138 WRITE(6,593)
140 WRITE(6,594)
142 WRITE(6,595)
143 WRITE(6,596)
144 WRITE(6,600)
145 KL=KK+1
(PI(I),I=KL,KK)
(U(I),I=KL,KK)
(RHO(I),I=KL,KX)
(GAM(I),I=KL,KK)
?(I),I=KL,KK)
(IP(I),I=KL,KX)
(I
WRITE(6,610)
150 DO 160 J=1,JTERM
IF(JNCXl.LT.l) GOTO 154
IF(J.NE.l) GOTO 154
P(J)=P(J+1)
RHO(J)=RHO(J+l)
U(J)=U(J+1)
GAM(J)=GAM(J+1)
GOTO 156
154 IF(JNCE2.LT.1) GOTO 156
IF(J.NE.JNl) GOTO 156
P(J+1)=P(J)
RHO(J+l)=RHO(J)
U(J+1)=U(J)
GAM(J+1)=GAM(J)
156 CALL CELLPR(P(J),RHO(J),U(J),GAM(J),P(J+l),RHO(J+l),
+ U(J+l),GAM(J+l))
GAMZ(J)=GAMC
MFR(J)=MFRS
MOR(J)=MORS
65
160 EFRCJ)EFRS
C
UFOAS1
C
DO 300 J2.JTERM
GATEMPGAMCJ)
GAMDGAMCJ)*MFSCJ)+DT*(GAHZCJl)*MFRCJ1)GAMZCJ)*
+ MFRCJ))*Y
MFS(J)MFS(J)+DT*(MFR(J1)MFRCJ))*Y
MTEMPMFSCJ)
STEMP=S(J)
IF(STEMP.LT.0.0001) STEMP=0.0001
RHOCJ )=MTEMP/STEMP
MOSCJ)MOS(J)+DT*(MOR(Jl)MOR(J))*Y
IFCMTEMP.LT.l.E9) MTEMP=l.E9
U(J)=M0S(J)/MTEMP
EFS(J)*EFS(J)+DT*(EFR(J1)EFR(J))*Y
P(J)((EFS(J)MOS(J)*U(J)/2.)/STEMP)*(GAM(J)1.)
Q(J)RHO(J)*D(J)*U(J)/288.
UF(J)=ABS(U(J))+SQRT(P(J)*GAM(J)/RHO(J))
GAM(J)=GAMD/MFS(J)
IF(GAM(J).LE.l.) GAM(J)=GATEMP
UFOAMAZ1(UF(J).UFO)
PI(J)(P(J)PAH)/144.
IQ1(Q1I(J)+Q(J))*DT/2.
IQ(J)IQ(J)+I
Q1I(J)0(J)
IlP(PlI(J)+PI(J))*DT/2.
IP(J)IP(J)+I1P
P1I(J)PI(J)
IFCJ.GE.JA1) GOTO 300
CALL FIGAM(P(J),RHO(J),GAM(J))
300 CONTINUE
C
DTDX/(TAD'UFO)
C
TT+DT
IF(T.GT.TLIM) GOTO 400
TPRT*1000.
GOTO 130
C
C
500 FORMAT(1H1,23X,24H* *PROGRAM BTUBE*****)
C
520 FORMAT(10X,4HDX= ,F10.4,5X,3HY= ,F10.4<,4X,4HS1= .E12.2/
+ 10X,4HS2= ,E12.2,3X,3HH= ,F10.4,4X,4HJD= .14/
+ 9X,5HJA1= ,14,9X,5HJA2= ,14,9X,5HJA3= ,14/
+ 10X,4HJN= ;14,7X,7HGAMA1= ,F8.4,3X,7HGAMA2= ,F8.4/
+ 8X,6HXMA1= ,F8.4,4X,6HXMA2= ,F8.4,4X,6HGAMT= ,F8.4/
+ 9X,5HPA1= ,F8.4,6X,4HTA= ,F10.4,3X,5HEPE= .E14.5/
+ 8X,6HTLIM= ,F8.4,5X,5HTAD= ,F8.4,4X,6HKOUT= ,14/
+ 7X,7HJNCK1= ,I4,7X,7HJNCK2= ,I4,11X,3HW= .F10.4//)
C
580 FORMAT(/6X,3HT .F8.2.4HMSEC)
C
590 FORMATCIX,2HP,IX,10(IX,F9.2))
C
591 FORMAT(IX,2HU=,IX,IOC IX,F9.2))
C
592 FORMATCIX,2HR,IX,10(IX,F9.6))
C
593 FORMATCIX,2HG=,IX,10(IX,F9.6))

Full Text 
PAGE 1
i COMPUTER SIMULATION OF SHOCK TUBES by John Randall Burris B.S., University of Colorado, 1983 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Department of Mechanical Engineering 1986
PAGE 2
This thesis for the Master of Science degree by John Randall Burris has been approved for the Department of Mechanical Engineering by Da te_3____::_;;k=...:...._[
PAGE 3
iii Burris, John Randall (M.S., Mechanical Engineering) Computer Simulation of Shock Tubes Thesis directed by Professor William Clohessy The use of computational analysis to model the flow in straight shock tubes has been made in the past, but variable diameter shock tubes. and nonideal shock tubes have not been considered in any detail. This study modifies an existing onedimensional, Eulerian, finite difference code for a straight shock tube, to match data obtained from an expanding diametershock tube and a nonideal shock tube. The finite difference scheme employed to solve the shock tube equation is called the Godunov method, which integrates the Euler equations for each cell using Green's theorm. A description of the nonideal precursor wave is also presented. Results obtained from the numerical computation agree well with the experimental data.
PAGE 4
iv ACKNOWLEDGMENTS I would like to take this opportunity to express my appreciation to Roy Heyman, for his valuable advice and guidance which made this project possible. And to the Martin Marietta Denver Aerospace company for it's financial support of this effort. Special thanks are also due to my committee chairman Dr. William Clohessy for his guidance and encouragement, and to Dr. Ken Ortega who has given me the inspiration to continue learning. I am indebted to Grady Romine for his understanding and support through this effort. Finally I wish to thank Sue, Beau, and Adrienne for their endless love and patience, and for making everything in life so worthwhile.
PAGE 5
v CONTENTS CHAPTER I. INTRODUCTION. . . . . . 1 II. PHYSICAL MODEL. . . . . . 4 Computer Code . . . . . 9 Expanding Nozzle Shock Tube....... 15 Precursor Shock'Tube.. .. . .. 18 III. GODUNOV METHOD...................... 24 IV. RESULTS. . . . . . . . 34 REFERENCES . . . . . . . . 4 2 APPENDIX A. Computer Listing and Subroutines.... 43 B. Expanding Nozzle Code............... 54 C. Precursor Code. . . . . . 60
PAGE 6
vi FIGURES Figure 2.1. Initial conditions pressure driven shock tube.................. 4 2.2. Flow in a shock tube after the diaphram is broken................. 5 2.3. Cell boundary...................... 12 2.4. Thunder tube facility.............. 16 2.5. Computer model of expanding cell... 17 2.6. Precursor wave development......... 19 2.7. Detailed precursor wave............ 21 2.8. Model for precursor shock tube..... 22 2.9. Configuration for precursor shock tube. . . . . . . 23 3.1. One dimensional model for computation. . . . . . 2 6 3.2. Diagram of solution................ 26 4.1. Experimental data from the Thunder tube. . . . . . 35 4.2. Experimental data from the nonideal shock tube............... 36 4.3. Computer generated data for the expanding shock tube........... 37 4.4. Computer generated data for the nonideal shock tube........... 38 4.5. Figure 4.3 imposed upon figure 4.1. 40 4.6. Figure 4.4 imposed upon figure 4.2. 41
PAGE 7
CHAPTER I INTRODUCTION Once thought to be doomed bythe computer, it is becoming evident that shock tubes are just as important as ever before. Someday will do almost all aerodynamic predicting. However,: computational fluid mechanics involves solution6f some version of the NavierStokes equation (in our case the Euler equation) This still requires input of some physical parameters which must be obtained from experiment. The passage of a high strength shock wave over an object produces loading on the object. Surface pressures and loads caused by the blast wave are highly transient and generate severe pressure gradients. Successful predictions of these loads is necessary if the designer of any blast resistant structure is to account for them. Such problems are encountered in the design of missle bases, space launch facilities, nuclear generating plants, and explosive or propellent manufacturing facilities. In the past, the study of shockinduced flow over bodies has been primarily done by experiment. The experiments involve large shock
PAGE 8
2 tubes or blast tubes and are extremely expensive to perform. To reduce test costs a computer program has been developed to simulate the shockinduced flow. The analysis can be employed to determine the explosive configuration, the optimum gas pressures, densities, and molecular weights, and the shock tube configuration to achieve a desired test condition. Although this idea of calculating blast waves of large amplitude in a shock tube is not new; variable area ducts and nonideal wave predictions have not previously been considered in any detail. This study then develops computer codes to model these two new applications for shock tube testing as well as general shock tube applications. The first application deals with the expansion in a shock tube. In order to accomadate a large scale model a shock tube located at Sandia Corporation in Albequerque was enlarged from a 6 foot diameter to a 19 foot diameter. The driver section, contairting the shock produ6ing explosives, has remained at a 6 foot diameter to reduce costs and explosive handling problems, while the test section has a 19 foot diameter. This configuration then introduces a 40 foot expansion section from the 6 foot driver section to the 19 foot test section. The second new application is in trying to model the precursor phenomena. A precursor or non ideal blast wave is formed when a
PAGE 9
3 nuclear weapon is exploded at some height over a thermally absorbing surface. The thermal radiation forms a hotair layer above the earth surface through which the blast front propagates at a higher velocity. The precersor wave or nonideal wave apparatus is a straight tube but it is unique in that various gases are placed in the tube in order to develop a simulated nonideal wave. Such facilities have been constructed at SRI International and at Kirkland Air Force Base. Test data for both of these tests showed large descrepencies with the qriginal computer predictions. The computer code required extensive modification to model. these new configurations and yield accurate predictions for future test data.
PAGE 10
CHAPTER II PHYSICAL MODEL Before we examine our method of solution let us develop the model for the physical situation of a one dimensional shock tube. Our blast tube facility can be idealized by the classical pressure driven shock tube shown in figure 2.1.a below. Driver section Driven section High Pressure Low Pressure G) G) (a) Pressure .... Distance (b) Figure 2.1 Initial conditions in a pressure driven shock tube A seperates the region of highpressure "driver" gas on the left (region 4) from the region of lowpressure
PAGE 11
5 gas on the right (region 1) The pressure distribution is shown in figure 2.1.b. The gases in regions 1 and 4 will generally have different temperatures and molecular weights. The pressure discontinuity that exists after rupture of the diaphram propagates into the test gas as a strong shock and into the driver gas as a centered rarefaction. in other words, when the diaphram in figure 2.1.a is broken a shock wave propagates into section 1 and an expansion wave propagates into section 4. This picture is sketched in figure 2.2. Expar:tsion wave propagating to the left Pressure Ill Contact surface moving at velocity of the gas behind the shock wave r G) Distance Normal shock wave propagating to the I right with velocity W .. 0 Figure 2.2 Flow in a shock tube after the diaphram is broken As the normal shock wave moves to the right with velocity
PAGE 12
6 W, it compresses the gas considerably thus resulting in high pressures and temperatures of the gas behind it (region 2) and induces a mass motion with velocity Up. The contact surface is defined as the interface between the driver and the driven gasses, moves with velocity Up. that the pressure and velocity are preserved across the contact surface: P 3 =P 2 and u 3=u2=Up, but the entropy changes discontinuously. From our figure it is seen that the expansion wave propagates to the left, smoothly and continuously decreasing the pressure in region 4 to the lower value P 3 behind the expansion wave. The flowfield in the tube after the diaphram is broken (figure 2. 2) is completely determined by the conditions in regions 1 and 4 before the diaphram is broken. For our blast tube there is no diap.hram and the driver section contains explosives whose sudden release of energy creates the high pressure of region 4 in figure 2.1.b. Besides this difference the mechanics of the flow field produced in the blast tube is similar to that produced in the pressure driven shock tube. Returning now to figure 2.1 we define the ratio P 4 /P 1 as the diaphram pressure ratio. Along with the initial conditions of the driver and driven. gas, P4/P1 determines the strength of the shock wave developed when the diaphram is removed (or as in our case a large
PAGE 13
7 pressure gradient is produced) As we know u 3=u2=up and P2=P3 across the contact surface. And it can be shown from the RankineHugoniot equations that the particle velocity induced behind the shock can be given as Also, since we are dealing with a simple centered expansion wave we know from the Riemann invarients 2a u + '= constant through the expansion wave. {.:_1 We can evaluate this constant by applying the equation to region 4 2a4 2a4 u4 +.= 0 +.=constant {1 {1 from which we obtain because a= V7JiT
PAGE 14
8 and because the flow is isentropic then If we now apply this equation between the head and tail of the expansion wave we have Solving for u3 gives Since u2=u3 we can equate these equations giving Which can be rearranged to yield
PAGE 15
9 P4 = P2l1 (1'41)(ada4) (P2/P1 1) l2y4h4 1 P1 P1 .J2n [211 + (1'1 + 1) (P2/P1 1)] This is the shock tube' equation which is used as a fundamental part of the Godunov method, as will be shown. This equation now gives the incident shock strength P 2/P1 as an implicit function of the diaphram pressure ratio P 4/P1 Because a= v0Jfi = V!(R/M) T the speed of sound in a light gas is faster than in a heavy gas. Since P 2 /P 1 will be stronger if a 1 I a 4 is made smaller, to maximize the incident shock strength for a given P 4/P1 the driver gas should be a low molecular weight gas at high temperature, and the driven gas should be a high molecular weight gas at low temperature. We are now ready to examine the method of solution. Computer Code The computer program used in this study is called BTUBE, which utilizes subroutines CELLPR, FIGAM, and AIR. BTUBE is written in Fortran IV and uses the finite difference technique called the Godunov Method (chapter III} The basis for the Godunov method is to integrate the Euler Equations for each cell using Green's theorm This will
PAGE 16
10 relate the volume changes in a cell of mass, momentum, and energy, to the flux of these quantities across the cell boundary. The shock tube equations are then used to evaluate the fluxes and determine the properties in each cell. This process is performed for each cell for each time step. Further insight is provided if we examine the path of the logic in the computer program BTUBE. Starting in the main program the dimensions and initial conditions of the cells are set up. Theprogram has the option of initializing cells containing different gases with varying properties. Another option is for the tube boundary conditions of either an opened or closed end. If the end is opened the gas is allowed to escape. For the closed end any incident wave is reflected back into the tube. Driver pressures are initialized by the weight of explosive in each cell. A time step is then calculated based on the length of the cell and the sonic speed of the gas. BTUBE then makes a call to subroutine CELLPR. This subroutine first determines which condition has occured in the cell; it determines if w1 and w2 are shock waves or expansion waves. Given the pressure, density, velocity, and values of gama across a cell boundary CELLPR then finds the mass, momentum, and energy fluxes across the cell boundary by solving_ the associated centered wave problem using the Godunov method. These boundary
PAGE 17
11 conditions are then used to solve the continuity, momentum, and energy equations for each cell by a finite difference scheme. Starting with the mass, momentum, and energy we have M a.s.s Momentum Energy Equations of State op 8 + (pu) = 0 {}t ox {} {} at (pu) +ox (p + pu2) = 0 {} {} 8t (pE) +ox (u (pE + p)) = 0 E = e + u2 /2 e=Pfp(tl) Now let = t2 t1 and make sufficiently small that the quantities p,u,P,E over S, the height of the cell, may be considered as constant. Consider the integral of these equations over and S. Recall that Greens theorm for a plane gives f ({}N oM) dx dy = f N dy + M dx } 5 oX 8Y Jc Where M and N are continuous singled valued functions over S bounded by the curve C and C is positive in the counter clockwise direction. Starting with equation 2.2.1 we have
PAGE 18
12 1 l f dt dz dy "' 1 l (p ),, (p ),, "' (pS),, (pS),, and taking S as constant over jr I 1t2 : (pu) dt dx dy jr I D.t!_ (pu) dx dy Js t1 a: }5 8x Then applying Greens theorm to obtain . lis (pu) dx = D.t L pu dy Figure 2. 3 shows the boundary C being made up of line segments 6 5 7 4 Yl s Y2 8 lr lr 3 1 2 Figure 2.3
PAGE 19
13 The boundary conditions are: along 12 and 56 dy=O.O along 23, 45, 67, and 81 U=O.O along 34, PU= p U= p U), and dy is positive along 78 P U), and dy is negative Then f pu dy = (pu )t Y2 (pu )t y1 Jc 2 1 And finally j ls J.:' (1) dt d" dy "' (pS),, (pS),, + Llt ((pu )2 y2 (pu )1 y1 ) = 0 or in another form Equations 2. 2. 2 and 2. 2. 3 can be integrated in similar fashion to obtain
PAGE 20
14 These are the finite difference equations which are then used by the main program. The conditions in each cell are now calculated for that time. The time is then increased the time step and the whole process is started over again. In summary CELLPR determines which shock condition has occured, solves the shock tube equation, and returns the fluxes to the main program. Subroutine FIGAM iterates with subroutine AIR to .. obtain the ratio of specific heats (gama) for the driver cells. Subroutine AIR provides a "semiphysical fit" to the equation of state which is generated from tabulated values. The explosive in the driver cells is modeled by inputing the energy of combustion and mass of the explosive into FIGAM which which then iterates with AIR to find the thermodynamic properties of the hot air. Thus, the amount of energy which is put into the gas by the explosive, is used to determine gama which is used for the equation of state. We can then use these values to determine the other propert.ies in the driver cells. Subroutines FIGAM and AIR are.used only for the driver cells, because the driver gets so hot. Gama is considered constant for the "driven" cells of our model. A complete computer listing of the program and its subroutines is given in appendix A.
PAGE 21
15 Expanding nozzle shock tube Now that we have described the general model of a shock tube we are ready to look at two applications which modify this model. First consider the Sandia Labs Thunderpipe facility in New Mexico which is illustrated in figure 2. 4. Both ends of the tube are open to the atmosphere. The tests are powered by strands of Iremite explosive placed inside the 6 ft. diameter "driver" section which is approximatly 156 ft. long. As a result of the energy release associated with the ignition of the explosive, a blast wave is created. The blast wave first travels thru the 40 ft. long expansion nozzle which increases the diameter from 6 ft. .to 19 ft. The shock then moves toward the test section located 156 ft. downstream. Although the tube is 19 ft. in diameter after expansion it contains a flat concrete roadbed which decreases the cross sectional area. An equivalent diameter, due to this road bed, was found to be 17.35 ft .. In terms of the physical model the main deviation from the simple shock tube is the smooth expansion nozzle, which is handled in a straight forward manner. The computer modeling of this geometry was accomplished by increasing the in steps. Through the expansion section each cells height is incrementally increased, eventually I
PAGE 22
Figure 2.4 Thunder tube facility
PAGE 23
17 stopping at the height of the driven section radius. So the computer code was rewritten to handle two situations, the straight sections and the expanding sections. The uniqueness of this analysis was in the expanding section. This approach explores the application of a onedimensional method to a two dimensional problem. From the finite difference equation for the continuity equation we have, MFS(J) MFS (J) + DT (MFR(J 1)MFR(J)) YI This is for the straight section of tube since the difference term is multiplied by a constant height of YI to obtain the mass flow. For the expanding section the mass of the previou. s cell (MFR (J1)) is multiplied by the height YI, and the mass flux of the present cell (MFR(J)) is multiplied by YI+DELY. This scheme is pictured in figure 2.5. EELY/2 r YI+DELY FLY/21 Figure 2. 5 Expanding cell geometry
PAGE 24
18 The momentum and energy equations were handled in exactly the same manner. Although this approach seems rather simple, the results agreed well with the experimental data as shown later. Precursor shock tube The second application we shall consider is the precursor phenomena. A precursor is the term used for an auxiliary blast wave caused by a nuclear explosion, which propogates ahead of the main shock wave due to the absorption of heat energy by a nonideal surface. If the intense thermal radiation from a nuclear detonation impinges on a thermally nonideal (heat absorbing) surface, a hot thermal layer is formed near the surface. This layer often includes smoke, dust, and other particulate matter, and is formed ahead of the blast since the thermal radiation travels faster than the airblast. The enhanced sound speed in the Thermal layer near the ground surface causes that portion of the incident blast wave which enters the Thermal layer to be accelerated. Thus, the shock wave propogating within the Thermal layer rapidly outruns the incident blast wave resulting in the formation of an aboveambient pressure region in the Thermal layer that still has unshocked, ambient pressure air lying above it. The thermal layer and it's interaction with the blast
PAGE 25
Dust Reflected blast wave Incident blast wave Preshock thermal layer Early Development Dust 19 Incident blast wave or Preshock thermal layer Late Development Figure 2. 6 Precursor wave development wave is shown in figure 2.6. Within the precursor region considerable modification of the usual blast wave characteristics may occur. From figure 2. 6 we see that after the precursor forms the main shock.front usually no longer extends to the ground. Between the bottom edge of the main shock wave and the ground is a gap, through which the energy that feeds the precursor may flow. Ahead of the main shock front, the blast energy in the precursor is free to follow the rapidly moving shock. front in the thermal layer and to propagate upward into the undisturbed air ahead of the main shock front. This diverging flow pattern within the precursor tends to weaken it, while the
PAGE 26
20 energy which is continually fed into the precursor from the main blast wave tends to strengthen the precursor shock front. Interest in the precursor wave comes about because the development of a precursor usually results in lower peak overpressures, increased rise times to peak overpressure, and increased dynamic pressure. The increased dynamic pressure can be attributed to higher particle velocities and/or greater air densities due to dust and other debris picked up by the blast wave as it passes over the ground. Another interesting phenomena which is created by the precursor (but not considered in this analysis) is the creation of vorticity. If we write the voticity equation for inviscid, compressible flow we have an ( ) 1 at + v x n x v P2 P x v P where Since the pressure gradient associated with the blast wave is nearly orthogonal to the density gradient in the Thermal layer underlying the ambient air, in general \7Px\7pf.O
PAGE 27
21 Thus, the formation of a vortex centered along the top of the Thermal layer with a strength which is p":i:oportional to both the shock strength and the Thermal layer density gradient is possible. This process is shown in figure 2.7. below, and is worthy of further study Reflect shock Vortex Incident shock Mach stem _) Hot thermal _Lr Figure 2. 7 Detailed precursor wave Two test models have been developed to simulate the dynamic and static pressure profiles of the precursor wave. The first employs a bag of helium which is placed on the floor of the shock tube to simulate the thermal layer, as shown in figure 2.8. The helium then acts as the hot thermal layer. This approach can accomodate a 1/30 scale test, in a 2 0 ft. diameter shock tube. But since the results from response tests are more accurate, and instrumentation is easier on a larger model, it is
PAGE 28
22 Shock wave .. I Helium Figure 2.8 Precursor shock tube desirable to increase the scale of our test. The second approach, which is used in this study, permits a 1/6 scale configuration in a 20 ft. diameter shock tube. This model loads the cells in the shock tube with varying gases having different molecular weights, starting with the heaviest and decreasing down the tube. This formation produces the desired one dimensional simulation of the precursor wave. A satisfactory configuration was found to be as shown on the next page in figure 2. 9. The first model develops the actual configuration of the precursor wave but requires an enormous diameter to accomodate larger scale testing. The second model does not physically develop the precursor wave, but it does simulate the dynamic and static pressure profiles experienced by the precursor. We now turn our attention to the finite difference method which was employed.
PAGE 29
Distributed explosive density SF 6 cells co 2 Test section He 23 Driven cells Figure 2. 9 Precursor shock tube
PAGE 30
CHAPTER III GODUNOV METHOD The problem considered here is the modification of a onedimensional finite difference program in order to simulate two applications of a shock tube. The finite difference method used is an Eulerian differencing scheme called the Godunov method. Godunov' s method uses the differential equations in. conservative form; as a result, the integration of these equations over the volume of a cell results in equating the changes in mass, momentum, and energy in the cell to the fluxes across the boundary of these quantities during the time step. Pressure will be considered as part of the momentum flux. At the start of a time step the flow properties are taken to be uniform in the cells on either side of a cell wall. As a result, the fluxes are given by the solution of the centered wave problem. Basically each cell can be visualized as a mini shock tube, and our problem narrows down to solving the shock tube equation for each cell. Our centered wave then has an imaginary diaphram at the cell boundary with conditions in the cell on that side of the boundary, and
PAGE 31
25 conditions on the other side of the diaphram equal to the conditions in the cell on the other side of the boundary, as shown in figure 3.1. At the instant that the diaphram is removed, two centered waves, w1 and w2 and a contact surface leave the origin as depicted in figure 3.2. The waves may be both shock waves, both expansion waves, or one a shock wave and one an expansion wave. The Eulerian boundary of each cell will lie in one of the four regions shown in figure 3. 2 for the entire timestep. Also, the values of the flow variables in these regions are taken as the values of these flow variables on the boundaries of. the cell. The main idea of the Godunov method is that of obtaining the fluxes on the boundary of a cell by solving the centered wave problem for that cell (as shown in figure 3. 2) Flow properties in the adjacent cells (i1/2, j+1/2) and (i+1/2, j+1/2) are known, and the flow properties on the boundary of the cell correspond to the properties on the vertical line X=O.O Now let us examine the equations of motion for our cell. We shall start with wave w1 Assuming it is a shock wave and applying conservation of mass and momentum we obtain (pV)I = (pV)u 3.1
PAGE 32
I pil J '+l 2' 2 pil j+l 2' 2 p+1 + 1 t 2 ,J 2 X uil i+ 1 2' 2 Figura 3.1 One dimensional model for computation II Contact surface J \ \ Pes \ Ucs \ \ III T W2 p.+1 '+1 t 2 2 Figura 3. 2 Diagram of solution IV 26
PAGE 33
27 3.2 Note that the subscripts denote regions in figure 3.2 V is the velocity of the fluid relative to the wave, P is the pressure, and P is the density. For simplicity we set the mass flux crossing wave w1 to bs so that equations 3.1 and 3.2 become 3.3 3.4 If Vw1 is the velocity of wave w1 and qn is the gas velocity normal to the cell boundary in figure 3. 2 we would have 3.5 Substitution of equation 3.5 into equation 3.4 yields 3.6 Through manipulation bs can be expressed in terms of P11
PAGE 34
and known values in region I b [ PI PI I ]1 12 "1/pil 1/PI we can now write and using the relationship 1 1 1 [ ( 1 + 1) PI + ( 1 1) PI I l PIIPI=PI (l1)Pr+(Fl)PII1 we have 1 [ 2 ( Pr PII) l = pI (I 1) Pr + (! + 1) PII 3.7 Substituting this back into equation 3.7 we have bs = [('Y + 1) PII + ('Y1) P1 ]112 2/pr 3.7a 28 Equations 3. 6 and 3. 7a are the final equations when considering wave w1 as a shock wave. Now let us suppose that wave w1 is an expansion wave; then the equations of unsteady isentropic flow will apply across the wave w1 It would be convenient to have
PAGE 35
29 these equations in the same form as equations 3.6 and 3.7a because this allows us to express the equality of pressure and velocity across the contact surface of our shocktube model. Since w1 is now a simple expansion wave propagating upstream we know that the Riemann invariant J+ is constant across w1 So we have Where C is the speed of sound substitute but 2 ,P c =p ( p == RT) p we have 2ytr [(p11 ') 1/2 (p1 ') 1/2] qnr qnrr 1 1PII, PI, (_1 ) 1/2 .. (:I ) 1/2y (_.!_) 1/2 PII, II, PI, 3. 8 If we 3.9 holds for an isentropic expansion. Substituting this into equation 3.9 yields = 2ytr [(Pn ) 1/2 ( p1 ) 1/2y (PI' ) 1/2] qnr qnrr 1 p IPI II. PI, 3.10
PAGE 36
30 We can now transform this equaiton into the same form as equation 3.5 by multiplying and dividing the righthand side of the equation by PrrPr and rearranging to give 3.11 Now if we define and substituting this relation into equation 3.11 gives 3.7b Which is the same form as equaiton 3.5. Thus equation 3.5 holds across wave w1 withb=bs (equation 3.7a) for a shock wave and with b=be (equation 3.7b) for an expansion wave. For wave w2 from figure 3.2, we obtain a similar set if. equations only with a change of sign. We can explain this change in sign in the following way. If wave w2 is a shock wave our equation becomes 3.12
PAGE 37
31 Where as is used for the mass flux crossing wave w2 If vw2 is the velocity of wave w2 3.13 Substitution of this value into equation 3.12 gives 3.14 Note that the sign has changed from equation 3.5 because wave w2 propogates downstream where wave w1 was propogating upstream. As seen from equation as is analogous to bs which is given from equation 3.7a as a,= [(I+ 1) Pu1 + (/1) P1vl112 2/PIV 3.15a Now if wave w2 is as expansion wave, the sign change occurs because the wave is a simple expansion wave propogating downstream, and the J_ Riemann invariant is constant across w2 Therefore the equation coriesponding to equation 3.8 is 3.16
PAGE 38
32 Also ae can be derived in the same way as be was to obtain ,1 ae=2! 1Prv [ 1Pur/ Prv l 1/ PIV 1(PuJ/ Piv)ryl/2"Y 3.1Sb Similar to w1 equation 3.14 always holds across wave w2 and equation 3. 1Sa holds if w2 is a shock wave while equation 3.1Sb holds if w2 is an expansion wave. If we now let and and solve equations 3.5 and 3.14 we obtain 3.17a 3.17b with a and b being the appropriate choice of as or ae and bs or be. These equations along with the applicable values for a and b are the basic equations for the exact method of Godunov. This set of equations requires an iterative solution. For example; assume a value for Pes; if Pcs>PI wave w1 is a shock wave and equation 3.7a gives the value of b. If Pcs Piv wave W2 is a shock wave and equation and equation 3.1Sa gives the
PAGE 39
33 value of a. And if Pcs
PAGE 40
CHAPTER IV RESULTS The success of any computer program is determined by how well it matches with the experimental data. Figures 4 .1 and 4. 2 show the experimental data obtained from the test section cif the expanding shock tube test and the nonideal shock wave test respectively. The computer results for these two applications are shown in figures 4.3 and 4.4. Comparison of the computer predictions to the test data show good for both cases. Note that for both cases the computer predicts an earlier arrival time than is found test. A timing scheme for the explosive was developed to try and simulate the burning of the explosive and possible delay the arrival time. The driver section was "turned on" cell by cell starting at the last driver cell and marching backwards, at the burn velocity of the explosive, until the end of the tube is reached. This method did delay the arrival of the shock some but not enough to compensate for the whole time gap. We believe this time descrepency could possible be explained by test procedures such as when the time starts
PAGE 41
o = rn U') Ul c: = rn :z: 72.00 OIL 00 50.00 Ye.oo YO.OO 32.00 211.00 10.00 0.00 o.oo 6.00 :. .. ..I .. .................................................... .............................. l ............................. ............................. r ........................... ............................. r ...................... :: l 1 1 : : :oo l : : .1. ; : .............. ................. l ............................... l'', : : ; : .................... ........ ; .... ,.! .. : l i : I : I: :1: :r:l t:::.t: ; I : : : : : : : : ....... .;,i',:,,. 1 ., r 1 r ... l ............................ :.i ... ... .. .. .. .. .. ... .. .. . ............................. : .................... 0 0 .................. i ............................. I I 1 I I ... ........................... ..... .. .. .. .. .. .. .. .. ........................... : .. .. ....... 1 ..... ..................... ............... l i ... ....................... \ o.oo TIME IN HSEC Figure 4.1 Experimental data from the Thunder tube w 01
PAGE 42
39 20 p R E 15 s s u R 10 E s p s I e ..... s 1e 39 48 .l .I '14 ... \, \4. y / \ /'v 50 69 80 9B tee 1.19 TIHE, i HSEC I Figure 4.2 Experimental from the nonideal shock tube 120 13G w 0\
PAGE 43
ao.o 70.0 50.0 40.0 JO.O a: 2.0.0 ro.o 0.0 10.0 0.0 2$.0 1\ \ \ \ I\ \ ............... :( 1"'"50.0 75. 0 fOO.O 1.25. 0 150.0 175.0 200.0 225.0 TIM (SEC) Figure 4.3 Computer generated data for the Thunder tube 2!SO.O w ....]
PAGE 44
30.0 25.0 20.0 15.0 '(I) 10.0 ::J (I) (/) Q 5.0 0.0 5.0 10.0 30.0 i : ; ( ( 40.0 I '\ .. ""'....... . 50.0 60.0 70.0 80.0 90.0 100.0 110.0 Time (msec) Figure 4.4 Computer generated data from the nonideal shock tube 120.0 .. 1 ) 130.0 w (X)
PAGE 45
39 with regards to when the explosive is ignited. Figures 4.5 and 4.6 show the computer plots superimposed upon the test data where the computer data has been shifted in time to account for the time gap. These figures show how closely the computer data matches the test data. So the computer programs utilizing the Godunov finite difference method has been shown to be very versatile in modeling the two applications studied in this analysis. Further work is hecessary to obtain even more accurate predictions and to explain the timing descrepency. However, the present model can be used with confidence in predicting the pressure loads experienced in testing with the expanding shock tube and the nonideal shock tube.
PAGE 46
., ::0 I'TI CJ') CJ') c: ::0 I"T1 72.00 "" 00 SI!I.OO 1&1!1.00 "o.oo 3Z.OO le.oo e.oo o.oo i i ............................... 1 I I'! r : r I I .... : ..... ... ............................. .. i i t ............................... : ............................... : .. : : : : : : :"''"'''''''i'''''''''''''''''''''''''''''''l''''''''''''''''''''''''''''''' j 1 : : : : .............................. ................. .............. i l i ............... ............................... .............................. i : : ........................... : ............................... : ............................. .. i : : : ; i :. : ...... .............................. .. : : : i : I I i : : : : '''''''''!oooo I I !.,_ r r o.oo sb.oo 0.01) 200.00 ____ TIME IN HSEC Figure 4.5 Figure 4.3 imposed upon figure 4.1
PAGE 47
39 ..29 .P R E 15 s s u R 10 E s 1 j ,. "' "'" \1 M p s I e \1 './ s 19 39 49 se 6e ;re ae 99 1 ee 1. 1 e 120 139 TIHE HSEC Figure 4.6 Figure 4.4 imposed upon figure 4.2
PAGE 48
REFERENCES J.D., Jr., Modern Compressible flow. With Historical perspective, McgrawHill, New York, 1982 Courant, R. and K.O. Friedrichs, Supersonic Flow and Shock Wayes, Interscience, New York, 1948 Gayden, A.G. and I.R. Hurle, The Shock Tube in High Temperature Chemical Physics, Reinhold, New York, 1963 Glasstones, S. and P. Dolan (ed.), The Effects of Nuclear Weapons, U.S. Department of Defense, 1977 Liepmann, H.W. arid A. Roshko, Elements of Gasdynamics, Wiley, New York, 1957 McNamara, W., Flame Computer Code for the Axisymmetric Interaction of a Blast Wave with a Shock Layer on a Blunt Jour. of Spacecraft, Vol.4, No.6, June 1967 Shapiro, A.H., The Dynamics and Thermodynamics of Compressible Fluid Flow. Vol. I, Ronald, New York, 1953 Zucrow, M.J. and J.D. Hoffman, Gas Dynamics, Whiley, New York, 1977
PAGE 49
APPENDIX A.
PAGE 50
c c c PROGRAM BTUBE(INPUT,OUTPUT,STAT,DYN,STIMP,DYNIMP, + TAPE5aiNPUT,TAPE6=0UTPUT,TAPE7=STAT,TAPEB=DYN, + TAPE9=STIMP,TAPE10=DYNIMP) C PROGRAM EMPLOYS GODUNOV FINITE DIFFERENCE C METHOD (SEE W. MCNAMARAJ.SPACECRAFT AND C ROCKETSJUNE 1967PP.790795) TO COMPUTE C THE PRESSURES INSIDE A DUCT. c c ******************** c *************** C THIS VERSION EMPLOYS C VARIABLE GAMMA USING C OOANNICKLE SUBROUTINE c ... c ...... c c c c c c c c REAL MFR(600),MOR(600),MFS(600),MOS(600),MFRS,MORS,MTEMP, + DIMENSION EFR(6.00), P(600), PI(600) ,RH0(600), U(600), EFS(600), + S(600).GAM(600),GAMZ(600),Q(600),UF(600),W(200), + P1I(600),Q1I(600),II(10) COMMON/UNIBK/PA1,RH1,AS1,PAH COMMON /UNBBK/ J ARK, MFRS MORS, EFRS, GAMC COMMON/AIR1/LGAM,TEMP DATA + 49753.7978,3.141592654 / C READ(7,700) DX,Y,S1,S2,H C READ(7,710) JD,JA1,JA2,JA3,JN C READ(7,720) GAMA1,GAMA2,XMA1,XMA2,GAMT,EPE C READ(7,730) PA1,TA C READ( 7. 740) TLIM, TAD, KOUT. JNCK1 JNCK2 C READ(7,750,END=5) W c c c c NAMELIST /BANG/ DX,Y,S1,S2.H,JD,JA1,JA2,JA3,JN,GAMA1,GAMA2 + XMA1,XMA2,GAMT,EPE,PA1,TA,TLIM,TAD; KOUT,JNCK1,JNCK2,W,II READ(5,BANG) WRITE(6,BANG) c c c 5 PAHPA1*144. LGAM0 JTERM=JN1 JLP=(JN2)/10 TAR=TA+TR ASl=SQRT(GAMT*TAR*RG) RH1=PAH/(RG*TAR) EPA=PAH/((GAMTl.)*RHl) RGAl=RGO/XMAl 44
PAGE 51
c c c c c c RHA1PAH/(RGA1*TAR) RGA2RGO/XMA2 RHA2PAH/(RGA2*TAR) UFOAS1 S(1)S1 S(JN)S2 P(l)PAH P(JN)PAH RHO(l)RH1 RHO(JN)RH1 PI(l)0.0 PI(JN)0.0 U(l)=O.O Q(l)0.0 U(JN)O.O Q(JN)0.0 GAM( 1) GAMT GAM(JN)GAMT GAMZ(1)GAMT GAMZ(JN)=GAMT MFS(l )RHl"S(l) MFS(JN)=RH1*S(JN) MOS(l)0.0 MOS(JN)0.0 EFS(1)=P(l)*S(1)/(GAM(1)1.) EFS(JN)=P(JN)S(JN)/(GAM(JN)1.) MFR(l)=O.O MFR(JN)0.0 MOR(l)=O.O MOR(JN)0.0 EFR(l)0.0 EFR(JN)O.O IP(l)=O.O IQ(l)O.O Pli(l)0.0 Qli(l)0.0 IP(JN)=O.O IQ(JN)=O.O Pli(JN)O.O Qli(JN)0.0 AEXPI"H*H/4. WRITE(6,500) WRITE(6,520) DX.Y,Sl,S2,H.JD,JAl,JA2,JA3,JN, + GAMAl,GAMA2.XMA1,XMA2.GAMT.PAl. TA,EPE,TLIM. + TAD,KOUT,JNCKl,JNCK2,W(2) WRITE(6,600) 00 10 K=2, JTERM 10 S(K)=DX*Y 00 100 J=2,JTERM IF(J.GT.JD) GOTO 60 RHOEW(J)/(GC'AE'DX) EAVG=(EPA*RH1+EPE'RHOE)/RHO(J) RTEMPRHO(J)*0.51579185 GM=0.4 45
PAGE 52
c c c CALL AIR(ETEMP,RTEMP,GM) GAM(J)=l.+GM P(J)GM*EAVG*RHO(J) GAMZ(J)=GAM(J) PI(J)(P(J)PAH)/ 144 GOTO 90 60 IF(J.LT.JA1) GOTO 80 IF(J.GE.JA2) GOTO 75 P(J)=PAH PI(J)O.O GAM(J)=GAMA1 GAMZ(J)=GAMA1 GOTO 90 75 IF(J.GE.JA3) GOTO 80 P(J)=PAH PI(J)=O.O RHO(J)=RHA2 GAM(J)=GAMA2 GAMZ(J).;,GAMA2 GOTO 90 80 P(J)PAH PI(J)=O.O RHO(J)RH1 GAM(J)=GAMT GAMZ(J)=GAMT 90 MFS(J)=RHO(J)*S(J) MOS(J)=O.O EFS(J)=P(J)*S(J)/ (GAM(J)1.) Q(J)0.0 MFR(J)=O. O MOR(J)=O. O EFR(J)=O.O IP(J)=O.O IQ(J)=O.O Pli(J)=O.O Q1I(J)=O.O UF(J)=ABS(U(J))+SQRT(P(J)*GAM(J)/RHO(J)) UFO=AMAX1(UF(J),UFO) 100 CONTINUE DTDX/(TAD*UFO) T=O.O TPR=O.O KOUNT=O 130 KOUNT=KOUNT+1 WRITE(7 ,1000) (TPR, (PI(II(I)), I=l, 10)) WRITE(8,1000) (TPR,(Q(II(I)),I=l,lO)) WRITE(9,1000) (TPR,(IP(II(I)),I=l,10)) WRITE(l0,1000) (TPR,(IQ(II(I)),I=l,10)) IF(TPR.LT.0.0001) GOTO 132 IF(KOUNT.LT.KOUT) GOTO 150 132 WRITE(6,580) TPR KOUNT=O KL=1 DO 145 K1, JLP KK=10*K WRITE(6,590) (PI(I),I=KL,KK) 134 WRITE(6,591) (U(I), I =KL,KK) 46
PAGE 53
c c 136 WRITE(6,592) (RHO(I),I=KL,KK) 138 WRITE(6,593) (GAM(I),I=KL,KK) 140 WRITE(6,594) (Q(I),I=KL,KK) 142 WRITE(6,595) (IP(I),I=KL,KK) 143 WRITE(6,596) (IQ(I),I=KL,KK) 144 WRITE(6,600) 145 KL=KK+1 WRITE(6,610) 150 DO 160 J=1,JTERM IF(JNCK1.LT.1) GOTO 154 IF(J.NE.1) GOTO 154 P(J)=P(J+1) RHO(J)RHO(J+1) U(J)=U(J+1) GAM(J)=GAM(J+1) GOTO 156 154 IF(JNCK2.LT.1) GOTO 156 IF(J.NE.JN1) GOTO 156 P(J+1)=P(J) RHO(J+1)=RHO(J) U(J+1)U(J) GAM(J+1)=GAM(J) 156 CALL CELLPR(P(J),RHO(J),U(J),GAM(J),P(J+1),RHO(J+1), + U(J+1),GAM(J+1)) GAMZ(J)=GAMC MFR(J)MFRS MOR(J)=MORS 160 EFR(J)=EFRS c c UFO=AS1 DO 300 J=2,JTERM GATEMP=GAM(J) GAMD=GAM(J)*MFS(J)+DT*(GAMZ(J1)*MFR(J1)GAMZ(J)* + MFR(J))Y MFS(J)=MFS(J)+DT*(MFR(J1)MFR(J))*Y MTEMP=MFS(J) STEMP=S(J) IF(STEMP.LT.0.0001) STEMP=0.0001 RHO(J)=MTEMPISTEMP MOS(J)=MOS(J)+DT*(MOR(J1)MOR(J))*Y IF(MTEMP.LT.1.E9) MTEMP=1.E9 U(J)=MOS(J)/MTEMP EFS(J)=EFS(J)+DT*(EFR(J1)EFR(J))*Y Q(J)=RHO(J)*U(J)*U(J)/288. UF(J)=ABS(U(J))+SQRT(P(J)*GAM(J) / RHO(J)) GAM(J)=GAMD/MFS(J) IF(GAM(J).LE.l.) GAM(J)=GATEMP UFO=AMAX1(UF(J),UFO) PI(J)=(P(J)PAH)/144. IQ1=(Q1I(J)+Q(J))*DT/2. IQ(J)=IQ(J)+IQ1 Q1I(J)=Q(J) I1P=(P1I(J)+PI(J))*DT/2. IP(J)=IP(J)+IlP P1I(J)=PI(J) IF(J.GE.JA1) GOTO 300 CALL FIGAM(P(J),RHO(J),GAM(J)) 300 CONTINUE 47
PAGE 54
c c c c TmT+DT IF(T.GT.TLIM) GOTO 400 TPR=T*1000. GOTO 130 500 FORMAT(1H1,23X,24H*****PROGRAM BTUBE*****) c 520 FORMAT(10X,4HDX= ,F10.4,SX,3HY= ,F10.4,4X,4HS1= ,E12.2/ + 10X,4HS2= ,E12.2,3X,3HH= ,F10.4,4X,4HJD= ,I4/ + 9X,SHJA1= .I4,9X,SHJA2= .I4,9X,SHJA3= ,I4/ + 10X,4HJN= ,I4,7X,7HGAMAl= ,F8.4,3X.7HGAMA2= ,F8.4/ + 8X,6HXMA1= ,F8.4,4X,6HXMA2= ,F8.4,4X,6HGAMT= ,F8.4/ + 9X,SHPA1= .F8.4,6X,4HTA= .F10.4.3X,5HEPE= ,El4.5/ + 8X,6HTLIM= .F8.4.SX,SHTAD= ,F8.4.4X,6HKOUT= ,I4/ + 7X,7HJNCK1= ,I4,7X,7HJNCK2= ,I4,11X,3HW= ,F10.4//) c 580 FORMAT(/6X,3HT= ,F8.2,4HMSEC) c C 590 .. 10(1X,F9.2)) 591 FORMAT(lX,2HU=,lX,10(1X,F9.2).) c 592 FORMAT(lX,2HR=,lX,lO(lX,F9.6)) c 593 FORMAT(lX,2HG=,1X,10(1X,F9.6)) c 594 FORMAT(lX, 2HQ= ,1X,l0(1X,F9. 2)) c 595 FORMAT(1X,3HIP=,l0(1X,F9.2)) c 596 FORMAT(1X,3HIQ=,l0(1X,F9.2)) c 600 FORMAT(//) c 610 FORMAT(////) c 700 FORMAT(2F12.4,2El2.2,Fl2.4) c 710 FORMAT(SI4) c 720 FORMAT(SF10.4,El4.5) c 730 FORMAT(2F10.4) c 740 FORMAT(2F10.4,3I4) c 750 FORMAT(7F10.4) c 1000 FORMAT (11(1X,E10.4)) c c c c 400 WRITE(6,600) END SUBROUTINE CELLPR(PA,RHOA,UA,GAA,PB,RHOB,UB,GAB) REAL MA,MB,MFRS.MORS COMMON/UNIBK/PAl,RHl,ASl,PAH 48
PAGE 55
c c c FHL(PTE)mPHHPTE/((1.+SWT*BEQ*(UHUL)/AHRBEP* + (PTE1.)*SQRT(1./(BLH*PRR*(ALP*PTE+1.))))**(1./BEP)) PCK2=ABS(PBPAH) EA=PA/((GAA1.)*RHOA) EB=PB/((GAB1.)*RHOB) UAA=ABS(UA) UAB=ABS ( UB) IF(PCKl.LT .. 01.AND. PCK2.LT .. O,l.AND. UAA. LT .. OLAND. + UAB.LT .. 01) GOTO 100 PAB=PAIPB IF(PAB.GE.1.0) GOTO 10 PH=PB RHOHRHOB UH=UB GAMH=GAB SWT=1.0 PL=PA RHOLRHOA UL=UA GAMLGAA PHH. /PAB GOTO 20 10 RHOH=RHOA UH=UA GAMH=GAA SWT=l.O PLPB RHOL=RHOB GAML=GAB PHH=PAB 20 PRR=PH*RHOL/(PL*RHOH) ALR=SQRT(GAML*PL/RHOL) AHR=SQRT(GAMH*PH/RHOH) ALP=(GAML+1.)/(GAML1.) BEP=(GAMH1.)/(2.*GAMH) BEQ=GAMH*BEP BLH=(GAMLl.)/(2.*GAMH) IF(PRR.LT.200.) GOTO 21 PF2=PHH*(1.+SWT*BEQ*(UHUL)/AHR) GOTO 26 21 IF(PRR.GT .. 001) GOTO 22 PF2=1.0001 GOTO 26 22 IP=O PFl=PHH PF2=PHH/2. 23 FP1=FHL(PF1) FP2FHL(PF2) IF(ABS(FP2).LE .. 02) GOTO 26 IF(ABS(FP2FPl).LT.l.E9) GOTO 26 IP=IP+l IF(IP.LT.20) GOTO 24 WRITE(6,800) FPl,FP2 800 FORMAT(l0X,20HNO CONVERGENCEFPl= ,F12.3,2X,5HFP2= ,Fl2.3) GOTO 26 49
PAGE 56
24 PF3=PF2FP2*(PF2PFl)/(FP2FPl) PFl=PF2 PF2PF3 GOTO 23 26 PSFPF2 IF(PSF.LT .l.E9) PSF=l.E9 PII=PSF*PL IF(SWT.GE.O .O) GOTO 30 PI2=PII IF(PII/PB.GT.0.999) PI2=0.999*PB + *BEP) GOTO 45 30 PI2=PII IF(PII/PA.GT.0. 999) PI2=0.999*PA + **BEP) 45 URA=UA URB=UB VA=UAMA/RHOA VB=UB+MB/RHOB PCSc(MA*PB+MB*PA+MA*MB*(URAURB))/(MA+MB) IF(PCS.LT.l.E9) PCS=l.E9 UCS=(PAPB+MA*URA+MB*URB)/(MA+MB) RHOA2RHOA*((GAA+l.)*PCS+(GAAl.)*PA)/((GAAl.)*PCS+ + (GAA+l.)*PA) RHOB3RHOB*((GAB+l.)*PCS+(GABl.)*PB)/((GABl.)*PCS+ + (GAB+l.)*PB) IF(VB.GT.d.O) GOTO 60 PCPB UC=UB GAMC=GAB JAHKl GOTO 90 60 IF(UCS.GT.O.O) GOTO 70 PC=PCS RHOC=RHOB3 UC=UCS IF(EB.LE.l.E+9.0R.GAB.GT.l.4) GOTO 65 GAMR=GAB CALL FIGAM(PC.RHOC,GAMR) GAMC=GAMR GOTO 67 65 GAMC=GAB 67 JAHK=l GOTO 90 70 IF(VA;GE.O .O) GOTO 80 PC=PCS UC=UCS IF(EA.LE.l.E+9.0R.GAA.GT.l.4) GOTO 75 GAMR=GAA CALL FIGAM(PC,RHOC,GAMR) GAMCGAMR GOTO 77 75 GAMC=GAA 77 JAHK=2 GOTO 90 80 PC=PA 50
PAGE 57
c c c c c RHOC=RHOA UC=UA GAMC=GAA JAHK=2 90 MFRS=RHOC*UC MORS=PC+MFRS*UC EFRS=(GAMC*PC/(GAMC1.)+MFRS*UC/2.)*UC RETURN 100 MFRS0.0 MORS=O.O EFRS=O.O JAHK=O GAMC=GAA RETURN END SUBROUTINE FIGAM( PG, RHOG, GAMG) C THIS SUBROUTINE ASSUMES THAT PRESSURE CALL IS IN C UNITSOF LBF/FT2 AND DENSITY CALL IS IN UNITS. OF C LBFSEC2/FT4. IT THEN ITERATES TO FIND THE VALUES C OF GAMMA AND TEMPERATURE APPROPRIATE TO THOSE C CONDITIONS. c PMPG*478.8 RH01=RHOG*0.51579185 RHOM=RHOl GM=GAMG1.0 GMOLD=GM DO 10 I=1,40 GM=(GM+GMOLD)/ 2. GMOLD=GM EE=PM/(GM*RHOM) IF(EE.LE.2.E+13) GOTO 5 EE=2.E+13 RHOM=PM/(GM*EE) 5 CALL AIR(EE.RHOM.GM) IF(ABS((GMGMOLD)/GM).LT.1.E5) GOTO 20 RHOMRH01 IF(I.LT.40) GOTO 10 WRITE(6,600) GM,GMOLD GOTO 20 10 CONTINUE c c c c c c c c 20 GAMG=GM+1.0 RETURN 600 FORMAT(10X,22HNO CONVERGENCEGMOLD= ,F9 6,2X, + 4HGM= ,F9.6/ ) END SUBROUTINE AIR(EEE,RRR,GMONE) C THIS SUBROUTINE ADAPTED FROM L R .DOAN AND G .H.NICKEL C "A SUBROUTINE FOR THE EQUATION OF STATE OF AIR" 51
PAGE 58
C AFWL RTD(WLR)TM632.MAY 1963WITH UPDATES BY C C.NEEDHAM,ET.AL.IN THE DNA 1KT STANDARD c C EEE IS ENERGY(ERGSIGM) C RRR IS DENS!TY(GM/CM3) C GMONE IS c c c c c c c c c c c c c c COMMON/AIR1/LGAM,TEMP E1=(8.5E)/,975 IF(ABS(E1).LT.5.0) GOTO 20 IF(E1.GT.O.O) GOTO 10 FO=O.O WS=O.O GOTO 30 10 FOEXP(E/4.46) WS=l.O GOTO 30 20 DE10.975*RH0**0.05 EE1=8.5+0.155042*ALOG(RHO) E1=(EE1E)/DE1 WS=1./(EXP(E1)+1.) FOEXP(E/4.46)*WS FONEXP(E/6.63)*(1.WS) 30 BETA=O.O IF(E.LT.1.0) GOTO 40 E2=(E40.)/3.0 IF(ABS(E2).LT.5.) GOTO 60 IF(E2.GT.O.O) GOTO 50 40 FN=O.O WS=O.O GOTO 70 50 FN=EXP(E/25.5) WS=1.0 GOTO 70 60 DE2=4.*RH0**0.085 EE2=45.*RH0**0.0157 E2=(EEE2)/DE2 WS1./(EXP(E2)+1.) 70 E3=(E160.)/6.0 FE=O.O IF(E3.GT.(5.0)) FE=l./(EXP(E3)+1.) GMONE=(0.161+0.255*F0+0.28*FON+0.137*FN+0 05*FE)* 52
PAGE 59
c c c c + RHO* *BETA IF(I.GAM.LT.1) RETURN RHOLNaALQG(RHO) RHOLS=RHOLN*RHOLN GMGMONE PaGM*EEE*RRR IF(E.GT.120.) GOTO 100 F=0.9722 .7143E3*RHOLN+6.582549E5*RHOLS G0.026456.418873E4*RHOLN2.338785E5*RHOLS H9.21E5+5.971549E6*RHOLN+3.923123E7*RHOLS CON1=3480.*GM*E .CON2=F+G+H*E*E TEMP=CON1/CON2 GOTO 200 100 ALR=ALOG10(RRR) C1=1.69081E5+2.99265E7*ALR C2=0.769813+3.8618E3*ALR TEMP=Cl*EEE**C2+102.6275735 200 BETT=(E3) / 0 66 IF(BETT.GT 10 ) RETURN TEMP=C0Nl/(CON2*(1.SWIT)+SWIT) IF(TEMP.GT. O.O) RETURN ALRALOG10(RRR) C1=1.69081E5+2.99265E7*ALR C2a0.769813+3 8618E3*ALR TEMPaC1 *EEE* *C2+ 102 6275735 RETURN END 53
PAGE 60
APPENDIX B.
PAGE 61
c c c c c c c c c c c c c c c c c c c c c c c c c c PROGRAM TBTUBE PROGRAM EMPLOYS GODUNOV FINITE DIFFERENCE METHOD (SEE W. MCNAMARAJ.SPACECRAFT AND ROCKETSJUNE 1967PP.790795) TO COMPUTE THE PRESSURES INSIDE A DUCT. THIS VERSION EMPLOYS VARIABLE GAMMA USING DOANNICXLE SUBROUTINE REAL MFR(600),MOR(600),MFS(600),MOS(600),MFRS,MORS,MTEMP, R IP(600),IQ(600),I1P,IQ1 DIMENSION EFR(600),P(600),PI(600),RH0(600),U(600),EFS(600), D S(600),GAM(600),GAMZ(600),Q(600),UF(600),W(140), D Pli(600),Qli(600),II(lO) COMMON/UNIBK/PAl,RHl,ASl,PAH COMMON/UNBBK/JAHK,MFRS,MORS,EFRS,GAMC COMMON/AIRl/LGAM,TEMP DATA GC,TR,RG,RGO,XPI/32.2,459.67,1717.6027, D 49753.79741,3.141592654/ READ(7,700) DX,S1,S2,H READ(7,710) JD,JA1,JA2,JA3,JN READ(7,720) GAMA1,GAMA2,XMA1,XMA2,GAMT,EPE READ(7,730) PA1,TA READ( 7, 740) TLIM, TAD, KOUT, JNCK1 JNCK2 READ(7,750) Y1,YF,JEXPS,JEXPI READ(7,760,END) W NAMELIST /BANG/DX,S1,S2,H,JN,JD,JA1,JA2,JA3,GAMA1,GAMA2, N XMA1,XMA2,GAMT,EPE,PA1,TA,TLIM,TAD,KOUT,JNCK1,JNCK2, N Y1,YF,JEXPS,JEXPI,W,II READ(5,BANG) WRITE ( 6, BANG) 5 PAHPA1. l.GAMO JTERMJN1 JLP(JN2)/10 TARzTA+TR ASlzSQRT(GAMT*TAR*RG) RH1PAH/(RG*TAR) 55
PAGE 62
EPAPAH/((GAMT1.)*RH1) UFOaAS1 S(l)81 S(Jl0S2 P(l) .. PAH P(JN)PAH RHO(l)RR1 PI(l)0.0 PI(JN)O.O U(l)0.0 Q(l)0.0 U(JN)O.O Q(JN)0.0 GAM( 1) GAMT GAM( JN) GAMT GAMZ( 1) GAMT GAMZ ( JN) GAMT MFS(l)RH1*S(l) MFS(JN)RH1*S(JN) MOS(l)0.0 MOS(JN)O.O EFS(1)P(1)*S(1)/(GAM(l)1.) EFS(JN)P(JN)*S(JN)/(GAM(JN)1.) MFR(l)0.0 MFR(JN)O.O MOR(l)0.0 EFR(l)0.0 EFR(JN)O.O IP(l)0.0 IQ(l)O.O P1I(l)O.O Qli(1)o.o IP(JN)EO.O IQ(JN)O.O Pli(JN)O.O Qli(JN)0.0 AEXPI*H*H/4. WRITE(6,500) WRITE(6,510) Y1 WRITE(6,5ll) YF WRITE(6,512) JEXPS WRITE(6,513) JEXPI WRITE(6,520) DX,S1,S2,H,JD,JA1,JA2,JA3;JN, + GAMA1,GAMA2,XMA1,XMA2,GAMT,PA1,TA,EPE,TLIM, + TAD,KOUT,JNCK1,JNCK2,W(2) WRITE(6,600) YIYl JEXPS1JEXPS1 DO 10 K=2,JEXPS1 10 S(K)DX*YI YYI DELY(YFYI)/JEXPI 56
PAGE 63
JEXPE1JEXPE1 00 15 KJEXPS,JEXPE1 YY+DELY S(K)DX*Y 15 CONTINUE 00 20 KJEXPE, JTERM 20 S(K)DX*Y 00 100 J,JTERM IF(J.GT.JD) GOTO 80 RHOEW(J)/(GC*AE*DX) RHO(J)RHOE+RH1 EAVG(EPA*RH1+EPE*RHOE)/RHO(J) ETEMPEAVG*929.02267 RTEMPRHO(J)*0.51579185 GM.4 CALL AIR(ETEMP,RTEMP,GM) GAM(J).+GM P(J)GM*EAVG*RHO(J) GAMZ(J)GAM(J) PI(J)(P(J)PAH)/144. GOTO 90 80 P(J)PAH PI(J)O.O RHO(J)RH1 GAM( J) GAMT GAMZ ( J)GAMT 90 MFS(J)RHO(J)*S(J) MOS(J)O.O EFS(J)P.(J)*S(J)/(GAM(J)1.) U(J)O.O Q(J)0.0 MFR(J)O.O MOR(J)O.O EFR(J)O.O IP(J)O.O IQ(J)O.O Pli(J)O.O Q1I(J)O.O UF(J)ABS(U(J))+SQRT(P(J)*GAM(J)/RHO(J)) UFQAMAX1(UF(J),UFO) 100 CONTINUE DTDX/(TAD*UFO) TO.O TPRO.O KOUNTO 130 KOUNTKOUNT+1 YIYl WRITE(7,1000) (TPR,(PI(II(I)),I,10)) WRITE(8, 1000) (TPR, (Q(II(I)), I, 10)) WRITE(9,1000) (TPR,(IP(II(I)),I1,10)) WRITE(10,1000) (TPR,(IQ(II(I)),I1,10)) IF(TPR.LT.0.0001) GOTO 132 IF(KOUNT.LT.KOUT) GOTO 150 132 WRITE(6,580) TPR KOUNTO KL 00 145 K,JLP 57
PAGE 64
XK=10*K WRITE(6,590) (PI.(I),I=KL,KK) WRITE(6,591) (U(I),IaKL,KK) WRITE(6,592) (RHO(I),I=KL,KK) 138 WRITE(6, 593) (GAM( I), I=KL, KK) 140 WRITE(6,594) (Q(I),I=KL,KK) 142 WRITE(6,595) (IP(I),I=KL,KK) 143 WRITE(6, 596) (IQ(I), WRITE(6,597) (S(I),I=KL,KK) 144 WRITE(6, 600) 145 KLKK+ 1 WRITE(6,610) 150 DO 160 J1,JTERM IF(JNCK1.LT.1) GOTO 154 IF(J.NE.1) GOTO 154 RHO(J)RHO(J+1) U(J)U(J+1) GAM(J)GAM(J+1) GOTO 156 154 IF(JNCK2.LT.1) GOTO 156 IF(J.NE.JN1) GOTO 156 P(J+1)P(J) RHO(J+1)":"RHO(J) U(J+1)U(J) GAM(J+1)GAM(J) 156 CALL CELLPR(P(J),RHO(J),U(J),GAM(J),P(J+1),RHO(J+1), + U(J+l) ,GAM(J+1)) GAMZ(J)=GAMC MFR(J)MFRS MOR(J)MORS 160 EFR(J)EFRS UFOAS1 JSTART=2 JENDJEXPS1 210 DO 200 JJSTART, JEND GATEMP=GAM(J) I. GAMDGAM(J)*MFS(J)+DT*(GAMZ(J1)*MFR(J1)GAMZ(J)* + MFR(J))*YI MFS(J)MFS(J)+DT*(MFR(J1)MFR(J))*YI IF(STEMP.LT.0.0001) STEMP=O.OOOl RHO(J)MTEMP/STEMP MOS(J)MOS(J)+DT*(MOR(J1)MOR(J))*YI IF(MTEMP.LT.l.E9) MTEMPl.E9 U(J)MOS(J)/MTEMP EFS(J)EFS(J)+DT*(EFR(J1)EFR(J))*YI P(J)((EFS(J)MOS(J)*U(J)/2.)/STEMP)*(GAM(J)1.) Q(J)RHO(J)*U(J)*U(J)/288. UF(J)ABS(U(J))+SQRT(P(J)*GAM(J)/RHO(J)) GAM(J)=GAMD/MFS(J) IF(GAM(J).LE.l.) GAM(J)=GATEMP PI(J)(_P(J)PAH)/144.
PAGE 65
IQ1m(Q1I(J)+Q(J))*DT/2. IQ(J)2IQ(J)+IQ1 Qli(J)Q(J) I1Ps(P1I(J)+PI(J))*DT/2. IP(J)IP(J)+IlP Pli(J)=PI(J) IF(J.GT.JD) GO TO 200 CALL FIGAM(P(J),RHO(J),GAM(J)) 200 CONTINUE IF(J.GE.JTERM) GO TO 220 DO 300 J=JEXPS,JEXPE GATEMP=GAM(J) GAMDGAM(J)*MFS(J)+DT*(GAMZ(J1)*MFR(J1)*YIGAMZ(J)* + MFR(J)*(YI+DELY)) MFS(J)=MFS(J)+DT*(MFR(Jl)*YIMFR(J)*(YI+DELY)) MTEMP=MFS(J) STEMPS(J) IF(STEMP.LT.O.OOOl) STEMP20.0001 RHO(J)aMTEMP/STEMP MOS(J)MOS(J)+DT*(MOR(J1)*YIMOR(J)*(YI+DELY)) IF(MTEMP.LT.l.E9) MTEMPzl.E9 U(J)MOS(J)/MTEMP EFS(J)EFS(J)+DT*(EFR(J1)*YIEFR(J)*(YI+DELY)) P(J)((EFS(J)MOS(J)*U(J)/2,)/STEMP)*(GAM(J)1.) UF(J)ABS(U(J))+SQRT(P(J)*GAM(J)/RHO(J)) GAM(J)GAMD/MFS(J) IF(GAM(J).LE.1.) GAM(J)GATEMP UFOAMAX1(UF(J),UFO) PI(J)m(P(J)PAH)/144. IQl(Q1I(J)+Q(J))*DT/2. IQ(J)=IQ(J)+IQ1 Qli(J)Q(J) IlP=(Pli(J)+PI(J))*DT/2. IP(J)=IP(J)+IlP P1I(J)PI(J) YIYI+DELY IF(J.GT.JD) GO TO 300 CALL FIGAM(P(J),RHO(J),GAM(J)) 300 CONTINUE JSTART=JEXPE+1 JEND=JTERM GO TO 210 220 CONTINUE DT=DX/(TAD*UFO) T=T+DT IF(T.GT.TLIM) GOTO 400 TPRsT*1000. GOTO 130 500 BTUBE) 510 FORMAT(10X,18HINITIAL ,F8.2/) 511 FORMAT(10X,16HFINAL DIAMETER= ,F8 .2/) 512 FORMAT(l0X,26HEXPANSION STARTS AT CELL= ,I4/) 513 FORMAT(10X,30HNUMBER OF STEPS IN EXPANSION= ,I4/) 59
PAGE 66
APPENDIX C.
PAGE 67
c c C PROGRAM BTUBE c c C PROGRAM EMPLOYS GODUNOV FINITE DIFFERENCE C METHOD (SEE W. MCNAMARAJ. SPACECRAFT AND C ROCKETSJUNE 1967PP.790795) TO COMPUTE C THE PRESSURES INSIDE A DUCT. c c c C THIS VERSION EMPLOYS C VARIABLE GAMMA USING C DOANNICKLE SUBROUTINE c c REAL MFR(600),MOR(600),MFS(600),MOS(600),MFRS,MORS,MTEMP, R IP(600),IQ(600),I1P,IQ1 DIMENSION EFR(600),P(600),PI(600),RH0(600),U(600),EFS(600), D S(600),GAM(600),GAMZ(600),Q(600),UF(600),W(300), D P1I(600),Q1I(600),II(10) COMMON/UNIBK/PA1,RH1,AS1, PAH COMMON/UNBBK/JAHK,MFRS,MORS,EFRS,GAMC COMMON/AIR1/LGAM,TEMP DATA GC, TR,RG,RGO.XPI/32.2,459.67,1717 6027, D 49753.7978,3.141592654/ NAMELIST /BANG/DX,Y,S1,S2,H,JD,JA1,JA2,JA3,JA4,JN,GAMA1,GAMA2, N GAMA3,XMA1,XMA2,XMA3,GAMT,EPE,PA1,TA.TLIM,VEXP,TAD, N KOUT,JNCK1,JNCK2,W,II READ( 5, BANG) WRITE(6,BANG) 5 LGAM=O JTERMJN1 JLP=(JN2)/10 TAR=TA+TR AS1SQRT(GAMT*TAR*RG) RH1=PAH/(RG*TAR) EPAPAH/((GAMTl.)*RHl) RGA1RGO/XMA1 RHA1=PAH/(RGAl*TAR) RGA2aRGO/XMA2 RHA2=PAH/(RGA2*TAR) RGA3=RGO/XMA3 RHA3=PAH/(RGA3*TAR) UFOcAS1 S(l)=S1 61 l
PAGE 68
. P(JN)PAH RH0(1)RH1 RHO(JN)RH1 PI(l)0. 0 U(l)0.0 Q(l)=O.O U(JN)0.0 GAM( 1) =GAMT GAM(JN)GAMT GAMZ(1)=GAMT GAMZ (JN) =GAMT MFS(l)=RH1*S(l) MOS(l)=O.O MOS(JN)O.O EFS(1)P(1)*S(1)/(GAM(1)1.) EFS(JN)=P(JN)*S(JN)/(GAM(JN)1.) MFR(l)=O.O MFR(JN)O.O MOR(1).0 MOR(JN)=O.O EFR(l)O.O EFR(JN)0.0 IP(l)O. 0 IQ(l) .. O.O P1I(l)O.O Qli(l)0.0 IP(JN)=O.O IQ(JN)0.0 P1I(JN)=O.O Q1I(JN)=O. 0 AEXPI*H*H/4. WRITE(6,500) WRITE(6,S20) DX,Y,S1,S2,H,JD,JA1,JA2,JA3,JN, + GAMA1,GAMA2,XMA1;XMA2,GAMT,PA1,TA,EPE,TLIM, + TAD,KOUT, JNCK1, JNCK2, W(2). WRITE(6,600) 00 10 K=2, JTERM 10 S(K)=DX*Y DO 100 J=2,JTERM 60 IF(J.LT.JA1) GOTO 80 IF(J.GE.JA2) GOTO 75 P(J)PAH PI(J)O.O RHO(J)RHA1 GAM(J)GAMA1 GAMZ(J)=GAMA1 GOTO.eo 75 IF(J.GE.JA3) GOTO 77 P(J)=PAH 62
PAGE 69
PI(J)=O.O RHO(J)=RHA2 GAM(J) .. GAMA2 GAMZ(J)GAMA2 GOTO 90 77 IF(J.GE.JA4) GOTO 80 P(J)PAH PI(J)O.O GAM(J)GAMA3 GAMZ(J)GAMA3 GOTO 90 80 PI(J)O.O RHO(J)=RH1 GAM( J) GAMT GAMZ ( J) cGAMT 90 MFS(J)mRHO(J)*S(J) u(J)o.o Q(J)0.0 MFR(J)0.0 MOR(J)O.O EFR(J)O.O IP(J)O.O IQ(J)mO,O Pli(J)=O.O Q1I(J)=O.O UF(J)=ABS(U(J))+SQRT(P(J)*GAM(J)/RHO(J)) UFO=AMAX1(UF(J), UFO) 100 CONTINUE DT=DX/(TAD*UFO) TO.O TPR=O.O JSJD XREM=O.O KOUNT=O 130 KOUNT=KOUNT+l XST=(VEXP*DT)+XREM IF(JS.LE.2)GO TO 20 30 IF(XST.GE.DX)THEN IF(JS.LE.2)GO TO 20 RHO(JS)RHOE+RHl EAVG=(EPA*RHl+EPE*RHOE)/RHO(JS) ETEMPEAVG*929.02267 RTEMP=RHO(JS)*0.51579185 GM=0.4 CALL AIR(ETEMP,RTEMP,GM) GAM(JS)=l.+GM P(JS)=GM*EAVG*RHO(JS) GAMZ(JS)=GAM(JS) 63
PAGE 70
PI(JS)(P(JS)PAH)/144. MFS(JS)=RHO(JS)*S(JS) EFS(JS)=P(JS)S(JS)/(GAM(JS)1.) UF(JS)ABS(U(JS))+SQRT(P(JS)*GAM(JS) / RHO(JS)) JSEJS1 XSTXSTDX IF(XST.GT.O.O)THEN XREMmXST ELSE XREM=O.O' END IF ELSE XREMDXXST END IF IF(XREM.GE.DX) GO TO 30 20 CONTINUE WRITE(7,1000) (TPR,(PI(II(I)),I=1,10)) WRITE(8,1000) (TPR,(Q(II(I)),I=1,10)) WRITE(9,1000) (TPR,(IP(II(I)),I=1,10)) WRITE(10,1000) (TPR,(IQ(II(I)),I=1,10)) IF(TPR.LT.0.0001) GOTO 132 IF(KOUNT.LT.KOUT) GOTO 150 132 WRITE(6,580) TPR KOUNTO KL=1 DO 145 K,JLP KK*K WRITE(6,590) (PI(I),I=KL,KK) 134 WRITE(6,591) (U(I),I=KL.KK) 136 WRITE(6,592) (RHO(I),I=KL,KK) 138 WRITE(6,593) (GAM(I),I=KL,KK) 140 WRITE(6,594) (Q(I),I=KL,KK) 142 WRITE(6,595) (IP(I),I=KL,KK) 143 WRITE(6,596) (IQ(I).I=KL,KK) 144 WRITE(6,600) 145 WRITE(6,610) 150 DO 160 J,JTERM IF(c1NCK1. LT .1) GOTO 154 IF(J.NE.1) GOTO 154 P(J)P(J+1) RHO(J)RHO(J+1) GAM(J)=GAM(J+1) GOTO 156 154 IF(c1NCK2.LT.1) GOTO 156 IF(J.NE.JN1) GOTO 156 P(J+1)P(J) RHO(J+1)RHO(J) U(J+1)U(J) GAM(J+1).GAM(J) 156 CALL CELLPR(P(J).RHO(J),U(J),GAM(J),P(J+1),RHO(J+1), + U(J+1),GAM(J+1)) GAMZ(J)=GAMC MFR(J)MFRS MOR(J)=MORS 64
PAGE 71
I 160 EFR(J)EFRS c c c c c c UFOASl DO 300 J,JTERM GATEMPGAM(J) GAMDzGAM(J)*MFS(J)+DT*(GAMZ(J1)*MFR(J1)GAMZ(J) + MFR(J))*Y MFS(J)MFS(J)+DT*(MFR(J1)MFR(J))*Y MTEMPMFS(J) STEMP=S(J) IF(STEMP.LT.0.0001) STEMP=0.0001 RHO(J)zMTEMP/STEMP MOS(J)MOS(J)+DT*(MOR(J1)MOR(J))*Y IF(MTEMP.LT.1.E9) MTEMP=1.E9 U(J)=MOS(J)/MTEMP EFS(J)EFS(J)+DT*(EFR(J1)EFR(J))*Y P(J)((EFS(J)MOS(J) *U(J)/2. ) I STEMP) (GAM(J)1.) Q(J)RHO(J) *U(J) *U(J) /288. OF(J)=ABS(U(J))+SQRT(P(J)*GAM(J)/RHO(J)) GAM(J)GAMD/MFS(J) IF(GAM(J).LE.1.) GAM(J)=GATEMP OFOAMAX1(0F(J),UFO) PI(J)(P(J)PAH)/144. IQ1(Q1I(J)+Q(J))*DT/2. IQ(J)IQ(J)+IQ1 Qli(J)Q(J) I1P(P1I(J)+PI(J))*DT/2. IP(J)IP(J)+IlP Pli(J)PI(J) IF(J.GE.JAl) GOTO 300 CALL FIGAM(P(J),RHO(J),GAM(J)) 300 CONTINUE DTDX/(TAD*UFO) TT+DT IF(T.GT;TLIM) GOTO 400 TPRT*lOOO. GOTO 130 500 FORMAT(1Hl,23X,24H*****PROGRAM BTUBE*****) c 520 FORMAT(lOX,4HDX= ,Fl0.4,5X,3HY= :,Fl0.4,4X,4HS1= ,El2.2/ + 10X,4HS2= El2.2,3X,3HH= ,Fl0.4,4X,4HJD= ,!4/ + 9X,5HJA1= ,I4,9X,5HJA2= ,I4,9X,5HJA3= ,I4/ + ;I4,7X,7HGAMA1= ,F8.4,3X,7HGAMA2= ,F8.4/ + 8X,6HXMA1= ,F8.4,4X,6HXMA2= ,F8.4,4X,6HGAMT= ,F8.4/ + ,F8.4,6X,4HTA= ,Fl0.4,3X,5HEPE= ,El4.5/ + 8X,6HTLIM.;. ,F8.4,5X,5HTAD= ,F8.4,4X,6HKOUT= ,I4/ c + 7X,7HJNCK1= ,I4,7X,7HJNCK2= ,I4,11X,3HW= ,Fl0.4//) 580 FORMAT(/6X,3HT ,F8.2,4HMSEC) c 590 FORMAT(lX,2HP,lX.lO(lX,F9.2)) c 59l FORMAT(lX,2HU=,lX,lO(lX,F9. 2)) c 592 FORMAT(lX,2HR,lX, 10(1X,F9.6)) c 593 FORMAT(lX,2HG=,lX,l0(1X,F9.6)) 65

