- Email: [email protected]
An expansion equation of state subroutine
COMPUTER PHYSICS COMMUNICATIONS 5 (1973) 64—68. NORTH-HOLLAND PUBLISHING COMPANY
AN EXPANSION EQUATION OF STATE SUBROUTINE K. MORGAN* UKAFA, A WRE, Aldermaston, Berkshire, UK
Received 1 September 1972
PROGRAM SUMMARY Title of program: EOSEXP Catalogue number: ACSA Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland (see application form in this issue). Computer: ICL ATLAS 2, Installation: UKAEA, AWRE, Aldermaston Operating system: ATLAS Supervisor Programming language used: FORTRAN High speed storage required: 2600 words. No. of bits in a word: 48 Is the program overlaid? No No. of magnetic tapes required: None What other peripherals are used? Card Reader; Printer No. of cards in combined program and test deck: 424 Card punching code: BCD Keywords: Fluid Dynamics, Equation of State, Change of State. Nature of the physical problem Hydrodynamic calculations of physical systems require suitable equations of state for the materials of interest. EOSEXP is an expansion equation of state subroutine which can be used to give a realistic description of the single-phase and two-phase behaviour of an expanding material, provided that conditions along the saturation curve for this material can be tabulated.
*
Present address: Department of Mathematics, University of Exeter, Exeter, Devon, UK.
Method of solution In the single-phase region a Gruneisen type equation of state is used to obtain the pressure given the specific internal energy and the specific volume. In the two-phase region an implicit equation is solved for the temperature, and the pressure is then obtained by using the relationship between pressure and temperature valid along the saturation curve. Typical running time On the AWRE ATLAS loading and compilation of the test deck takes 17 sec. The running time depends on the number of times the equation of state subroutine is entered. Typically the calculation of 10 state points (5 single-phase and 5 two-phase) requires a CPU time of 4 sec.
K. Morgan, An expansion equation of state subroutine
65
LONG WRITE-UP I Introduction
(v
.
In many hydrodynamic calculations it is necessary to consider the behaviour of materials which from an initial cold state (at specific volume p = v0) are subject to increase in internal energy and expand. The pressure—volume (p, v) plane for v > v0 is divided into two by the saturation curve (see fig. 1). Above this curve the material exists as a single-phase, either liquid or gas, whilst below the saturation curve both phases co-exist. This paper describes an expansion equation of state subroutine EOSEXP which can be used to cover that region in the (p,v) plane for which vthe > v0, and which provides a realistic description of single-phase and two-phase behaviour of the material. The form used for the equation of state is based on a method first proposed by Morris [1].
, _______z~___.________.___._._.____._
—
V~)I(Vg—
v2)
=
(e
—
VQ)I(Cg
—
e2).
(3)
This is an implicit equation defining p (or 7) in terms of v and e, which is just the type of equation needed for hydrodynamic calculations. 3. Behaviour in the single-phase region In the single-phase region the behaviour of the material may be described by an equation of state of GrUneisen type. This is of the form 0~ [e er(V)] Iv, (4) Pwhere = PrO’) ~r er are assumed or determined functions Pr +and along some reference curve. With some prescribed variation for ~r this allows analytical expansion away from this reference curve to cover the remainder of the —
(p,v) plane in question. vopo,,
I
/ /~‘~~t.d hq.d
4. Computer The computer subroutine subroutine EOSEXP calculates a valuespecific and for the internal pressure,energy, p, given e. It specific is clearvolume, from the v,
~
SPECIFIC VOWME,v
Fig.!. Subdivision of the pressure—volume plane for v >
vo.
2. Behaviour in the two-phase region In the two-phase region the specific volume v of the material may be written as V =
av~+ (1
—
(1)
a)v2,
and the specific internal energy as e=
~g + (1
—
a)eQ,
(2)
a is the mass fraction of vapour the mixture, v2 and Vg are the specific volumes and e2 and eg the specific internal energies of saturated liquid and gas at the vapour pressure p and temperature T, i.e., they are conditions on the saturation curve. Eliminating a from eqs. (1) and (2) gives .~
preceding sections that if conditions along the saturation curve can be determined then the behaviour of the material in the region v > v 0 can be completely described, for eq. (4) can be used with the saturation curve as reference curve in the single-phase region and eq. (3) can be solved in the two-phase region. Corresponding values of V1.j, Pri’ en, r’ri and Tn (1 ~ I ~ KTBL) are tabulated at a series of points along the saturation curve and input as data to the computer. For each data point on the liquid side of the saturation curve the corresponding point at the same temperature and pressure on the vapour side must also be input. The method used at AWRE for inputting this data on cards is illustrated in table l,which shows the saturation curve data for water. The first card specifies the number (KTBL) of data points on the saturation curve. There then follows KTBL cards, read in order of increasing ~‘ri’each one specifying vri, Pri’ en, Fri and Tn (1 e~i ~ KTBL). This data
66
K. Morgan, An expansion equation of state subroutine Table 1 The method of inputting saturation curve data on cards (shown for water) Specific volume 3/g) (cm
CARD 1. CARD 2 CARD 3 CARD 4 CARD 5 CARD 6 CARD 7 CARD 8 CARD 9 CART) to CARD 11 CARD 12 CARD CARD CARD CARD CARD CARD CARD CARD CARl)
13 j4 15 16 17 18 19 20 21.
Pressure (Mb)
Specific internal energy
Gamma
Temperature (°K)
(cm”Mb/g)
20 ~.,00OO0~.E00 t,043200~ 00 ~..U732ooE 00 1,127300E 00 3.,216300E 00 1,4~21ooE 00 1,611500E 00 1,769000E 00 2,039000E (JO ~,80O000E 00 i 1470000E 00 ~,8530ooE 00 8,3790~0E 00 1,1491.50E 01 1,802210E 01. 6,66163001 ~.,94~273E 02 O,060732E 02 t,694043E 03 1,29201O~ 05
1.300000E—08 0, 10,000000E—07 4,18O0O0~—03 3,000000EwO6 5,61.0000E03 i0,000000E06 7,620000E—o3 3.000000E!,05 1,004D00~.02 1o,oooooo~~05 1,393000E—o2 i.400000E04 1,~48O0QE.o2 1.,700000E—04 1,660000E—02 2.000000E—04 1,7B6000E..o2 2,209000E—04 1,959000E—02
i.981000E—01 4.410000E.Oi. 5.110000E01 5.760000EO1 ~~5~0Q0~E—0i 4.483O00~.o1 2.811Q00~—01 4.027000Em01 2.755000E’Ol 1.82300OE’~O1
2.209000E’04 2.000000E04 1.,700000E—o4 1,40UO00E~04 1o,000000E~05 3,000000E—05 1(l,000000EO6 3,O0000OE~06 lo.000000E—07 1.3OU000E’08
i.89li00O~N01 6,470000E 02 2.732000E”Ol 6,410000E 02 3.263~Q~E.0~.6,27OOO~E 02 2.829000E~01 6,110000E 02 2.530000E01 5,860000E 02 3,000000E’.Ol. 5,o80000E 02 3.000000B—01. 4,540000E 02 3.000000E01. 4,070000E 02 3.000000E01 3.730000E 02 3.000000E01. 2,803000E 02
2,0~8000E—02 2,294000Ee02 2,4o56ooE~o2 2,476700E—02 2.544600E-o2 2,602400E—o2 2,583900~—o2 2,545400E—a2 2.,.505500E.02 2,385600E—02
2.803000E 02 3,730000F 02 4,070000E 02 4,540000F 02 5.080000E 02 5,860000E 02 6.1i.0000E 02 6,270000E 02 6~41.o0ooE 02 6,470000E 02
Table 2 The method of storing the data of table! in array EOSTAB for use with subroutine EOSEXP N
1 2
3 4 5
6 7 S 9
10
11 12
13 14 15
16 17 18
19 20 21
22 23 24
25
EOSTAB(N)
N
EOSTAB(N)
N
EOSTAB(N)
2.000000E 01 1.~O~)O00E00 1.i143200E 00 1.~?3200E 00
3.oc0000E—05 i.000000E—04 1.400000E—04 1,700000E—04 2.o00000E—04
00
26 27 28 29 30 31 32
51 52 53 54 55 56 57
00
33
00
34
2.11~90OUE 00 2.8~’’)O0OE 00 3.4?0000E 00 5.~53000E 00 8.3?9000E 00 l,149150f~01 1.8fl2210E 01 ~.661630E 01 1.943273E 02
35
1.9SQ000E—02 76 2.u~8000E—02 77 2.254000E—02 78 2.405600E02 79 2.476700E—02 80 2.544600E—02 81. 2.oo2400E—02 82 2.5e~900E—02 83 2.545400E02 84 2,5C5500E—02 85 2.365600E02 86 1.9~1000F—O1 87 4.410000E01 88 5.il0000EOl 89 5,7e~OO00E—(11 90 5.5(i00U0F~01 91 4.483000E—01 92 2.811000E-01 93 4.fl27000E—01 94 2.7~5000E01 95 1.8231J00E01 96 1.8Y)(100E01 97 2.7~20(i0E01 98 3~263000E01 99 2.R29000E—0l 100 101.
i.~~7300E 1.216300E 1.452100E a.611SOOE 1.7~9000E
00
00
36 37 38 39 40 41 42 43 6.’60732E 02 44 i.6c404~3P 03 45 1.292010E 05 46 1.~0000E—08 47 i.nC0000E—06 48 3.flfJ0000E—06 49 1.i~t)000OE—05 50
2.2ij9000E—04 2.2c~.0o0E—04 2.n00000E—04
1.700000E—04 1.400000E04 1.000000E—04 3.000000E—05 1.000000E—05 3.000000E—06 1.o00000E—06 1.300000E—08 0. 4.:l80000E—03 5.610000E—03
7.620000E—03 1.oo4000E—02 1.3c3000E—o2 1.548000E—02 1.6e0000E—02 1.786000F—02
56
59 60 61. 62 63 64 65 66 67 68 69 70 71 72 73 74
75
N
EOSTAB(N) 2.530000E01 3.1100000E01 3.~t~0000E01
3.~~~000E01 3.)IJ0000E01 3~UO0O0E—01 2.,3u3000E 02 3.730000E 02 4.(37’lOOOE 02 4.540000E
02
5.OS0000E 02 5.8t~.0000E 02 6.110000E 02 6.270000E 02 6.41(I000E 02 h.470Ii00E 02 ~.47’1000E 02 6.410000E 02 6.270000fE 02 6.110000E 02 5.860000E 02 5.e0000E 02 4.540000E 02 4.070000~ 02 3.7~0O00E 02 2.803000E 02
K. Morgan, An expansion equation of state subroutine
67
must be fed into the array EOSTAB in the form required for use with the expansion equation of state
cific volume v by interpolating linearly along lines of constant pressure, i.e., for each data point!, a1 is cal-
subroutine. The first store, EOSTAB (IV), in the equation of state table for a particular material contains the number KTBL for that material. Stores EOSTAB (N + 1) to EOSTAB (N + 5KTBL) are then allocated as follows: EOSTABt’ri’ (N + 1) to EOSTAB (N + KTBL) contam the EOSTAB (N + KTBL + 1) to EOSTAB (N + 2KTBL)
culated, where = e + (e. e.) (v v .)/(v. v.) (7) / / I J £ j Here i = n + 1 j, and n is the number of data points input. This table is then searched to determine the interval (~,~_!‘ which e is with given (v,e) thenej)liesinbetween thelocated. isobarsThe (andpoint isotherms) / and!, i.e., p. ‘~p~p.andT. ~
containthepri, EOSTAB(N+2KTBL+l)toEOSTAB(N+ 3KTBL) contain the en, EOSTAB (N + 3KTBL + I) to EOSTAB (N + 4KTBL) contain the Fri, EOSTAB (N + 4KTBL + 1) to EOSTAB (N + 5KTBL) contain the Tn.
—
—
—
—
—
‘
j-1
/
(6) The quantities eQ, eg,vQ and vg used in eq. (3) are expressed in the form e (7) = e. + (e. e. —
/
~
j-!
eg(T) = e~+ 1+ (e~ e,~1)Z,
(8)
—
This process is illustrated in table2which shows how the data of table 1 is stored in array EOSTAB. Values of quantities between the tabulated points are obtained by interpolation. The method of interpolation used in the subroutine assumes that on the liquid side of the saturation curve Tr, er and Fr vary quadratically with ~r’and that Pr may be obtained from a relationship of the form 1 = A ‘T + B (5\ °~P~
1, r
“ ‘
!
On the vapour side of the saturation curve en and Tr are taken to vary linearly with Tr, and Vr and T~are assumed to be related by logv~A2/T~+B2.
(T) = V11 + (v1 v1 1)Z, v (T) = exp [log v. + A log v/v .÷ VQ
—
g
l
1
1
where Z = (T T,_ ~)/(~‘j 7~~) and A = (liT l/T1)/( lIT1 l/T1_ i)~ In this way eq. (3) may be regarded as an equation for Z and the solution for Z in the range 0 ~ Z ~ 1 is determined by the method described by Peters and Wilkinson [2]. A print of the iterations resulting from the application of this method can be obtained by setting the value 1.0 to the common variable EXIM. When Z (and hence T), has been determined, the pressure is obtained by using eq. (5). —
—
—
—
(6)
Eq. (5) is again applied to relate pressure and tempera-
5. Test runs
tune. The method of calculation in the subroutine may be split into sections (identifiable on the listing) as follows: (1) The tabulated specific volumes are searched to determine the interval in which v is located, (2) The interpolated quantitiesPr’ en, Fr and Tr are then determined as described above. (3) A check is made on the sign of(e er). (4) If(e en) > 0 then the material exists as a single-phase and p is determined by applying eq. (4). (5) If(e er) <0 then both phases co-exist, and eq. (3) is solved to determine p. The saturation curve data is used to set up a table of energy values at spe-
This expansion equation of state subroutine is used in the TOPAL hydrodynamic code which was developed at AWRE to investigate the problem of molten metal/liquid thermal interactions. The subroutine has performed satisfactorily for a variety of materials and Hoskin and Morgan [3] bave described results obtainèd using the subroutine for calculations on water and also on sodium. Special driver routines have been provided to enable test runs to be perfqrmed. Subroutine PRELUDE allocates storage for the 4ommon variables and the main routine TEST reads in the saturation curve data, a specific volume anda specific internal energy.
—
—
—
68
K. Morgan, An expansion equation of state subroutine
Subroutine EOSEXP then calculates the corresponding pressure and, depending on whether the calculated state point lies in the single or two-phase region, sets the value I or 2 to the common variable IMAR. The allocation of storage in this deck is peculiar to ATLAS 2. If the program is to be run on another machine then storage allocation should be specified in a manner suitable to that installation. This may entail modifying or removing PRELUDE and altering or replacing the COMMON and DIMENSION statements in the other routines*.
References Ii] E. Morris, private communication (1969). 121 G. Peters and J.H. Wilkinson, Computer J. 12 (1969) 398. 13] N.E. Hoskin and K. Morgan, Proc. Intern. Conf. Engineening of fast reactors for safe and reliable operation, Karlsruhe (1972), to be published.
TEST RUN OUTPUT
EXFA~SIO~EnS DATA FO~ WATER 20 E~Tt~IESFC!~ TARLE VCL(ME 1.0000uOE 10 1.1J43200E 00 1.07321OE 00 1.12?3uOE 00 1.2i.~3i0E (it) 1.4~21~’OE00 1.6i.15tDE 00 1.7690uOE 00 2.o3c000E 00 2.80001iOE 00 3.4701!(OE 00 5.8530t’OE JO 8.37couoE 00 1.14c1~0E ul 1.802210E 01, 6.6616~OE (Ii. 1.943273E 02 6.060732E 02 1.694043E 03 1.292010E ;5
V
2.78100Th 2.c2s000E 3.09000Th 3.274000’ 3.470000E 3.98100Th 6.31000Th 1.00000Th 1.00000Th 1.00000Th
01
Oil Or; 00
on 00 00 (11
02 03
OF’
CONDITIONS ON SATU8ATICN CLRVE
FI~ESSURE ).300000E—08 1(.000000E—07 ~-000000E—06 i0.000000E—06 ~.000QOOE—05 lU.000000E—115 1.400000E—fl4 1.700000E04 2.000000E—04 2.209000E—04 2.209000E—04 2.000000E—04 1.700000E—04 1.400000E—04 ~o-C00000E—05 ~.000000E05 1C.000000E—06 ~.O00OQ0E—06 io.IJ00000E—07 1.3o0000E—0~
E 2.075700E—02 2.071600E—02 2.067100E—02 2.062500E—02 2.(158000E—02 2.047300E—02 2.OoelOOE-02 1.964900E02 1.727600E02 1.5i~500E—02
ENFRC~Y 0. 4.18~jOQflE—03 5.hl0000F—03 20000E.03 7.6 1.OO4000E—02 1.393000E—02 1.548000F—02 1.660000F02 1.786000E—02 1.959000E—02 2.058000E—02 2.294000E—02 2.405600E—02 2.476700E—02 2.544i500E—02 2.602400E—02 2.583Q00E—02 2.545400E—02 2.505500E—02 2.385600E—02
P PHASE 2.998971E—04 1 2.799685E—04 I 2.600073E—04 I 2.4001.60E—04 1 2.209000E-’04 1 2.004294E—i)4 2 1.428224E—04 2 9.739631E—O~ 2 1.032795E—05 2 8.910063E—07 2
fiAMhA 1.981000E—Q1. 4.4j,U000E—tIl. 5.ll0000E—01. 5.76(’OOOE—Ol. ~.5.J’C0U0E—OI. 4.48300UE—01 2.811000E01 4.027000E01 2.755000E01 1.823000E—01 ~.89i:00UE~O1 2.732000E—01 3.263000E—O1 2.829OnuE—01 2.530000E—01 3.Ou0000E—01, ~.O0~000E—D1 ~.0flhb0tt0E.01 3.Ou0000E—01 30000’u’OE—01
T~ItPENATURE 2.Sfl3000E 02 3,73000Th
02
4.070000E 4,5400000202 ~.O~30000E 02 ~.8SO000E 02 6,11000Th 02 6,2?0U00E 02 6.41000Th 02 6.470000E 02 6.470000E 02 6.410000E 02 6.270000E 02 6.110000E 02 ~.860O00E 02 ~.080000E 02 ~,~4O00OE 02 4.07000Th 02 3,730000E 02 2.8)3000E 02
AN EXPANSION EQUATION OF STATE SUBROUTINE K. MORGAN* UKAFA, A WRE, Aldermaston, Berkshire, UK
Received 1 September 1972
PROGRAM SUMMARY Title of program: EOSEXP Catalogue number: ACSA Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland (see application form in this issue). Computer: ICL ATLAS 2, Installation: UKAEA, AWRE, Aldermaston Operating system: ATLAS Supervisor Programming language used: FORTRAN High speed storage required: 2600 words. No. of bits in a word: 48 Is the program overlaid? No No. of magnetic tapes required: None What other peripherals are used? Card Reader; Printer No. of cards in combined program and test deck: 424 Card punching code: BCD Keywords: Fluid Dynamics, Equation of State, Change of State. Nature of the physical problem Hydrodynamic calculations of physical systems require suitable equations of state for the materials of interest. EOSEXP is an expansion equation of state subroutine which can be used to give a realistic description of the single-phase and two-phase behaviour of an expanding material, provided that conditions along the saturation curve for this material can be tabulated.
*
Present address: Department of Mathematics, University of Exeter, Exeter, Devon, UK.
Method of solution In the single-phase region a Gruneisen type equation of state is used to obtain the pressure given the specific internal energy and the specific volume. In the two-phase region an implicit equation is solved for the temperature, and the pressure is then obtained by using the relationship between pressure and temperature valid along the saturation curve. Typical running time On the AWRE ATLAS loading and compilation of the test deck takes 17 sec. The running time depends on the number of times the equation of state subroutine is entered. Typically the calculation of 10 state points (5 single-phase and 5 two-phase) requires a CPU time of 4 sec.
K. Morgan, An expansion equation of state subroutine
65
LONG WRITE-UP I Introduction
(v
.
In many hydrodynamic calculations it is necessary to consider the behaviour of materials which from an initial cold state (at specific volume p = v0) are subject to increase in internal energy and expand. The pressure—volume (p, v) plane for v > v0 is divided into two by the saturation curve (see fig. 1). Above this curve the material exists as a single-phase, either liquid or gas, whilst below the saturation curve both phases co-exist. This paper describes an expansion equation of state subroutine EOSEXP which can be used to cover that region in the (p,v) plane for which vthe > v0, and which provides a realistic description of single-phase and two-phase behaviour of the material. The form used for the equation of state is based on a method first proposed by Morris [1].
, _______z~___.________.___._._.____._
—
V~)I(Vg—
v2)
=
(e
—
VQ)I(Cg
—
e2).
(3)
This is an implicit equation defining p (or 7) in terms of v and e, which is just the type of equation needed for hydrodynamic calculations. 3. Behaviour in the single-phase region In the single-phase region the behaviour of the material may be described by an equation of state of GrUneisen type. This is of the form 0~ [e er(V)] Iv, (4) Pwhere = PrO’) ~r er are assumed or determined functions Pr +and along some reference curve. With some prescribed variation for ~r this allows analytical expansion away from this reference curve to cover the remainder of the —
(p,v) plane in question. vopo,,
I
/ /~‘~~t.d hq.d
4. Computer The computer subroutine subroutine EOSEXP calculates a valuespecific and for the internal pressure,energy, p, given e. It specific is clearvolume, from the v,
~
SPECIFIC VOWME,v
Fig.!. Subdivision of the pressure—volume plane for v >
vo.
2. Behaviour in the two-phase region In the two-phase region the specific volume v of the material may be written as V =
av~+ (1
—
(1)
a)v2,
and the specific internal energy as e=
~g + (1
—
a)eQ,
(2)
a is the mass fraction of vapour the mixture, v2 and Vg are the specific volumes and e2 and eg the specific internal energies of saturated liquid and gas at the vapour pressure p and temperature T, i.e., they are conditions on the saturation curve. Eliminating a from eqs. (1) and (2) gives .~
preceding sections that if conditions along the saturation curve can be determined then the behaviour of the material in the region v > v 0 can be completely described, for eq. (4) can be used with the saturation curve as reference curve in the single-phase region and eq. (3) can be solved in the two-phase region. Corresponding values of V1.j, Pri’ en, r’ri and Tn (1 ~ I ~ KTBL) are tabulated at a series of points along the saturation curve and input as data to the computer. For each data point on the liquid side of the saturation curve the corresponding point at the same temperature and pressure on the vapour side must also be input. The method used at AWRE for inputting this data on cards is illustrated in table l,which shows the saturation curve data for water. The first card specifies the number (KTBL) of data points on the saturation curve. There then follows KTBL cards, read in order of increasing ~‘ri’each one specifying vri, Pri’ en, Fri and Tn (1 e~i ~ KTBL). This data
66
K. Morgan, An expansion equation of state subroutine Table 1 The method of inputting saturation curve data on cards (shown for water) Specific volume 3/g) (cm
CARD 1. CARD 2 CARD 3 CARD 4 CARD 5 CARD 6 CARD 7 CARD 8 CARD 9 CART) to CARD 11 CARD 12 CARD CARD CARD CARD CARD CARD CARD CARD CARl)
13 j4 15 16 17 18 19 20 21.
Pressure (Mb)
Specific internal energy
Gamma
Temperature (°K)
(cm”Mb/g)
20 ~.,00OO0~.E00 t,043200~ 00 ~..U732ooE 00 1,127300E 00 3.,216300E 00 1,4~21ooE 00 1,611500E 00 1,769000E 00 2,039000E (JO ~,80O000E 00 i 1470000E 00 ~,8530ooE 00 8,3790~0E 00 1,1491.50E 01 1,802210E 01. 6,66163001 ~.,94~273E 02 O,060732E 02 t,694043E 03 1,29201O~ 05
1.300000E—08 0, 10,000000E—07 4,18O0O0~—03 3,000000EwO6 5,61.0000E03 i0,000000E06 7,620000E—o3 3.000000E!,05 1,004D00~.02 1o,oooooo~~05 1,393000E—o2 i.400000E04 1,~48O0QE.o2 1.,700000E—04 1,660000E—02 2.000000E—04 1,7B6000E..o2 2,209000E—04 1,959000E—02
i.981000E—01 4.410000E.Oi. 5.110000E01 5.760000EO1 ~~5~0Q0~E—0i 4.483O00~.o1 2.811Q00~—01 4.027000Em01 2.755000E’Ol 1.82300OE’~O1
2.209000E’04 2.000000E04 1.,700000E—o4 1,40UO00E~04 1o,000000E~05 3,000000E—05 1(l,000000EO6 3,O0000OE~06 lo.000000E—07 1.3OU000E’08
i.89li00O~N01 6,470000E 02 2.732000E”Ol 6,410000E 02 3.263~Q~E.0~.6,27OOO~E 02 2.829000E~01 6,110000E 02 2.530000E01 5,860000E 02 3,000000E’.Ol. 5,o80000E 02 3.000000B—01. 4,540000E 02 3.000000E01. 4,070000E 02 3.000000E01 3.730000E 02 3.000000E01. 2,803000E 02
2,0~8000E—02 2,294000Ee02 2,4o56ooE~o2 2,476700E—02 2.544600E-o2 2,602400E—o2 2,583900~—o2 2,545400E—a2 2.,.505500E.02 2,385600E—02
2.803000E 02 3,730000F 02 4,070000E 02 4,540000F 02 5.080000E 02 5,860000E 02 6.1i.0000E 02 6,270000E 02 6~41.o0ooE 02 6,470000E 02
Table 2 The method of storing the data of table! in array EOSTAB for use with subroutine EOSEXP N
1 2
3 4 5
6 7 S 9
10
11 12
13 14 15
16 17 18
19 20 21
22 23 24
25
EOSTAB(N)
N
EOSTAB(N)
N
EOSTAB(N)
2.000000E 01 1.~O~)O00E00 1.i143200E 00 1.~?3200E 00
3.oc0000E—05 i.000000E—04 1.400000E—04 1,700000E—04 2.o00000E—04
00
26 27 28 29 30 31 32
51 52 53 54 55 56 57
00
33
00
34
2.11~90OUE 00 2.8~’’)O0OE 00 3.4?0000E 00 5.~53000E 00 8.3?9000E 00 l,149150f~01 1.8fl2210E 01 ~.661630E 01 1.943273E 02
35
1.9SQ000E—02 76 2.u~8000E—02 77 2.254000E—02 78 2.405600E02 79 2.476700E—02 80 2.544600E—02 81. 2.oo2400E—02 82 2.5e~900E—02 83 2.545400E02 84 2,5C5500E—02 85 2.365600E02 86 1.9~1000F—O1 87 4.410000E01 88 5.il0000EOl 89 5,7e~OO00E—(11 90 5.5(i00U0F~01 91 4.483000E—01 92 2.811000E-01 93 4.fl27000E—01 94 2.7~5000E01 95 1.8231J00E01 96 1.8Y)(100E01 97 2.7~20(i0E01 98 3~263000E01 99 2.R29000E—0l 100 101.
i.~~7300E 1.216300E 1.452100E a.611SOOE 1.7~9000E
00
00
36 37 38 39 40 41 42 43 6.’60732E 02 44 i.6c404~3P 03 45 1.292010E 05 46 1.~0000E—08 47 i.nC0000E—06 48 3.flfJ0000E—06 49 1.i~t)000OE—05 50
2.2ij9000E—04 2.2c~.0o0E—04 2.n00000E—04
1.700000E—04 1.400000E04 1.000000E—04 3.000000E—05 1.000000E—05 3.000000E—06 1.o00000E—06 1.300000E—08 0. 4.:l80000E—03 5.610000E—03
7.620000E—03 1.oo4000E—02 1.3c3000E—o2 1.548000E—02 1.6e0000E—02 1.786000F—02
56
59 60 61. 62 63 64 65 66 67 68 69 70 71 72 73 74
75
N
EOSTAB(N) 2.530000E01 3.1100000E01 3.~t~0000E01
3.~~~000E01 3.)IJ0000E01 3~UO0O0E—01 2.,3u3000E 02 3.730000E 02 4.(37’lOOOE 02 4.540000E
02
5.OS0000E 02 5.8t~.0000E 02 6.110000E 02 6.270000E 02 6.41(I000E 02 h.470Ii00E 02 ~.47’1000E 02 6.410000E 02 6.270000fE 02 6.110000E 02 5.860000E 02 5.e0000E 02 4.540000E 02 4.070000~ 02 3.7~0O00E 02 2.803000E 02
K. Morgan, An expansion equation of state subroutine
67
must be fed into the array EOSTAB in the form required for use with the expansion equation of state
cific volume v by interpolating linearly along lines of constant pressure, i.e., for each data point!, a1 is cal-
subroutine. The first store, EOSTAB (IV), in the equation of state table for a particular material contains the number KTBL for that material. Stores EOSTAB (N + 1) to EOSTAB (N + 5KTBL) are then allocated as follows: EOSTABt’ri’ (N + 1) to EOSTAB (N + KTBL) contam the EOSTAB (N + KTBL + 1) to EOSTAB (N + 2KTBL)
culated, where = e + (e. e.) (v v .)/(v. v.) (7) / / I J £ j Here i = n + 1 j, and n is the number of data points input. This table is then searched to determine the interval (~,~_!‘ which e is with given (v,e) thenej)liesinbetween thelocated. isobarsThe (andpoint isotherms) / and!, i.e., p. ‘~p~p.andT. ~
containthepri, EOSTAB(N+2KTBL+l)toEOSTAB(N+ 3KTBL) contain the en, EOSTAB (N + 3KTBL + I) to EOSTAB (N + 4KTBL) contain the Fri, EOSTAB (N + 4KTBL + 1) to EOSTAB (N + 5KTBL) contain the Tn.
—
—
—
—
—
‘
j-1
/
(6) The quantities eQ, eg,vQ and vg used in eq. (3) are expressed in the form e (7) = e. + (e. e. —
/
~
j-!
eg(T) = e~+ 1+ (e~ e,~1)Z,
(8)
—
This process is illustrated in table2which shows how the data of table 1 is stored in array EOSTAB. Values of quantities between the tabulated points are obtained by interpolation. The method of interpolation used in the subroutine assumes that on the liquid side of the saturation curve Tr, er and Fr vary quadratically with ~r’and that Pr may be obtained from a relationship of the form 1 = A ‘T + B (5\ °~P~
1, r
“ ‘
!
On the vapour side of the saturation curve en and Tr are taken to vary linearly with Tr, and Vr and T~are assumed to be related by logv~A2/T~+B2.
(T) = V11 + (v1 v1 1)Z, v (T) = exp [log v. + A log v/v .÷ VQ
—
g
l
1
1
where Z = (T T,_ ~)/(~‘j 7~~) and A = (liT l/T1)/( lIT1 l/T1_ i)~ In this way eq. (3) may be regarded as an equation for Z and the solution for Z in the range 0 ~ Z ~ 1 is determined by the method described by Peters and Wilkinson [2]. A print of the iterations resulting from the application of this method can be obtained by setting the value 1.0 to the common variable EXIM. When Z (and hence T), has been determined, the pressure is obtained by using eq. (5). —
—
—
—
(6)
Eq. (5) is again applied to relate pressure and tempera-
5. Test runs
tune. The method of calculation in the subroutine may be split into sections (identifiable on the listing) as follows: (1) The tabulated specific volumes are searched to determine the interval in which v is located, (2) The interpolated quantitiesPr’ en, Fr and Tr are then determined as described above. (3) A check is made on the sign of(e er). (4) If(e en) > 0 then the material exists as a single-phase and p is determined by applying eq. (4). (5) If(e er) <0 then both phases co-exist, and eq. (3) is solved to determine p. The saturation curve data is used to set up a table of energy values at spe-
This expansion equation of state subroutine is used in the TOPAL hydrodynamic code which was developed at AWRE to investigate the problem of molten metal/liquid thermal interactions. The subroutine has performed satisfactorily for a variety of materials and Hoskin and Morgan [3] bave described results obtainèd using the subroutine for calculations on water and also on sodium. Special driver routines have been provided to enable test runs to be perfqrmed. Subroutine PRELUDE allocates storage for the 4ommon variables and the main routine TEST reads in the saturation curve data, a specific volume anda specific internal energy.
—
—
—
68
K. Morgan, An expansion equation of state subroutine
Subroutine EOSEXP then calculates the corresponding pressure and, depending on whether the calculated state point lies in the single or two-phase region, sets the value I or 2 to the common variable IMAR. The allocation of storage in this deck is peculiar to ATLAS 2. If the program is to be run on another machine then storage allocation should be specified in a manner suitable to that installation. This may entail modifying or removing PRELUDE and altering or replacing the COMMON and DIMENSION statements in the other routines*.
References Ii] E. Morris, private communication (1969). 121 G. Peters and J.H. Wilkinson, Computer J. 12 (1969) 398. 13] N.E. Hoskin and K. Morgan, Proc. Intern. Conf. Engineening of fast reactors for safe and reliable operation, Karlsruhe (1972), to be published.
TEST RUN OUTPUT
EXFA~SIO~EnS DATA FO~ WATER 20 E~Tt~IESFC!~ TARLE VCL(ME 1.0000uOE 10 1.1J43200E 00 1.07321OE 00 1.12?3uOE 00 1.2i.~3i0E (it) 1.4~21~’OE00 1.6i.15tDE 00 1.7690uOE 00 2.o3c000E 00 2.80001iOE 00 3.4701!(OE 00 5.8530t’OE JO 8.37couoE 00 1.14c1~0E ul 1.802210E 01, 6.6616~OE (Ii. 1.943273E 02 6.060732E 02 1.694043E 03 1.292010E ;5
V
2.78100Th 2.c2s000E 3.09000Th 3.274000’ 3.470000E 3.98100Th 6.31000Th 1.00000Th 1.00000Th 1.00000Th
01
Oil Or; 00
on 00 00 (11
02 03
OF’
CONDITIONS ON SATU8ATICN CLRVE
FI~ESSURE ).300000E—08 1(.000000E—07 ~-000000E—06 i0.000000E—06 ~.000QOOE—05 lU.000000E—115 1.400000E—fl4 1.700000E04 2.000000E—04 2.209000E—04 2.209000E—04 2.000000E—04 1.700000E—04 1.400000E—04 ~o-C00000E—05 ~.000000E05 1C.000000E—06 ~.O00OQ0E—06 io.IJ00000E—07 1.3o0000E—0~
E 2.075700E—02 2.071600E—02 2.067100E—02 2.062500E—02 2.(158000E—02 2.047300E—02 2.OoelOOE-02 1.964900E02 1.727600E02 1.5i~500E—02
ENFRC~Y 0. 4.18~jOQflE—03 5.hl0000F—03 20000E.03 7.6 1.OO4000E—02 1.393000E—02 1.548000F—02 1.660000F02 1.786000E—02 1.959000E—02 2.058000E—02 2.294000E—02 2.405600E—02 2.476700E—02 2.544i500E—02 2.602400E—02 2.583Q00E—02 2.545400E—02 2.505500E—02 2.385600E—02
P PHASE 2.998971E—04 1 2.799685E—04 I 2.600073E—04 I 2.4001.60E—04 1 2.209000E-’04 1 2.004294E—i)4 2 1.428224E—04 2 9.739631E—O~ 2 1.032795E—05 2 8.910063E—07 2
fiAMhA 1.981000E—Q1. 4.4j,U000E—tIl. 5.ll0000E—01. 5.76(’OOOE—Ol. ~.5.J’C0U0E—OI. 4.48300UE—01 2.811000E01 4.027000E01 2.755000E01 1.823000E—01 ~.89i:00UE~O1 2.732000E—01 3.263000E—O1 2.829OnuE—01 2.530000E—01 3.Ou0000E—01, ~.O0~000E—D1 ~.0flhb0tt0E.01 3.Ou0000E—01 30000’u’OE—01
T~ItPENATURE 2.Sfl3000E 02 3,73000Th
02
4.070000E 4,5400000202 ~.O~30000E 02 ~.8SO000E 02 6,11000Th 02 6,2?0U00E 02 6.41000Th 02 6.470000E 02 6.470000E 02 6.410000E 02 6.270000E 02 6.110000E 02 ~.860O00E 02 ~.080000E 02 ~,~4O00OE 02 4.07000Th 02 3,730000E 02 2.8)3000E 02