NBSIR 75-814

A MODIFIED BENEDICT WEBB-RUBIN EQUATION OF STATE

FOR PARAHYDROGEN II

H. M. Roder

R. D. McCarty

Cryogenics Division

Institute for Basic Standards

National Bureau of Standards

Boulder, Colorado 80302

Final Report

June 1975

Prepared for:

NASA-Johnson Space Center

Houston, Texas 77058

NBSIR 75-814

A MODIFIED BENEDICT-WEBB-RUBIN EQUATION OF STATE

FOR PARAHYDROGEN-II

H. M . Roder

R. D. McCarty

Cryogenics Division

Institute for Basic Standards

National Bureau of Standards

Boulder, Colorado 80302

Final Report

June 1975

Prepared for:

NASA-Johnson Space Center

Houston, Texas 77058

U.S. DEPARTMENT OF COMMERCE, Rogers CB. Morton, Secretary

NATIONAL BUREAU OF STANDARDS Richard W Roberts Director

CONTENTS

Page

1. Introduction 1

2. The Sources of Data 2

3. Adjustment and Modification of the Data 3

3. 1 Temperature Scale Changes

.

4

3. 2 Adjustment of Saturated Liquid and Vapor PVT Near the Critical

Point 5 3.3 Estimation of Critical Density and Temperature 6 3.4 The Revised Saturation Boundary Equations 8 3.5 The Revised Vapor Pressure Equation 9

4. The Equation of State for Hydrogen 12

5. The Properties Deck 14

6. Discussion, Intercomparisons, and Errors 18 6.1 Density Deviations along Isotherms „ 18

6. 2 Departures in the P-T Plane 23

6.3 The MBWR Extrapolated to High Densities 30

7. Summary 3 2

8. References 33

Appendix A. Program Listings 35

Appendix B. Test Program and Sample Results 50

Appendix C. Conversion Factors 50

TABLES

Table 1. Data Used to Determine the Parameters of the Equation of State . . 3 Table 2. Temperatures Assigned to the PVT Data of Michels, et al 4 Table 3. Least Squares Estimates of the Coefficients for Equation (1) . . . . 6 Table 4. Adjusted Saturation Densities Near the Critical Point 6 Table 5. Coefficients for Equation (4) 8 Table 6. Coefficients for Equation (5) Table 7. Saturation Densities and Deviations from Equation (5) 10 Table 8. Coefficients for Equation (6) 9 Table 9. Vapor Pressures and Deviations 11 Table 10. Coefficients for the Equation of State (7) 13 Table 11. Listing of Subroutines, Functions, and Entry Points of the

Properties Deck 16

iii
FIGURES

Page

Figure 1. The Lorentz-Lorenz Function for Hydrogen 5 Figure 2. The Rectlinear Diameter for Parahydrogen 7 Figure 3. Schematic Flow Chart for the> Propertie s Deck 15 Figure 4. Density Deviations along Isotherms 19

Figure 5. Density Deviations in the P-T Plane 24

Figure 6. Enthalpy Deviations in the P-T Plane 26

Figure 7. C Deviations in the P-T Plane 28

P

Figure 8, The MBWR at High Densities 31

SYMBOLS

A

B

N or G

nP

R

= parameters for equation 1

= parameters for vapor pressure

equation

= specific heat at constant pressure

= specific heat at constant volume

= parameters for equation for

saturation densities

= parameters for equation of state

= index of refraction

= pressure

gas constant

Subscripts

:

c = critical point

g = gaseous phase

= liquid phase

t = triple point

= parameters for index of

refraction equation

= temperature, T. the Inter- DO

national Practical Temperature

Scale

= reduced temperature

= scaling law parameter

= energy derivative

= specific heat input

= density

UNITS

The primary variables in the computer programs are Pressure in atmospheres

Density in moles/liter, and

Temperature in Kelvin.

Conversions to other SI units and units normally used in applied problems are given in appendix C.

iv
A MODIFIED BENEDICT-WEBB -RUBIN EQUATION OF

STATE FOR PARAHYDROGEN - II

Hams M, Roder and Robert D. McCarty

A 32 term modified Benedict-Webb-Rubin equation of state has been

applied to data for parahydrogen. The adjustable parameters in the equation

of state were determined using 2665 points including very recent measurements at low temperatures and high pressures. The new values extend the

range of the PVT data sufficiently to warrant a refitting of the equation of state. Temperatures for the data range from the triple point to about 700

K with pressures reaching 3000 atmospheres near ambient temperatures.

The PVT data were adjusted to the T^g scale. In addition, extensive

modifications have been made to the previously accepted PVT surface in the re gion near the critical point. These adjustments have been made on

the basis of more recent refractive index data and the application of scaling

law equations. Detailed comparisons between experimental and calculated

values are given for density. Corresponding comparisons are made for

enthalpy and the specific heat at constant pressure.

Key words: Critical point; density; enthalpy; equation of state; hydrogen;

index of refraction PVT; saturation properties; scaling laws;

specific heat.

1, Introduction

Almost all engineering problems requiring thermodynamic data for the cryogenic

fluids are most easily solved by using an equation of state to describe the PVT surface of

the gas, and to calculate values of such variables as enthalpy and specific heat. Quite often some modification of the Benedict-Webb-Rubin (1940) equation of state (hereafter referred

to as MBWR) is preferred, because equations of this type are relatively easy to handle on a computer.

This report describes an accurate wide-range MBWR equation of state for parahydrogen . The study of the equation of state was actually completed in two phases. The first

phase was sponsored by the NASA-Lewis Research Center (P.O. C-32369-C) and is summarized in a report by McCarty (1974). The earlier work gives a review of the adjustments made to achieve the International Practical Temperature Scale, T. , as defined in

o o Metrologia (1969). Changes were made in the PVT data, in the vapor pressure curve, in the

critical parameters, and in the two phase envelope near the critical point. Discussed briefly

is the selection of the equation of state, 32 terms rather than 19, a description of the fitting

program, and the results of a preliminary fit using the data available at that time. To

provide a complete record, much of the earlier material is included into the present report.

The second, and present phase of the study is sponsored by NASA-Johnson Space

Center (P.O. T-6570C). The main thrust is that new experimental PVT measurements at

* Parahydrogen is the nuclear spin modification which is stable at low temperatures.
low temperatures and high pressures are now available (Weber, 1975). These values extend

the range of the existing PVT data (Goodwin, et al. , 1963 and Michels, et al. , 1959a)

sufficiently to warrant a refitting of the equation of state. The refit now includes the actual

experimental PVT data, adjusted to the T, scale, and calculated values of C , adjusted bo v from normal to parahydrogen, for the high temperature source (Michels, et al. , 1959b).

The representation of the PVT surface and derived thermodynamic functions is greatly

improved over previous versions. This fact is shown clearly, and for the first time by

detailed comparisons between experimental and calculated values of density, and by

corresponding intercomparisons for enthalpy and specific heat at constant pressure, C . P

2. The Sources of Data

The major sources from which experimental or calculated values are taken originate

in the laboratories of NBS and of the University of Amsterdam. A description of these basic

references follows. Weber, et al. (1962) measured vapor pressures, Roder, et al. (1963)

present critical parameters and densities along the two phase envelope. The bulk of the

PVT data is given in Goodwin, et al. (1963), while experimental heat capacity data is found

in the papers of Younglove and Diller (1962a, b). These papers cover pressures from 0 to 340 atmospheres with temperatures from the triple point to 100 K, and they are smoothed

and combined to yield calculated thermodynamic functions by Roder, et al. (1965). The very

recent measurements by Weber (1975) extend coverage to pressures up to 800 atmospheres

with temperatures up to 300 K. The second extensive set of PVT data is presented by

Michels, et al. (1959a), and these values are used to calculate thermodynamic properties in Michels, et al. ( 1959b).

A number of other references on hydrogen exist, see for example the survey by

Woolley, et al. (1948). We have omitted these sets of data because they do not cover a large range of pressure and temperature, because the various sets are mutually inconsistent, and because the experimental errors are estimated to be larger than the sources

chosen. Indirectly these sources are included because they were used by McCarty and

Weber (1972) in deriving parameters for the 17 term equation of state described in that

paper. We use 40 points generated from that equation in the fitting of the present equation.

These generated points insure that the present equation is not subject to undue oscillations

for temperatures from 423 to 2200 K with pressures up to 680 atmospheres.

The input to the fitting program, the "data" is taken from these sources as shown in

table 1. It is clear from the table that in addition to using PVT data the fitting procedure

uses higher order thermodynamic data such as C and the Gibbs constraint, in other words

the technique of simultaneous or multiproperty fitting.

Normal hydrogen is the equilibrium mixture at room temperature, 75% ortho,- 25% para.
Table 1. Data Used to Determine the Parameters of the Equation of State

Type of data number of points source comments

PVT 1218 Goodwin, et al. (1963)

377 Weber (1975)

18 Roder, et al. (1965) extrapolated beyond the

melting line

482 Michels, et al. (1959a)

40 McCarty and Weber (1972) calculated, temperatures

above 4 23 K

38 Roder, et al. (1963) saturated liquid, saturated

vapor adjusted by McCarty

(1974)

C

V

163 Younglove and Diller (1962)

15 Roder, et al. (1965) calculated values

295 Michels, et al. (1959b) calculated values

ulDDS

constraint 19 Roder, et al. (1965) saturated liquid - saturated

vapor

Total Points: 2665

3. Adjustment and Modification of the Data

Adjustment of the raw data was necessary for two reasons. First, international

agreement on a practical low temperature scale which is very close to the thermodynamic

scale was reached in 1968 (T,_ Metrologia 1969). Second, an index of refraction experiment

6 8 by Diller (1968) indicated that the definition of the PVT surface given by Roder, et al. (1965)

is in error as much as 7% in density in the region near the critical point. The fact that the

two-phase envelope might be in error was pointed out as early as 1963, [see for example

figure 4 in Roder, et al. (1963)], through analysis of errors in the intersection temperatures

of isochores - experimental runs - and the vapor pressure curve.

As a result, McCarty (1974) made adjustments in both temperature scale and density

in the region near the critical point. The changes are most noticeable in the critical

parameters, the vapor pressure curve, and the saturated liquid and vapor densities. The

sequence adopted was to first find values for the densities of saturated liquid and vapor

from the index of refraction experiment; next to estimate new critical parameters from

these densities using values close to the critical point and a mathematical representation

based on the scaling laws. The results were checked by looking at the rectilinear

diameter. Finally, the newly defined saturation boundaries near the critical point and

densities at lower temperatures were represented with empirical equations which include

scaling law terms. PVT values generated from these equations were finally used as "data"

3
for the equation of state, see line 6, table 1,

3. 1 Temperature Scale Changes

The NBS 55 temperature scale was used to determine the PVT data of Goodwin,

et al. (1963), the heat capacity data of Younglove and Diller (1962a, b), and the PVT data

of Weber (1975). For these references the' conversion of the experimental temperatures to

the T . 0 scale was straightforward.

bo

Conversion of the temperature scale used by the University of Amsterdam is based

on a similar adjustment made for the PVT data of Argon (see Gosman, et al. , 1969). The

assumption is made that the same thermometer and scale was used for hydrogen as was

used for argon. The temperatures given by Michels, et al. (1959a, b) and the ones chosen

for the present fitting are contrasted in table 2. Two bits of evidence indicate that the

assignment of temperatures for this set of data are reasonable. First, the authors'

experimental temperature scale (Levelt-Senger s, 1966; included a calibration point at the

temperature of sublimating CO^. Thus the experimental temperature scale except for the

three lowest and for the very highest temperatures is remarkably close to the T. scale

6 o as subsequently defined. This is of course exactly what the authors were trying to achieve;

to make the measurements as nearly on the thermodynamic scale as possible. Second, a

Table 2. Temperatures Assigned to the PVT Data of Michels, et al.

Reference, °C T. „, K this paper

bo Reference, °C T, Q , K this paper

o o -175 98. 1835 - 25 248. 147

-170 103. 1835 0 273. 15 -160 113. 173 25 298. 142

-150 123. 1625 50 323. 140

-135 138. 1585 75 348. 143

-120 153. 161 100 373. 15 -100 173. 166 125 398. 1595

- 75 198. 165 150 423. 170

- 50 223. 1555

separate analysis also indicated that the temperatures shown in table 2 are the most likely

set. The analysis consisted of fitting a surface to this set of PVT data, first with the

temperatures adjusted from IPTS48 to T. , and then by assuming that the temperatures

6 8 were measured on a local scale (Levelt-Senger s, 1966) with subsequent corrections to IPTS48 and then to T . A comparison of the sum of the residuals of the two fittings favors

6 8 the second procedure.

4
3.2 Adjustment of Saturated Liquid and Vapor PVT

Near the Critical Point

In addition to the temperature scale change, an adjustment was made in densities on the basis of an index of refraction experiment by Diller (1968). Diller measured the

index of refraction, temperature, and pressure. He obtained densities from the PVT

surface defined by Roder, et al. (1965) which on the saturation boundaries is identical to Roder, et al. (1963). The salient graph from Diller's paper is reproduced in figure 1, where it is seen that the saturation boundaries fall on two separate legs which do not meet

at the critical density. The Lorentz-Lorenz function of the saturation boundaries should

have the same general shape as do the isotherms. In particular , for temperatures between

28 K and critical the two-phase envelope should be very close to the 35 K isotherm. The

two legs of two-phase envelope are ascribed to incorrect densities in the critical region.

More precisely, as indicated earlier, the errors stem from the intersection temperatures

of experimental runs and the vapor pressure cur r e, which near the critical point are nearly

co -linear.

.050

1.040

Eu

CM 1.030

.020

1 1 1 1 1 1

y- SATURATED VAPOR

1

0.5 %

^-298 K j

^35 K / ^50 K 25 K-i

<*¥f^ * f Urn »W 11 ! u 0-CK£^

\ \ 1 00 K-^ / /

\ \ A 35 K-7

\ 60 K-^ /\

40 KJ

pc j ^SATURATED LIQUID

1

! ll 1 1 !

1 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

DENSITY, g/cm 3

Figure 1. The Loreatz-Lorenz Function for Hydrogen

To adjust the saturation densities in the critical region the function

N

K=l

AK

n K-l (1)

was fit to the 35 K isotherm, where P is density in g/cm and n is the index of refraction.
Equation (1) was fit to 24 experimental index of refraction points with densities ranging

3 from 0.003969 to 0.078833 g/cm . It is important that the data not be "overfit. " Since

statistical significance of the coefficients to eq (1) was lost when more than 4 terms were

used the 4 term equation was chosen. The coefficients obtained for eq (1) are given in

table 3.

Table 3. Least Squares Estimates of the Coefficients for Equation (1) A = -1.0880215243

A

2

= 1. 8280271481

A

3 = -1.0378774469

A

4

= 0. 29788205862

The measured index of refraction along the saturated liquid and vapor was then used

to calculate the adjusted densities using eq (1). The densities given by Roder, et al. (1963)

and the adjusted densities are contrasted in table 4.

Table 4. Adjusted Saturation Densities Near the Critical Point

Terr P,K Density of Liquid, g/cm3 Density of Vapor ,g/cm3

T

Index of Roder et al. Index of Roder et al. NBS-55 68 Refraction (1963) Eq (1) Refraction (1963) Eq (1) 28 0 28 . 0071 1 092881 . 058966 . 058998 1 011312 .007298 007299

29 0 29 . 0073 1 089174 . 056646 . 056674 1 013767 .008887 008876

30. 0 30 . 0076 1 084824 . 053930 . 05 3 944 1 016892 .010882 010881

31 0 31 . 0080 1 079479 . 050589 . 050586 1 021040 .013541 013539

31 6 31 . 6082 1 075401 . 048057 . 048021

32 0 32 . 0084 1 072075 . 045993 . 045927 1 027273 . 017498 017525

32 4 32 .4086 1 067671 . 043320 . 043152 1 031 177 .019934 020017

32 7 32 . 7088 1 062766 . 040412 . 040057 1 035711 . 022664 022906

32 9 32 . 9089 1 055690 . 036941 . 035586 1 042271 . 025490 027074

3.3 Estimation of Critical Density and Temperature

Extrapolation of the rectilinear diameter yields the best estimate for the critical

density. A plot of (p , + P )/2 against temperature, as shown in figure 2, limits the critical

* g 3 density to a value between 0.03122 and 0.03142 g/cm . The value of the critical density is seen to depend only slightly on temperature. Therefore, the next step is to extrapolate the

rectilinear diameter numerically.

To achieve this the equation given in the next section, eq (5) which describes the

saturation boundaries is used. The equation is truncated to 4 terms resulting in

P g + P A

= 2P

c + (Glg + G u

) (AT)P + (G

2g + C u) (AT)

1 "* + (G

3§ + G^) (AT) , . ( 2)

6
0.033 E

u

+_ o.Oc

0.031

•

28 29 30 31

TEMPERATURE, K

32 33

Figure 2. The Rectilinear Diameter for Parahydrogen

where AT = (T - T)/T . Note that the exponents of the last two terms are 1-a and 1 c c rather than 1 and 4/3. G . is usually assumed to be equal to -G and the second te

1* J-g of eq (2) vanishes. There is some doubt about the validity of this assumption, but in this

particular case at least, the assumption cannot be disproved on the basis of the available

data. Thus the problem of estimating P is eliminated. A least squares fit of the eight

3 pairs of liquid -vapor data in table 4 resulted in a critical density of 0.03136 g/cm for

values of a = 0. 1 and T = 32. 933 K. The fit was repeated several times with values of

c

Qf and T ranging from a = 0. 1 to a = 0. 25 and T = 32.938 to T = 32.95. The

c c c resulting estimates of p did not vary significantly (i.e. , maximum variation was less

3 C than ± . 00001 g/cm ).

The next step is to estimate the critical temperature. In this case eq (5) is

truncated to two terms, that is to the terms originating from the scaling laws. The

resulting equation, eq (3),

p = p + G. (AT)

P

(3) sat c J

is obviously valid only for densities very close to critical. With some trial and error
temperatures between 31.9 and 32.8 were found to be applicable. A value for is at

hand, a value for T can be estimated from the equation if pairs of P and T are

C S cl t S el t available. Since very few of the experimental densities are at temperatures between 31.9

and 32.8 K, a parametric interpolation function for n was used to provide the necessaryinput to eq (1). The function

r

3 n = r

x + r 2 (AT) , (4) .

was fit to the index of refraction data of the saturation boundary (liquid and vapor) between

28 and 3 2.9 K. In eq (4) n is the index of refraction, and AT = (T - T)/T . The

coefficients used in eq (4) are given in table 5. By combining equations (4) and (1) densities

could be calculated

Table 5. Coefficients for equation (4)

r = 1.0509586594

T = 32.93313976

c

r

2

= 0.091463402563

r

3

= 0.41043983745

every 0. 1 K between 31.9 and 3 2.8 K. The data so interpolated were finally used in a fit of

3 eq (3) with p = 0.03136 g/cm , establishing the critical temperature T as 32.938 K,

c c and also 3 = 0.3483, P,. . = 0. 3478.

gas liquid

The values of the 3's are quite close to the 0.35 predicted by the scaling laws.

The value of the critical temperature is virtually the same as that calculated previously by

Weber from index of refraction data (unpublished but cited by Goodwin, 1970, page 226),

3.4 Revised Saturation Boundary Equations

The equation used to represent the saturation boundaries is based in part on the

scaling laws but is otherwise empirical. The equation is

,[1 + d-D/3]

"(1+1)

1=1

P Sat = P c + Gi (AT) 0 + ^ G (AT)

(g)

using the values established in the previous section for T , p and the betas, eq (5) is c c used to represent the saturation boundaries. Input for the fitting are the saturated liquid

and vapor values of Roder, et al. (1963) with temperatures adjusted to the T, scale,

o o except that for temperatures of 28 K and above the adjusted densities given in table 4 are

used, and that for the vapor a few generated values were added to balance numbers between

liquid and vapor. The coefficients for eq (5) are given in table 6. Table 7 gives the

deviations between calculated and input densities.
Table 6. Coefficients for Equation (5)

Vapor Liquid

P

c 0.03136 g/cm3 0.03136 g/cm3

B 0.3483 0.3479

C1 -0.047501571529 0.048645813003

G3 3.4871213005xl0~ 2- 3 . 4779278186x1

0

-2

G3 -4. 1221290925x10" 1 4. 0776538192x10" 1 G4 1.5666598550 -1.1719787304

GE -2.8061427339 1.62139244

Gs 2.7105455626 -1.1531096683

G7 -1.3074773595 0.3382549203 9 G8 0.22921285922 0.0

3.5 Revised Vapor Pressure Equation

The vapor pressure data of Weber, et al. ( 1962) was converted to the T scale and

68

refit to the nonanalytical vapor pressure equation of Goodwin (1969). That equation is

In (P/P ) = BX X + BS X

3 + B3 X

3 + B4 X (l-X) Bs

, (6)

where X = ( 1 -T^/T) /( 1 -T /T )

, T is in kelvins. The coefficients to eq (6) are given in

table 8. The two data points given by Weber, et al. (1962) for T = 22 and 23 K were

omitted from the fit because their inclusion seriously degraded the representation of the

Table 8. Coefficients for equation (6)

T = 13.8 K B

2

= 2.80810925813

T = 32.938 K B„ = -0.655461216567

c 3 P = 0.0695 atm B

4

= 1. 59514439374

B, = 3.05300134164 B^ = 1.5814454428 .1 5

rest of the data. The T triple point temperature of hydrogen, 13.81 K, could not be

6 8 used because it also degraded the fit. Similarly the boiling point of parahydrogen could not be constrained to the T value of 20. 280 K, but is rather 20. 277 K. Since Goodwin's

6 8 equation is sensitive to the values chosen at the triple point and has been sucessfully used

to force thermodynamic consistency for several other gases, we interpret the departure

at the triple point and at the normal boiling point to indicate some remaining low level

inconsistency in the defined temperature scale, T . Table 9 gives the vapor pressures

6 8 from Weber, et al. (1962), the adjusted temperatures, and the deviations between the

experimental and calculated data points. The critical pressure is obtained by inserting the

critical temperature into eq (6).

9
ftX"8

of)

COOmoOinfOO^NvOOH

NNNNNNNNNNNN OOOOOOOOOOOO

(M CO

in inoOCONNHiHNPO^^C^NmvOOvONOO^ONOO r0O00r^rs-OOOOOOO'-H'-<rH(MfMr0Tt< '^'v£)t^00 i-Hi-iOOOOOOOOOOOOOOOOOOOOO ooooooooooooooooooooooo

OOOOOOOOOOOOOOOOOOOOOOOOi-H^-if-f'-fOOOOOOOO

cr"

* E

cn 69 sQ

oooooooooooooooooooooooooooooooooooo

CO o CO i-H 00 o m in m -—

i

in (\J r- m f- o- tr oo ON 00 o CO Cr o 00 oo i-H

1

o- r00 o CO CO in in 00 X> x> N i oj co 00 CO a cr N m N CO cr CT on rO oo m rH oo O cr- CO

co m GO -v CO co CO CT O in in CO CO 00 (M og oo -J- X> cr N CO vO o 00 o Cr ro 00 ^ 1X1 IT) un in •o X> -0 o oj o r- co o oo r- O o o o o o o i-H i-H 00 00 CO in f- 00

00 (M (M oj oO oj 00 OJ N oo —

i

.—

i

o o o o o o o o o o o o o o o o o o o o o

o o O O o o O O o o O O o p o o o o o o o o o o o o o o o o o o o o o o

OOOOO-OOOOOOOvO^OvOcOi-HrtOr^Oi^i-HT^oOOOO^OOvDCOOOcOi-ii-iCO oooooooooooooooooor^r^i^r^oir^ooinr^coooo^i^oooooocor^r-r~r-r- m^fin^r^oooo^HOOOOOOOOOOOoooooor^ooooooooo oooooooooocooocr-ocr-r-^oooooooooooooooooooooooooo

a

XW "

4-) ^

Cq

NNNNNNN^lNN^l(^^^lHo^^oocxll^lp^lf>^o^oo(^oOH(^J(<l^ln so^co(^ (<lrnfni«imrOi«lMI<lcOl<im(n(<lMNNNHHHHHHH[MNNNNNNNNNN

Is- o J r- -f 00 00 vD cr m a vO r- 00 oj vO r- o o o 00 Co- o CO o in CO CO

00 in un oo oo 00 o- a Cr in O 00 00 00 cr CO in 00 r- o t*- 00 00 00 CO vO

in o cr- o in cr o o o r- in o o- CO cr in O 00 o 1—

1

oo i—

i

o r- a 0- o 00 o cr- m o co in CO o CO o 30 00 r- r- r~ vD in in r- -J3 in un ~t CO oo o Co- oo 0- vO oo oo -X) CO -f un in m m CO CO CO m co CO CO co CO r- r- r- o- r~ r- r- r- vO o O in in o o o o o o o o o a o o O o o o o o o o o o o o o o o o o o o o o o o o

c

0)u p CU ft

Tt<c000xti^HrHO00c0c000"-iOi-i00c0'^<'*—i^h^-i^hOOi—I'-Hi-l-HOO'-iOO'-iOcOTt1 oooooooooooooooooooooooooooooooooooo

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOC5I

a1 CO o in o oo ro 0- 0- o 0- CO vO o i—

i

in r- CO —

1

co- o un o

r- r~ S3 00 00 x> 'oo o vO o oj 00 r~- 00 00 vO cr- o r- cr- r~ a CO- -)< 00 ca CO ,—

1

co vO

lD o cr o in X CT- o 00 m CO o r- CO CN un o 00 co- i—

i

pH —

i

o r- o 0- o 00 o cr vO m o CO in X o CO vJ3 00 CO r- r- 0- 0- vO sO m un 0- vO in un CO i—

i

o CO- 00 Is- v£) 00 CO vO

-t -I- in m un in CO cO CO CO co ro CO CO CO r~ 0- r- r~ r- 0- 0- 0- 0- vO vO vO in in O o o o o o o O o o o o o o O o o o o o o o o o o o o o o o o o o o o o

* u

ft in

C>c»sOTtl OOOX)co.-HOOOOOOOOOOl^O'-H'-H^CX)OOC>oo^cocXlcOi-H^o^ ooooc»cx)oocoi^t^r^ooooooooocor^oounr^cocx)C^r^c»cooocoi^ ooooooooocoxfunot^ooo^ot-ioa^ooooooo-ooooooooo ^I^^ONOOOOOCXCXJCOCOCXJOOOOCT-O^OOCr-OOOOOOOOOOOOOOOOO

OOOOOOOOr-i^OC^OOOOOOOOOOOOOOOOOOOOcOCOunvO

COcOCOCOCOCOcOOOOOcOcOCOCOCOCOCOCOCO'—l^Hr-H^H 00000—looco^unvoi^ooo

1-i'-iCOOO NNNNNNNNN

10
Table 9. Vapor Pressures and Deviations

Pressure, atm

Zxd er ;~ er i&l T68' K

Pressure, atm

Zc (6) Percent

Diff.

1.6124 22. 0088 1. 6143 -0. 12*

2. 0688 23 . 0086 2. 0712 -0.11*

1 . 0000 20. 2770 1 . 0000 -0. 00

3. 2462 23.00^8 3. 2469 -0. 02

3. 9826 26. 0073 3 . 9822 0.01

4. 8285 27. 0071 4. 8275 0. 02

5. 7920 28 . 0071 1 n l c O. 7918 0 . 00

O O L O

: . 2 ; _ ; on n n 1 / OQA "1 0 . oo4

f

n no U . UZ

8. 1162 30 . 0076 0 i i A n 8 . 116 9

n n l - U . U 1 8. 1169 o n n n "7 A 3U . UU to O 1 l(ln 0 . 1 io V

n nn - U . UU

8. 1171 on r\ n "7 £. -30. UU f O 8. 1 io9 n nn U . UU

8. 7873 30 . 5078 0 Tom 8 . 78 91 n n 0 - U . UZ

8. 7885 30. o078 8 . 78 91 n ni - U . U 1 8. 7886 30 . 5078 8.0 70m 7891 n a 1 - U . U

1

9. 5029 3 1 . 0080 9. 5010 n no U . UZ

9. 5023 *j t r\ ^ 0 r\ 31 . 0080 9. 5010 n

U . Um1 9. 5005 .si . 0080 9. 5010 n n l -0.01

9. 5003 31 . 0080 9 . 5010 n ni - U. Ul

10.2525 3 1 . 5082 1 0 . 2546 n no - U . UZ

10. 2535 3 1 . 5082 in a L. 1 U . Zo4o A A 1 - U . U 1 10. 2539 3 1 . 5082 1 f\ O Z. A A

1 U . Zo4o n n i - U . Ul

1 1 . 0502 o£ • UUo-± 11. UoZo — AU . A9U£

1 1 . 0516 O O AAO /I 32 . 0084 11 n:7Q 1 i . U0Z8 A A 1 - U . U 1 11. 0522 32. 0084 ii n;io A A 1 - u . u 1 11. 8988 32 . 5087 it Q n 0*5 11. 8 99Z - A A A U . U U 11. 8976 32. 5087 ii 0 n no

1 1 , 8 99Z A A 1 - U . U

1

1 1 . 8989 32 . 5087 11 2 Q Q 0

1 1 . £5 7 7Z A A A 12. 0749 :>2. o087 1 Z . O / 48 AU . AUUA 12. 0742 32. 6087 i 0 n 7 .4 Q 1 Z . U ( 40 - AU . AAU U 12. 0751 32. 6087 l 0 m< 0

1 Z . U (4o A A A U . U U 12. 2526 32. 7088 1 0 0 c O — 1Z. Z3Z f - A A A U . UU

12. 2520 32. 7088 12. 2527 - A A 1 U . U 1 12. 2536 32. 7088 12. 2527 A A 1 u . u 1 12.4326 32. 8089 12. 433U - A A A U . UU

12. 4330 32. 8089 ii y)i?n

1 Z . -i5 5\J AU . AUUA 12. 4352 32. 8089 1 0 /I ^ "2 A 1 Z . 7jj U A A?

12. 6168 32. 9089 19 A 1 A A 1 Z . O 1 0 U u. Ul

12. 6187 32. 9089 1 0 A 1 A A 1 Z . 0 1 O J A A?

12. 6183 32. 9089 1 1 A 1 A A li, O 1 0 U

0. 0778 1*5 O n "7 *"7 13 . 9977 A H77fl U . U / 10 0. 1327 15. 0020 n 1 "3. 0 *7 U . 1 -5Z

1

n m

0. 2129 16. 0051 n 717Q U . Z 1 Z V 0. 3250 17. 0071 n i7cn U

.

jZdU AU . ATU 1 0.4759 18.0084 0.4759 -0. 00

0.6726 19.0088 0.6727 -0. 00

0. 9228 20.0090 0. 9229 -0. 02

0.0695 13. 8000 0. 0695 -0. 00

32. 9380 12.6698**

* Points omitted from the least squares fit. ** Critical pressure

11
4. The Equation of State for Hydrogen

Since the major modification of the Benedict-Webb -Rugin (1940) equation of state by

Strobridge (1962), there have been many more. Each author claims his particular

modification to be the best of several he has tried for the particular fluid being correlated.

In some cases a given form was chosen because it worked well for a number of fluids.

Several of the MBWR's have been applied to hydrogen. Strobridge's equation (16 terms) was applied by Roder and Goodwin (1961) to parahydrogen. It was found that two sets of

coefficients, one for liquid and one for gas, were required to reproduce the experimental

PVT surface. In 1967 (see Roder, et al. , 1972) a 17 term equation was applied to values

above 50 K, thus omitting the two-phase region entirely. In this fit major discrepancies

remained at the junction of the two sets of experimental data near 100 K, and in addition the

enthalpies around ambient temperatures showed statistically significant deviations. A

subsequent refit of this equation by McCarty and Weber (1972) included values of C in the

data set. While the enthalpies near 300 K were improved the departures in PVT at 100 K

remained substantial (see figure 4h this report, the line labelled TN 617).

In phase I of this work McCarty (1974) studied both the 19 term version by

Bender (1970) and the 3 2 term version by Jacobsen (1972). He selected the 3 2 term equation

as being superior, and that equation of state is used here. Actually, if the term PRT is counted there are 33 terms, and if the coefficient of the exponential term y is counted there are 33 parameters. The equation of state is: P = pRT + p

2 (N: T + Njj T

1/2 + N3 + N4 /T + N5 /T2

) + p

3 (N6 T + N7 + Ns /T + N9 /T2

) + p

4 (N10 T + Nlt + N18 /T) + PB (Nia ) + P

S (Nl4 /T + N1B /T2

) + p

7 (Nie /T) + p

8 (N17 /T + N18 /T2

) + p

9 (N19 /T2

) + P

3 (N30 /Ta + Nai /T3

) exp (-YP

2

) (7) + P

5 (N33 /T2 + N33 /T4

) exp (- Y p

3

) + P

7 (N34 /T2 + N3B /T3

) exp (-yp

2

) + p

9 (N3S /T2 + N37 /T4

) exp (-yp

2

) + P 11 (N38 /T2 + N39 /T3

) exp (-yp

2

) + p 13 (N30 /T2 + N31 /T3 + N33 /T4

) exp (- Y p

2

) The coefficients for the equation of state were estimated from a least squares fit using the

data set indicated in table 1 with p in moles /liter, P in atmospheres, and T in kelvins.

For these units the values of R and y are given below.

12
Table 10. Coefficients for the Equation of State (7) R = 0.082056 16 4 . atm/mol. K. 0.0041

G( 1) = 4.61438775565437326033D-04

G( 2) = 4.233184556086770434400-02

G( 3) =-5.096556226503733321570-01

G( 4) = 2.923059738269586053460*00

G( 5) =-2.?87609i472il36g29049D*01

G( 6) = 1.38314860141Q70378866D-05

G( 7) =-1.322256954639226520670-03

G( 8) = 3. 016504431701892492910-01

G( 9) = 5.Q937Q556C851742825920+01

G(10) = 1.97382832491904714077D-07

G(ll) = 2. 858492039828227170630-04

G(12) =-2. 228 279239123480 57045D-02

G(13) =-2.2 57481136 7643 040 697 20-06

G(14> = 2.41427236974667590421D-05

G(15) =-1.69571339 8588410470130-03

G(16) =-5.393 6763912751931915lD-07

6(17) = 3. 998955244328083808620-09

G(18) = 1. 142457561274493541050-06

G(19) =-1.2 52 5662 25 8960 5271*1230-08

G(20) =-<. 917 8619348826398 82960 + 01

G(21) =-1.585666017368677796970*02

G(22) =-1.901602946272185543660-01

G(23) = 9.19802086250050278199D+00

G(24) =-3. 180455 51 88 10 44498 74 ID- 04

G(25) = 1.19105779192652709183D-03

G(26) =-3.791352773225991761320-07

G(27> =-3.983377599095395450920-05

G(28) =-1.234510854688972907080-10

G(29) = 1.950266293499069896810-09

G(30) =-2.380343917109169846870-13

G(31) =-4. 0 73 5766 0 8 192 893 866180-13

G(32) = 3.801354930777624867160-12

Several comments on the final fit of the equation of state are appropriate. The critical

point was constrained to the value P = 12.670 atm, p = 15. 556 moles per liter, 2 2 T = 32.938 K and the derivatives to dP/dp = 9 P/dp = 0. The thermodynamic conditions

for phase equilibrium for the coexisting liquid and vapor phases have been included as data

in the least squares estimating procedure. The fitting of the equation of state was attempted

with several combinations of the available data. Representation of the data improved in steps as the data set was changed. Initially we added Weber's (1975) new measurements.

Next, we replaced all but 40 points of the data generated from the PVT surface of McCarty

and Weber (1972) with the experimental PVT data of Michels, et al. (1959a). A significant

improvement occurred when we added the calculated values of by Michels, et al. (1959b).

A final improvement resulted when we added 15 calculated values of along the saturated

vapor line, and 18 generated PVT points near the melting curve. In both locations the

definition of C was improved considerably.

13
5. The Properties Deck

The "properties deck" is a collection of subroutines and functions designed to return

a wide variety of state variables, thermodynamic properties and derivatives. A listing of

the deck is given in appendix A, a test program and sample results in appendix B. For the user we classify the programs into initialization, basic programs, and second level

programs.

Initialization . The first step in a using program would be to call the data subroutine

DATAPH2. This routine is normally called only once to set the coefficients of the equation

of state, the vapor pressure curve, etc. It will be evident from the program listings given

in appendix A that an index N can be associated with this routine and certain of the ideal

gas functions. N identifies the species of hydrogen under consideration, and for the

purposes of this report is always assumed to be = 1, i.e., we are considering only

parahydrogen.

Basic programs . In most problems two of the variables in the set of pressuredensity-temperature are known, and the requirement is to find the third variable. Accordingly the three possibilities are

1) density and temperature known, find pressure -* subroutine PRESS(PR, D, T),

2) pressure and temperature known, find density -' function FINDD(P,T), and

3) pressure and density known, find temperature -» function FINDT(P, D)

Since the equation of state is explicit in pressure the subroutine PRESS is straightforward.

However, both FINDD and FINDT have to be based on an iterative solution of the equation of

state using some initial guess and appropriate derivatives of the PVT surface.

Second level programs . The common denominator for the second level programs is

that they require as input some combination of pressure, density, and temperature. More

often than not a basic program has to be called before a second level program can be used.

All of the remaining subroutines and functions are included in this grouping, that is all phase boundaries, all derivatives and integrals of the equation of state, all thermodynamic

functions, and all other properties such as transport properties or dielectric constant.

Flow chart and table of programs . The schematic flow chart, figure 3, presents

the logic a user has to apply to any given problem. Approximately 30 programs combining

about 50 possible entry points comprise the spagetti bowl which is loosely labelled

"properties deck." Each of the entry points corresponds to a single result, answer, output,

property returned. These 50 or so possible end points are shown in table 11, which is arranged according to the input these programs require. It is evident from table 11 that

density is the input most often required.

One peculiarity in the structure of the properties deck should be noted. The sub-
INPUT

zez

Call

2nd Level

Progra ms

Subroutines

Required

Result

Property Desired

OUTPUT

No

Initialization Required?

Yes

Index N,

Entries are optional

Indicate H2 Species

Determine

Properties Required

Determine

INPUT Needed

Note:

N=l, i.e. Parahydrogen

for this report

END

Pressure Density Temperature

PRESS FIND D FIND T

Figure 3. Schematic Flow Chart for Properties Deck

15
0uQ

;ed

Calcu]

RESULT, Property

ft>HQ

ft « Q «w - w 2 W

ft 05 > ft

ft &.

ft ft

HH

73

ft>u3O

Ui3n

<n

<uM

ft0

a.>U

U J! M3 M-l

OQ 0

01

OJ o

pr

ti\

<u

or *t a u

a) 1) > v

ft a

J3

ft °

73

si

V

-Cu

ft

cif

a 0h

V

On

"3 s m

yi 9) 0] n) d M 00 bo

•a

a

4> 732 2 « « a H H H < < <

33O

ftu

H

X ~ r P - 2 _ ft 2 D

ft ft > Q ><w ft

h

r

—

ft ft Z

a) <D

t. 3 3

rth

a)ftSH

H

«M .2 <M

waw

CO r,

U *> ft

i-i cy -i- +j !-< aw

•W +J 4J 4-> 4^

ff ° - 05 Q. s

OS

w

ft

w ~ U W 05 S S:

05

ft

16
17
routine PROPS is a multiple entry routine designed to reduce ttie number of operations

involving the coefficients and terms of the equation of state. Included in closed form are

derivatives and integrals of the equation of state which are required to obtain enthalpy and

entropy.

6. Discussion, Intercomparisons and Errors

Hydrogen may be the first, and perhaps the only case where we can determine, with

reasonable assurance, what the errors in the PVT surface and in the derived properties

actually are. This fortunate circumstance arises because for hydrogen there is a wealth

of data available which can be used both as input to the equation of state and to check the

quality of the MBWR, We show the results of extensive comparisons in two different ways.

The first set of graphs might be called "standard, " because it is the conventional way of

plotting density deviations for a set of selected isotherms. The second set of graphs gives an overview of density, enthalpy, and C errors in the P-T plane. In the last part of this

P

section we look at the MBWR extrapolation to high densities.

6. 1 Density Deviations Along Isotherms

Density deviations are plotted in 16 segments of figure 4 for isotherms of 26, 33,

60, 98/100, 150/153, 198/200, 298 /300 and 423 K. The deviations, expressed in percent, are the differences between values predicted by the MBWR on the one hand and the PVT data

of Goodwin, et al. (1963), Weber (1975), and Michels, et al. (1959a) and one prior correlation by McCarty and Weber (NBS Technical Note 617, 1972) on the other hand. In these

graphs the 32 term MBWR is the zero or reference line, and the departures for each

isotherm are plotted against both pressure and density. One of the reasons for plotting

against both variables is evident for the 26 K isotherm, figure 4a and 4b. The plot against

pressure is continuous, the plot against density shows both vapor and liquid segments of

the isotherm. Figure 4 illustrates that the differences between the data of Goodwin and the

correlation of TN 617 — a polynomial smoothing of the surface — are negligible. The plot

against pressure in figure 4c is typical of temperatures near critical. This plot shows the

inability of an analytic equation of state, the MBWR, to accurately represent a PVT surface near the critical point. In the corresponding plot 4d it is seen that these departures extend over quite a large range in density, and also that there is some difference between the

experimental data (Goodwin, et al. , 1963) and the polynomial smoothing (Roder, et al. , 1965). The lowest temperature at which we are able to intercompare all sources is 98/100 K, figure 4h. We note that the rather sharp change between Goodwin, et al. (1963)

and Weber (1975) first seen at 60 K is also evident at 100 K. At low densities Goodwin and

Michels disagree, at higher densities the agreement between Weber and Michels is very

satisfactory. We note the large deviation of TN 617. We are forced to conclude that the 17

18
CO

r- o

• Z

t-< o

it

c 0

rj o 0 °

0) w

0c

-J 5- W

£

Li0

f\] CO o 0 m

_ >

All SN3Q MI lN33d3d All SN3Q MI !N30d3d

-> r

or =»

to O

&0o

>-

cQ

<03

CO

A1ISN3Q NI iN33d3d A1ISN3Q NI !N33S3d

19

A1ISN3Q NI iN33H3d All SN3Q NI !N30a3d

21
A1ISN3Q NI !N3D83d AiISN30 NI lN3083d

T3

ac

oU

Al I SN3Q NI lN33«3d A1ISN3G NI lN33U3d

22
term MBWR used in TN 617 does not offer sufficient flexibility to represent an entire PVT

surface

.

Of the remaining plots figure 4n is particularly interesting. The general progression

of Michels 1 and "Weber's data is quite similar up to a density of 30 mol/1. Can the difference

be ascribed to different assumptions about the stretching of the PVT pipet? The drastic

change in Michels' data at densities above 30 mol/1 is not understood, particularly in view

of figure 4p. Could the change be a one-time nonelastic stretching of the pipet? The

behavior of the 17 term MBWR (TN 617) is also of interest in figure 4n. At low densities

this surface follows Michels' data, at higher densities it switches to represent the extrapolation isochores of Goodwin, et al. (1963). The latter can now be seen to be qualitatively

similar to the new measurements of Weber (1975).

6.2 Departures in the P-T Plane

These deviation plots involve density because density is the most basic variable

(see table 11), and enthalpy as well as C^ because these are the most important properties

for the engineer. The graphs are assembled in three sets of four plots each. The three

sets should be considered together since they are interdependent. Each set is obtained by

plotting the P-T locus of successively larger errors of the variable under consideration.

From these sets it is quite clear that the region near the critical point is the only area which is really troublesome.

The error plots for density, figure 5, are in percent without regard to sign, and

are taken directly against the two major sources of experimental data, NBS and Michels.

The average deviation in density is 0. 14%. This is slightly larger than the average

experimental uncertainty, and is about the best one could possibly expect.

The error plots for enthalpy, figure 6, are in J/mol for the simple reason that the

enthalpy values go through zero, from about -600 to 15000 J/mol. The comparison is between different methods of calculation, values are compared at the experimental points

of PVT. Values calculated with the 32-term equation are compared with those calculated

byRoder, et al. (1965), by Weber ( 1975) and by Michels, et al. (1959b). The graphs

illustrate clearly that given a density deviation there will be a corresponding departure in

enthalpy. At first sight departures of up to 30 J/mol at the higher temperatures, 300 - 423 K, are startling. A little reflection shows these errors to be almost negligible. We

recall that what is plotted are differences in enthalpy, not percent. A specific example

concerns the ideal gas value at 300 K which is about 8400 J/mol and the enthalpy value at 2200 atm and 300 K which is 11600 J/mol. Given these values a difference of 30 J/mol for

the isothermal increment of about 3200 J/mol is not severe. From figure 6 we estimate

the average enthalpy deviation to be around 3 J/mol.

23
100 TEMPERATURE, K

A, DEV. GREATER THAN 0. ! PERCENT

Figure 5,

in Density Deviations in the

P-T Plane a~

4'0 SO 80 ISO 200 TEMPERATURE, K

B, DEV. GREATER THAN 0.3 PERCENT

24
so eo too 20o TEMPERATURE, K

D, DEV. GREATER THAN 1.0 PERCENT

Figure 5. Continued

100 TEMPERATURE, K

-0, DEV. GREATER THAN 0.5 PERCENT

25
i

2

*—'—+-

SI 80 !05 2! TEMPERATURE, K

A, DEV. GREATER THAN L 0 J/MOL

-f

Figure 6. Enthalpy Deviations in th<£02." P-T Plane

S: 63 lie 2!! TEMPERATURE, K

B, DEV. GREATER THAN 3. 0 J/MOL

26
I'll

t

• • . I

;

•

;

:

2::

CO

co 5 -

.

a:Q_

A ' s: f:

1 ±::: 2: TEMPERATURE, K

D, DEV. GREATER THAN v'/MOL

1 :—' . t

Figure 6. Continued

'Empera":re, k

27
I ' I —

2::

£25-

a:

i::-

i: s: 8": iflo 21: TEMPERATURE, K

A, DEV. GREATER Th'Ar: :. C PERCE':

2::

Figure 7. q: Cp Deviations in the P-Tuj2 - Plane

a:

-t——1—1—»-

1 « so

' foo ' zT: TEMPERATURE, K

B, DEV. GREATER THAN 3. 0 PERCE? J

2 p.
I—I—I—I—t411 •

2::-

U1

en

7^ 80 100 200 TEMPERATURE, K

D, DEV. GREATER THAN !0. PERCENT

Figure 7. Continued

t: s: bo :oo 2: TEMPERATURE, K C DEV. GREATER THAN' 5. C PERCE'

2Q
The C plots figure 7, which are in percent, show that the departures at the higher

P temperatures are negligible, just as the enthalpy departures would be had we expressed them

in percent. The effect of adding data along the saturated vapor line (line 8, table 1) is

to reduce departures in from several percent to less than 1 percent. Similarly the

addition of extrapolated PVT data (line 3, table 1) improves the errors in C near the

P

melting line from greater than 10 percent to 3 percent or less. In addition, an inconsistency

in C noted in the earlier correlation (see isotherm 140°R, figure 5) of McCarty and

Weber (1972) is removed. From figure 7 we estimate the average departure in to be

somewhat less than 3 percent.

6.3 The MBWR Extrapolated to High Densities

One of the peculiarities of the MBWR is that occasionally the iteration to find density

yields an invalid result*. The situation can be understood by considering figure 8 where an isotherm of the MBWR is extrapolated to high densities. The range of valid PVT data is normally somewhere to the left of the maximum in pressure. It is easy to visualize an

iteration using the slope (3P/9p)^, of the surface and an initial density of p = 0 yielding a density to the right of the maximum in pressure, i.e. , an invalid result. A case in point

is the addition of extrapolated PVT data (line 3, table 1) mentioned in section 6. 2. Adding

these 15 generated points shifts the maximum in pressure for the 23 K isotherm from ~ 380 atm to ~ 1000 atm. It is the change in slope (BP/BP ) near the melting line ~ 200 to 350 atm that results in the desired change in C .

The present equation of state yields only negative pressures for densities beyond

the cutoff shown in figure 8. We note that this behavior is diametrically opposed to what

one might anticipate. In general one would expect pressure to increase as the density

increases unless a phase transition is encountered. At the possible phase transitions

liquid /solid and molecular /atomic the pressure would remain constant for the associated

change in density. Several attempts to force the 32 term MBWR to return positive

pressures at very high densities have not been sucessful.

Occasionally described as "the equation of state blows up."

30
ure 8. The MBWR at High Densities

31
7 . Su mmar y To obtain the desired equation of state 2665 points were used including very recent

measurements at low temperatures and high pressures. All PVT and C

v data had to be

changed to the T , _ scale, a non-analytical equation was applied to the vapor pressure, and

6 8 the results of an index of refraction experiment were applied to the critical parameters and

the two phase envelope near the critical point. Multi- proper ty fitting of PVT and data and

the imposition of Gibbs phase constraints insured a substantially improved representation of

the various thermodynamic quantities.

The resulting equation of state is valid for temperatures from the triple point, 13.8 K,

to the onset of dissociation ~ 1500 K; it includes pressures gradually increasing from 700

atmospheres at the melting line to 3000 atmospheres at room temperature. In practical

terms pressures up to 12000 psia are included, and the equation can be used for temperatures

up to 5000°R if appropriate arrangement is made for dissociation.

Two characteristics of this equation of state must be kept in mind. First, the equation

is analytic in nature, this means the critical region cannot be represented accurately.

Second, the limiting behavior at high densities does not correspond to our a priori expectations, this means care has to be exercised to stay within the valid range of the equation.

The equation developed here sets a standard of what can be achieved in the fitting of an equation of this type. The high quality of the surface, illustrated in extensive deviation

plots, is possible because ample data over a wide range of temperatures and pressures with

inherently high precision is available as input. The PVT surface defined by the equation

needs no further numerical treatment unless new experimental data become available, or the

international temperature scale is redefined.

32
8. References

Bender, E. , Equations' of State Exactly Representing the Phase Behavior of Pure

SuOstances, (Proc. 5th Symp. of Thermophysical Properties) 227-35, ASME,

N. Y. (1970).

Benedict, M. , Webb, G. B. , and Rubin, L. C., An Empirical Equation for Thermodynamic

Properties of Light Hydrocarbons and their Mixtures I. Methane, Ethane, Propane,

and n-Butane, J. Chem. Phys. 8_, 334-45 (1940).

Diller, D. E. , Refractive Index of Gaseous and Liquid Hydrogen, J. Chem. Phys. 49,

No. 7, 3096-3105 (Oct 1968).

Goodwin, R. D. , Nonanalytic Vapor Pressure Equation with Data for Nitrogen and Oxygen,

J. Res. Nat. Bur. Stand. (U.S.) 73A (Phys. and Chem.), 487 (1969).

Goodwin, R. D. , Estimation of Critical Constants T , P from, the P (T) and T(P) Relations

c c at Coexistence, J. Res. Nat. Bur. Stand. (U.S. ) 74A , No. 2 (Phys. and Chem.),

221 (1970).

Goodwin, R. D. , Diller, D. E., Roder, H. M. , and Weber, L. A., Pressure-DensityTemperature relations of Fluid Para-hydrogen from 15 to 100 K at Pressures to 350 atmospheres, J. Res. Nat. Bur. Stand. (U.S.) 67A (Phys. and Chem.), No. 2, 173-192 (Mar-Apr 1963).

Gosman, A. L. , McCarty, R. D. , and Hust, J. G. , Thermodynamic Properties of Argon

From the Triple Point to 300 K at Pressures to 1000 Atmospheres, Nat. Bur. Stand.

(U.S.), NSRDS-NBS 27, (Mar 1969).

Jacobsen, R. T. , The Thermodynamic Properties of Nitrogen from 65 K to 2000 K with

Pressures to 10,000 atm, Ph.D. Thesis, Washington State University (Jun 1972).

Levelt-Sengers, J. M. H. , private communication to J. G. Hust (1966).

McCarty, R. D. , A Modified Benedict-Webb-Rubin Equation of State for Par ahydrogen,

Nat. Bur. Stand. (U.S.), NBSIR 74-357, (Feb 1974). Available National Technical

Information Service, Springfield, Va. , Order No. COM-74-10551 - Price: Paper

Copy $6; 75 - Microfiche $2.25.

McCarty, R. D. , and Weber, L. A., Thermophysical Propertie s of Parahydrogen from

the Freezing Liquid Line to 5000 R for Pressures to 10,000 psia, Nat. Bur. Stand.

(U.S.), Tech. Note 617 (Apr 1972).

Metrologia, The International practical Temperature Scale of 1968, Metrologia 5, No. 2, 35 (Apr 1969).

Michels, A., de Graaff, W. , Wassenaar, T., Levelt, J. M. H. , and Louwerse, P.,

Compressibility Isotherms of Hydrogen and Deuterium at Temperatures between

-175 = C and +150°C (at densities up to 9b0 Amagat), Physica 25, 25-42 (1959a).

33
Michels, A., de Graaff, W. , and Wolkers, G. J., Thermodynamic Properties of Hydrogen

and Deuterium at Temperatures between -175°C and 150°C and at Densities up to 840 Amagat, Physica _25, 1097- 1 1 24 ( 1959b).

Roder, H. M. , Diller, D. E., Weber, L. A., and Goodwin, R. D. , The Orthobaric

Densities of Parahydrogen, Derived Heats of Vaporization, and Critical Constants,

Cryogenics 3, No. 1, 16-22 (Mar 1963).

Roder, H. M. , and Goodwin, R. D. , Provisional Thermodynamic Functions for Parahydrogt

Nat. Bur. Stand. (U.S.), Tech. Note 130 (Dec 1961).

Roder, H. M. , McCarty, R. D. , and Hall, W. J., Computer Programs for Thermodynamic

and Transport Properties of Hydrogen (Tabcode-II), Nat. Bur. Stand. (U.S.), Tech.

Note 625 (Oct 1972).

Px.oder, H. M. , Weber, L. A., and Goodwin, R. D. , Thermodynamic and Related

Properties of Parahydrogen from the Triple Point to 100 K at Pressures to 340

Atmospheres, Nat. Bur. Stand. (U.S.), Monogr, 94, 112 pages (Aug 1965).

Strobridge, T. R., The Thermodynamic Properties of Nitrogen from 64 to 300 K between

0.1 and 200 Atmospheres, Nat. Bur. Stand. (U.S.), Tech. Note 129 (Jan 1962).

Weber, L. A., Thermodynamic and Related Properties of Parahydrogen from the Triple

Point to 300 K at Pressures to 1000 Bar, Nat. Bur. Stand. (U.S.), NBSIR 74-374

(Mar 1975), also NASA Spec. Publ. No. SP-3088, 105 pages (Mar 1975). Available

National Technical Information Service, Springfield, Va. , Order No. NASA

SP-3088 - Price: Paper Copy $5. 25.

Weber, L. A., Diller, D. E., Roder, H. M. , and Goodwin, R. D. , The Vapour Pressure

of 20 K Equilibrium Hydrogen, Cryogenics^, No. 4, 236-238 (Jun 1962).

Woolley, H. W. , Scott, R. B. , and Brickwedde, F. G. , Compilation of Thermal Properties

of Hydrogen in Its Various Isotopic and Ortho-Para Modifications, J. Res. Nat. Bur.

Stand. (U.S.) 41, 379-475 (Nov 1948).

Younglove, B. A., and Diller, D. E., The Specific Heat of Saturated Liquid Parahydrogen from 15 to 32 K, Cryogenics 2, No. 5, 283-287 (Sep 1962a).

Younglove, B. A., and Diller, D. E. , The Specific Heat at Constant Volume of Parahydrogen at Temperatures from 15 to 90 K and Pressures to 340 Atmospheres,

Cryogenics 2, No. 6, 348-352 (Dec 1962b)

34
Appendix A. Program Listings

SUBROUTINE DaTa P HZ DIMENSION 5(32), vP(8),GI(ll)

01 HENS ION GV (9) ,GT(9) , FV (C> ,FT»4J ,EVI8I , ET(4) COMHOn/CPID/GI

CONHQM/CRIT/ EM, EOK, KM, TC, DC, X , PC, SIG COHHQN/DAT A/G GAMMA, V?,DTP COmMOn/OaTai/GV, GT,F V, Fr,EV,ET

TYPE OOUbL EG, K, GAMMA

COMMON /PARA PERCEN

COMMON/ISP/N

N = l GO TO 1 ENTRY N H2

H«2

GO TO 1 ENTRY 0 M2

N = 3 GO TO 1 ENTRY E H2

N*4

GO TO 1 ENTRY F H2

N*5

CONTINUE

EH32.015y*

GAHHAs-. Odti

R*. 08205616

G( 1) = 4.614387755654373260330-04

G( 2) = 4. 233l8455bO8b7704344OD-Q2

G( 3) =-5.09b55b2265Q373332l57D-Qi

Gl 4) ~ 2. 923Q597382b95860§34oO4Q0

G( 5) =-2.987609l4721136029849D*01

G( b) = BBiiwosjui 07C 375bobD-j5

G( 7) 32225o95463922652Q67D-Q3

G( 8) = 3. 31650-+43170 189249291O-Q1

G( 9) = 5. G93/'055b085174282592D-*-0i

G(lU) = 1.973828324919J4714O77O-07

G(ll) - 2. 85o492Q39628227170o3O-04

3 (12) -~Z. 22 827 9239 12 3* 8 0 57 0-*5O- 0 2 G (13) =-2.

2

57h8 113b 7b t3 0W a o 9720- Ub

G(14) - 2. 4l427236974bb7590421D-05

G(15) =-1.6 9571339358 8*1 3*7013 0-0

3

G(lb) =-5. 3 93o7b39127 519319151D-0 7 G(17) = 3,998955244328:83*08820-09

G(18) = 1. 1424575bi274-»9354l050-06

G(19) = -1. 2 525be>225 69 bO 5 2 741230- 0 8 G (20) = -4. 9 1 78b 1931* 83 25 39 8 8 29b 0* 01

G(21) =-1. 585bbb0173b8=>7779b97D*02

G(22) =-1. 9 0160 2946 27 21 9 5 54 3660- 01

G(23) = 9. 198C2Q86250 0 5 0 2 7 8 1990* 0

0

G (24) = -3. 18045 5513 61 U'*4498 741 0-04

G(25) = 1. 191057791926527091330-03

G(26» = -3. 7 9 135 277 3 22 5 9 917 6 132D- 07

G(27) =-3. 98 33 7 7699 095395450920-05

6(28) = - 1. Z 3451 i 854 68 8 97 2 90 70 80- 10

G(29) - 1.950266293499369896810-09

G(30) =-2. 3 o 034391 7 10 91 6 9 84 687 0-1

3

G(31) =-*. 07 35 7 66 0819299386 618 0-13

G(32) = 8. 8.1354933 777524867160-12

VP (1> = 3. 0530 0134164

VP(2)=2 . 8081i>9258l3

VP (3) =-0. 5 55*6121056

7

VP(4)=1. 59514439374

VP (5)=1. 5314H54423

VP(7)=13.8

VP (b) =0. 0695

VP (b) = 32 .933

DTP=1./26.17*1000.

RETURN

ENO

35
FUNCTION FINOTtMHJBSI

CHANGED TC WORK FOR HVERQ^p, STARTS AT T-CRITICAL

COHMW/UAi • ' i:ifc''»H«AtWP»DTP

OIHEHSZON G«32H »WP(»I

TYPE OOUBLE&vRtGANNA

T=VP«8J(

00 7 1*1,38

P=VPN(T)

IF(ABS <P-POBS)-.000 001»POBS) 8,8,6

9 CONTINUE

CORR=(POBS-P) /OPOTVP(T)

7 T=T*CORR

3 CONTINUE

FINOTV=T

RETURN

ENO

FUNCTION PMELT(P)

FINOS TEMPERATURE FOR AM IMPUT MELTING PRESSURE

2 DIMENSION PP (77) ,TT ( 77

)

DATA (NTR=1)

k IF(NTR.EQ.2) GO TO 12

5 NTR=2

a T = l<+.000

7 DO 11 1=1,77

9 PP(I)=PRESSM(T)

lu TT|I)=T

11 T=T*Q. 120

12 00 Ik 1=1,77

13 IF (PP(I)-P) 1^,15,15

i<* CONTINUE

15 TAPP=TT(I)

la 00 23 1=1,10

17 T=TAPP

18 PM=PRESSM(T)

19 FUNC=PM-P

20 T=TAPP*0.QQ1

21 PM=PRESSM(T)

22 FUNCP= (PM-FUNC-PI/O. 0Q1

23 TAPP=TAPP-FUNC/FUNCP

Zk PMELT=TAPP

25 RETURN

END

FUNCTION OIE(OP)

DI = DP

CM=. 99575- 0.09069*01*1. 1227 »0I»*2

CM=1./CM

DICM=DI*CM

EP=-(i.*2. »DICM) /(OICM-1.)

DIE=EP

RETURN

END

FUNCTION VPN(TT)

COMMON/DAT A/G»R»GAMMA» V? DIMENSION G(32> , VP(B)

TYPE DOUBLE G,R, GAMMA

T = TT

X=(1.-VP(7)/T»/(1.-VP(7>/VP<8>>

VPN=VP(6)*EXP (VP (II *X*VP(2)»X»X*VP(3)*X**3*VP(<»)*X»(1.-X)»*VP(5)>

RETURN

END

FUNCTION DPOTVP(TT)

COMMON/OATA/G,R,GAMMA, VP DIMENSION 5(32), VP(8)

TYPE DOUBLE G,R, GAMMA

T = TT

IF(TT.GT.VP(8)>G0 TO 1 X=(l.-VP(7)/T)/(i.-VP(7)/VP(8)>

0X0T=(VP(7)/T**2>/(1.-V :>(7)/VP(8»)

DPDT=VP(1)*DXDT*2.»VP(2) * X*OXOT*V P ( 3 ) *3 . *X**2»0X0T* VP ik ) * 1 ( (i.-X) **VP( 5) ) DXOT+VP(i»> X*((l.-X)»*(VP(5)-i.))*VP<5) *(-DX0T) PDT=DPDT* VPN (T) P0TVP=DPDr

RETURN

1 QPOTVP=0

RETURN

END

36
FUNCTION P*ESSM(T>

C CALCULATES MELTING PRESSURE FROM AN INPUT TEMPERATURE

21 PS=.0695+< T-13.8 03) *3Q .3312»EXP (-5.6 93/T)*(T-13.80 3)*2.G*T/3.0

19 PRESS*=PS

RETURN

END

FUNCTION OSATV(T)

DIMENSION GV(8),GL(7)

DATA(&L=0. 04 6645 8130 03 , - 3.

4

7 792 781 86E-2, 4.0 7 7653 8192E-1,

1-1. 171 97 87 304,1. 621392 4:«, -1.15 31 0 96683, 0.3382 54920 39)

OA T A ( RHOC = 0 • 0 313 6),(BETAL=.3478 60 27325),(TC=32.938)

DATA(GV=-Q.0475Q15 71529,3.4871213005E-2,-4.1221290925E-l,

11.56665985 5 , -2. 8 0611*27339, 2. 710 54 55626,-1. 30 74 7 7 359,

2 0.229212859 2 2) , ( BET A V = .

3

k 8 31 23762 5 ) ,< FACT=496. 0 4651)

A= (TC-T)/T3

3V = RH0C*G\/ <1> *A* *BET AY OO 1 1=1,7

1 l)V=QV*GV(I*l)*A*Mi. + (I-l> /3.)

OV=DV*FACT

OSATV=DV

RETURN

ENTRY OSATL

A= (TC-T) /TC

0V = RH3C*GL( 1

) *A**BET AL 00 2 1=1,6

2 OV=OV*GL (1*1 ) »A*» (1. ( 1-1) /3.

)

OV=DV*FACT

OSATV=DV

RETURN

E NO

FUNCTION C?I(T)

COMMON/CPIQ/G(ll) COMMON/ISP/N

IF (N.NE. 0) 30 TU 5 K=l

1 U=G(9)/T

EU=EXP (U)

TS = 1./T**<+

GO TO (2,3,h) ,<

2 CPI=G(8) *U»U*EU/ (EU-1. ) **2

00 10 1=1,7

TS=TS*T

1U CPI=CPI+G 1 1

) *TS

CPI=CPI*8. 31434

RETURN

5 CPI=CPO(T,N)

RETURN

ENTRY SI

IF(N.NE.0)GO TO 6 K=2

GO TO 1

i CPI=G(6)*(U/(EU-l.)-ALOG(l.-l./EU)) _ #

i

l-G(l)»TS»T/3.-G(2)*TS»T*T/2.-G<3)/T*G(4)*AL0G(T)+G(5)»T-»-G(6)»T»T/2

2.+GU) *T*»3/3.

CPI=CPI*8. 31434+G <li)

RETURN

a CPI=CPOS(T,N>

RETURN

ENTRY HI

IF(N.NE.0)GO TU 7 K=3

GO TO 1

4 CPI=G(6)*U»T/<EU-i.)-G<l>/<2.*T»T)-G<2)/T*G<3>»AL0G<T)*G<4)*T

1+G (5) +1*1/2 . *G(6)*T**3/3. G < 7) T* *4/4

.

CPI=CPI*8.31h3^+G(10)

RETURN

7 CPI = CPOH(I",N)

RETURN

E NO

FUNCTION JILV(T)

COMMON/0 AT Al/G\/,GT,FV, FT ,£V,ET

J I MENS ION GV ( 9) , GT (9 ) , F\l ) * FT (4) , EV< 8) , ET(4)

SUM = 0 TF = T** ( l ,/3 . » TFF=T»»<-4./3.)

00 10 1=1,9

T FF = TFF* TF

10 iUM = SJM + G\/(I) *TFF

OIL V = SUM*1« l . . RETURN

ENTRY DILT

TF = T*» (1./3 . I 37
TFF=T»»<-4./3.)

SUM»0

00 20 1*1,9

TFF«TFF»TF

3 SUM=SUM+GTCI)*TFF

DILV=SUM

RETURN

END

FUNCTION CPO(TI,N)

DIMENSION T 1 58) * CPP i 58 > , CPNf 58 » ,C PO ( 58) , CPE 158

1

COHMON/PARA/ PER^EN

CALCULATES IOESL GAS SPECIFIC HEAT FOR H2 BY INTERPOLATING

OATA TAKEN F ROM RP 1932, UNITS OF THE TABLES ARE CAL/MOL DEG 5.

UNITS OF OJTPUT ARE JOULE S/MOL OEG K. THE INDEX N DETERMINES THE

SPECIES, FOR N*l, PARAHY0R0GEN«N«2 NORMAL » N=3 ORTHO. N*4 EQUILIB N=5,S0ME ORThO-PARA MIXTURE SPECIFIED BY COMMON /PARA/ , PERCENT

RANGE OF T E MP IS FROM 10 TO 5000K.

DATA (T*

1 10. 0, 12. U, 14.0, Ib.O, 18.0, 20.0, 25.0, 30.0, 35.0, 40. M, 45, U,

2 50.0, 55.0* 6U.0* 55. 0, 70.0, 75.0, 80.0, 85.0 , 90.0, 95.0*100.0*

3105.0,110. 0 , 115. 0,12 0. 0,1 25. 0*130. 3,135.0*140. 0,14$. 0,150 .9, 160. 9,

417 0. 0 * 180. 0. 19U. 0*20 0.0, 21 0.0, 220. 0,2 30. 0,240.0, 259.0,260.1). 270. 0^

5280.0,290. 0 ,3.iu. 0,35 li. 0,40 0.0 ,500 . 0,600. ,70 0. ,1000. ,1500. ,20 00.

63000. .WOOD. , 5000.)

OATA ( ( CPE ( I ) ,1 = 1 ,58) =4.968 .4.96884.4.97647,5.01153,5.07451,5.21)6

11 ,5.83508, 5. 8128 2,7. 07 989, 8 .60613 ,9. 0 0 231. 9. 0 80 05, 8 .93278 ,8. 658 94,

28. 3360 3. 8. 012G7, 7. 71 00 9, 7. 4416, 7. 2110 9, 7 .01858, 6. 85857, 6. 72557,6.6

32 055,6. 53555, 6.46904, 6.420 03, 6. 38 40 3, 6. 36151 ,6„ 34602,6. 33753, 6. 340

401,6.34577 ,6. 372 76, 6. 4 13, 6. 45925,6. 50 975, 6. 5605, 6. 6395, 6. 65724,

7 6. 692,6. 734, 6. 7.71,6. 604, 6 .8 32,

56. 856, 6. 87 7, o. 895, 5. 95 C , S . 9/4, 6. 9 93 , 7 . 2 0 9,7 . 336 , 7. 219 , 7 . 720, 8 .199

,

68.859,9.342* 9.748)

DATA(CP0«

14.968,4.93 8, 4.968,4. 96 8, 4 . 468 ,4. 968 , 4 . 96 8,4 . 968 , 4. 968 ,4. 968 ,4. 961,

24.966,4.968,4.96 8,4. 969,4.972,4.975,4.981,4.990,5.082,5.018,5.039,

35. 06*, 5. 094, 5. i29, 5. 16 9,5.213,5.261,5.313,5.369,5.427,5.487,5.612,

45.741,5.86 8 ,5.992,6. 10 9, 6 . 219,6. 320 , 6 . 411 ,6 . 493 , 6 . 566, 6 . 629 , 6. 684

,

5 6.732,6.77 3 ,6. 80 8, 6. 91 7, S. 962,6.99 3, 7. 0 09, 7. 036, 7. 219, 7. 720, 8. 195,

68.859,9.342,9.743)

DATA(CPP=

14.968,4.96 8,4.96 8,4. 96 8, 4 . 968 96 3 , 4 . 9o 6 ,4 . 966 , 4 . 969 *4 . 972, 4.983

,

25. 0 06, 5. 0*8,5.114, 5. 20 7, 5. 326, 5. 475, 5. 646, 5. 835, 6. S36, 6. 245, 6. 454,

3b, 659, 6. 85"*, 7 .. 37 ,7. 20 3, 7. 351, 7. 480, 7. 590, 7.681, 7. 753,7. 807, 7. 8 70,

47.883,7.85 8, 7.00 8,7.74 2,7.66 7,7.5 91,7.516,7.445,7.380,7.322,7.2 70,

57 .225* 7.1fto* 7.15 2,7. 05 0,7. 010,6.996* 7. 010 ,7.037, 7.219, 7 .720* 8.159,

68.859,9.342 , 9.740)

OATA(CPN=

14.968,4.90 8 ,-+.96 0, 4. 96 6,4. 96 8*h.96 8* 4. 96 8,4. 96d ,«. 966. 4.969 ,4. 972,

24.977,4.988, b. 00 5, 5. 029,5. 061,5.100,5.147,5.201 ,5.261,5.325,5.393,

35.463,5.534,5.o0 6,5. 67 7 , 5 . 74 8 ,5. 816 , 5 • 882 ,5 . 947 , 6« Qtid , 6. 067, 6. 177,

46.276,6.36 6 ,6.446,6. 517,0.581,6.638*6.687,6.731 ,6.769*6. 802,6.831*

5 6. 855, 6. 8? 6 *b.894,6. 95 0 , 6 . 974 , 6 . 9 93 , 7 . 0 0 9*7 . 036 , 7. 219 , 7. 720 , 8. 195

,

68.859,9.342,9.748)

60 T0(l,2,3*4,5) ,N

1 CPO=ATKINT(TI,CPP,T,58, 6 * NES , . 01 ) *4 . 18*

RETURN

I CPO = ATKINT ( TI ,CPN»T, 59 , 6 , NES , . 01 ) *4 . 184

RETURN

3 CP0=ATKINT(TI,CPQ,T,58, 6 , NES, . 01) *4 . 18«r RETURN

4 CP0=ATKINT ( TI ,CPE,T, 50 , 6 , NES, . 01 ) *4 . 164

RETURN

5 TUP=TI*.5

TDN=TI-,5

HUP = CPOH <T JP , 5) HCiN=CP0H(T3N,5)

CPO=(HUP-HON)

RETURN

ENO

FUNCTION CPOH(TI,N)

DIMENSION T (5 8) , HP 15 8) V HN(58 ) ,H0(58) ,HE (58) COMMON/PARA/ PERCEN

CALCULATES THE ENTHALPY OF THE IDEAL GAS FOR H2 BY INTERPOLATION

OATA TAKEN FROM RP 1932, UNITS OF TABLES ARE CAL/MOL

UNITS OF OJTPUT ARE JOU.E/MOL. THE INDEX N DETERMINES THE SPECIES

SPECIES, FOR N=l,PARAHYDR05EN*N-2 NORMAL, N=3 ORTHO, N=4 EQUILI8

N=5,S0ME ORTHO-PARA MIXTURE SPECIFIED BY COMMON /PARA/ , PERCENT

RANGE OF TEMP IS FROM 10 TO 5000K.

DATA ( T=

1 10. u, 12.0, 14. Q, Ib.O, 18.0, 20 . 0 , 25.0, 30.0 , 35.0, 40.0, 45. 0,

2 50.0, 55.0, 60. Q, 65.0, 70.0, 75.0, 80.0, 85.0, 90.0, 95.0,100.0,

3105. 0,110. 0,1 15. 0,120. 0,125. 0,130 . 0,1 35. 0,140.0, 145 . 0,15 0.0, 16 0.0, 4170. 0, 180. 0, 190. 0, 20 0.0, 210.0, 220. 0, 230. 0,240.0, 250. 0,26 0.0 ,270. U,

38
5280. «, 290. 3, 30.. 3, 353. 0,400. 0,500. 0,600. ,70 0. ,1030. ,15 00. ,2000.,

63 0 0 0 . ,40 03 . DATA (HN=

1 303.57

1 427.86

3 628.02

f 8*1. t0

51107.22

61626.93

72730.94

823230.9

0 ATA (HP

5000 .

)

49.68

173.88

378.22

628.14

969.04

61586.36

72729.17

Q23230.9

D AT A ( HO

1 388.33

2 512.53

3 711.29

4 912.49

51153.27

6i6h3 . to 72731.52

823230.9

DATA (HE

1 49.68

2 195.77

3 538.05

4 815.54

51102.11

61626.75

72730. 94

823230.9

0 TO(l

313. 60

452.71

653. 6t

86 9 . 61

1169.49

16 3 4 • 44

3429. h6 .51

323-+5.,-rl895 .) 59.61, 69.55

198.73, 223.61

4Q0.OI, 434.71

663 . 75, 703 . 14

1047.84,1126.58

1550 .49, 1733 .99

3429 . 24, 4129 .48

32345. ,41395 .) 338

537

736

938

1210

1735

3429,

2'

37

la

12

D4

76

.3 32345.,

5^.62

236 .90

575.93

323

477

679

898

1232

1762

129

.54

.56

.51

.17

.71

.30

408

562

761

9o4

1268

1771

4129

1895

69

280

612

879

1233

1762

4129

647. 76

1166 . 03

1594.32

342 9 • to

32345. ,41695

2,3,4,5) ,N

CPOH=ATKlNT (TI,riP,T,

RETURN

CPOH=ATKINT

(

TI,hn,T,

RETURN

CPOH=ATKInT (TI,rlQ,T,

RETURN

CPOH=AT<INT (TI,hE,T,

RETURN

PERCENT=PERCEN /100.

CPOH= I A TKl ,T ITI.HO.T

1ATKINT (T I, MP , T,po, o, RETURN

END

.20

.21

.11

.18

.09

.74

.52

.1.57

.97

.56

.77

.39

.23

.51

.)58 , 58,

58 , 58 ,

,58

NES

333.48

502.43

705 .66

927 .08

1296.78

1 330 .,o

4 8 31 .66

79.49

248.58

464.38

737,23

1204.93

1806 .95

4831.65

418 .14

587.05

786.09

990.37

1327.39

1338*32

4 8 31 • 6 a

79.56

326.24

646 . 13

311.63

1295.23

13 30 .43

4831 .66

343.41

527. 34

732.13

956.33

1361. oi

1898 . 91

6 965 .23

83.42

273.71

495.09

774.92

1282.69

1879.42

6966.23

428.07

611.89

811.14

1015 . 60

1387. 91

1905 .41

6 966*23

89 . 66

371. 33

682 . 82

943.40

1363.58

1896.36

6966 .23

353.35

552. 32

758.92

985. 91

1427.10

1967.57

10697.2

99.36

299. 11

525.64

613.10

1359.74

1951.47

10697.2

438.31

636. 73

636.28

1043.51

1449.56

1972.94

10697.2

99.96

415.35

716. 78

975.11

1426.43

1967. 55

10697.2

378.19

577.40

786. 16

1015.30

1493.20

2036.43

14079.2

124.20

324.90

559.62

851.69

1436.03

2023.15

14679.2

462.85

661.57

861.54

1070.50

1512.26

2040.86

14679.2

127.50

437.86

750.14

10G6.80

1492.76

2036.42

14679.2

6,NES,.01)*4.184

6, NES,. 31) *4. 184

5, NES, . 01 ) »4 . 104

i , NES, . 01) *4 . 184

6,NES,.01)*( 1 .-PERCENT)

4

. 01) ^PERCENT) *4 . 184

433.03 , 602.62,

813.55,

1045 .99,

1559. 84,

2382.74

,

143.04,

351.22,

593.41,

690 .60 , 1511. 56 , 2377.83,

487.69 , boo .42 , 806.93 , 1097 .78 , 1575.93,

2384. 33

,

159.12,

493 .74 , 733.03,

1038.52,

1559.55,

2332 . 74 ,

FUNCTION CPOS(TI,N>

DIMENSION T(6G),SP(5C),SN(bO),SO(60),SE(60)

COMMON/PARA/ PERCEN

C CALCULATES ThE ENTROPY Or TnE IDEAL GAS FOR H2 3Y INTERPOLATING

C OAKEN FROM RP 1932, UNITS OF ThE TABLES ARE CAL/HOL OEG K. 1

C UNITS OF OJTfUT ARE JOUlES/mOL DEG K. TmE INDEX N DETERMINES TnE

C SPECIES, FOR N=l,PARAnTDR03EN,N=2 NORMAL, n=3 ORTnO, N=4 E0JIlI3

C N=5,S0ME ORTHO-PARA MIXTJRE SPECIFIED 3Y COMMON /PARA/ , FERCE <T C RANGE OF TEMP IS FROM 1] TO 5303<.

0 AT A I T =

1 10.3, 12.0, 14.0, 16.3, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.3,

2 50.0, 55.3, bu.O, 55.0, 73.0, 75.3, 80.0, 85.3, 93.3, 95.0,100.3,

3105. 0,110. 0,1 13. 0, 12 CO, 125. 0,130. 0,1 35 . 0,143. 3, 145 .0,153.0,163.3,

4170. u, 180. 3, 19. . 3, 231. C, 210. 3,220. 0,230. 3, 243. 3, 253. 0,26i.. 0,273. 3,

5260. 3, 290. 3, 303. 3, 35C. 0,400. 0,450. 0,500. 3, 553. 3, 600.0,700. ,1000.,

61500. ,2003., 3300. ,43CQ.,5Q0C.)

0 AT A ( i P = Ul 1.214,12. 12.',12. 3 35, 13. 549, 14.135

116. 1Q2, 15.5oo, 19. 21 4, 19. 593, 23. 13 5. 02 2. 017 , 22. 356, 22. 633, 2 3. CI* , 23. 334

42 4. 043 , 25. 132, 25. 413,25. 551, 25. 94 5. 526. 231, 28. 577, 28. 932, 29. 253, 29. 586

630.963,31.212,32.305,33.244,34.069,

739. 701, 42. 723 , 45. 337 , 51. 221, 53. 833

OATA(SO=

115. 551, 16. 486, 17. 252, 17. 916, 13. 50 1 22 2.46o, 23. 353, 23. 576, 24. 050,24.48 2. 326. 215 , 26. 50 0 , 26. 771, 27. 029,27. 275

428.375,28. 575,28.769,28.958,29.143

530.308,31. 109,31.401,31. 584,31.95 3 633.216,33.446,34.535,35.432,36.253

741. 884 t 44. 90 3, 47. 190, 5 1.221, 53. 839

14.653, 15 . 766 , 16 . 672 , 17. 433,

20. 546, 20. 935, 21. 310 « 21. 669,

23. o-* 3,23. 957 , 24. 260 , 24. 3 57,

26. 451, 26. 92 9 , 27. 379, 27. 832,

29. 869, 33. 177, 33. 452, 30. 716,

34. 606,35. 473, 36.332,37.165,

55.969)

13. 024, 20. 133 , 21. 033, 21. 33-*,

24. 880, 25. 248 , 25 . 591 , 25. 912,

27.512,27.739,27.95 8,28.170,

29. 502, 29. 84 6, 33. 177, 30. 498,

32.225,32.43 4,32.73 5,32.979,

36.989,37.55 6,38.265,39.3 48,

55 . 969)

39
DATAISN*

115.607,16. S12,17.27S,17.9<*l«ld.527,19.050»20.159t|l.a$4> ( 21.|3&,

222. %9<» s 23. 0 79,23.603 0 78 , 2«|. 51 3 . 2% . 91%, 25 . 28 8 tiS »l?.f6§,

326. 263, 26. 562, 26. 668, 2 7. li»3, 27\ %0 7, 27. 663,27. 911, 28.151 ,26. 36%,

"2 8. 611 ,28. 832 ,29. 8%7 ,29. 256, 29. %6 1,29. 856, 30. 23%, 30.595, 30. 9%2,

531. 27%, 31. 59H,Ji.901 ,32. 197, 32. %8 3, 32. 758, 33. 025, 33. 282 ,33. 531,

733. 772, 3%. 0 0 5, 35. 073, 36. 0 0 3, 36. 825, 37. 561. 38. 228 ,38. 836, 39.920,

7<,2.%55,%5.<»75,*7.762,51.221,53.83 9,55.969J

DATA(SE=

111.215,12.120 ,12.887 ,13. 55%, 1%. 1<*9 ,1% .692,15. 918 ,17. |69,1|. 196,

219.29^,20.331,21.285,22.1%5,22.911,2 3.592,2%.19B,l%.7%0,25*i|ffl

325. 67<», 26. 3 8 G, 26. %55, 26. 6 0%, 27. 12 9, 27. %35,27. 72%, 27. 999,28. 2&0.

%28. 510, 28. 75 a, 28. 980, 29. 2 03, 29. %1 8, 29. §28, 30. 216, 36 .56%, 30 .93%,

531. 269 , 31. 59 1,31. 899, 32. 196, 32. %8 2, 32 .758,33.62%, 33. 262, 33. 531, 633.772,3%. 00 5 ,35 . 0 73 , 36. 0 0 3 , 36 . 625 , 37 .561.38 . 22 8 , 36 .636,39.920

,

7%2.%55.%5.%75,%7.762,51.221,53.839,55.969)

GO T0(I,2,3,%,5> ,N

1 GPOS=ATKINT(TI»SP»T»60» 6 , NES, . 01 » % . 18%

RETURN

2 CPOS=ATKINT(TI,SN,T,60, 6, NES, • 01 ) •%• 18%

RETURN

3 CPOS=ATKINT<TI,SO,T,60, 6, NES, .01)*%. 18%

RETURN

<i CPOS=ATKINT (TI,SE,T,60 , 6 NES , . 01 ) *< . 18%

RETURN

5 PERG£NT=PERCEN /100. S PER=1 .-PERCENT

CPQS= (ATKINT(Tl,SO,T,60,6.NES,0.0il»PER»ATKINT«ri,SP.T,60,6,N£S,.0

lIl»PERCENTf *%.18%-C8.317* (PERCENT»ALOGf PERCENT) PER'A LOG tPER) >

)

RETURN

ENO

FUNCTION FIND T(P,0)

NEM FEB 1975

CONHON/OATA/G ,R, GAMMA* VP, OTP

OIMENSION G(32),VP(8)

TYPE OOUBLEG,R, GAMMA

PP = P 00 = 0

USES A FIRST GUESS IN TEMPERATURE FROM Tl

CALL T1(PP,00,TA)

TT = TA DO 10 1=1,10

CALL PRESS(PP,00,TT)

P2=PP

IF(ABS (P-P2 ) -l.E-7*P) 23 ,20,1

1 CALL QPDT ( PP , DO, TT

)

OP = PP

CORR= ( P2-P ) / OP

IF(A8S (C0*R)-l.E-5 )20,20,10

10 TT=TT-CORR

20 FIND T=TT

RETURN

END

FUNCTION FIND D(P,T)

CHANGED, T*IAL OENSITY VIA SUBROUTINE RH01/T1, 2% FEB 1975 COMMON/OAT A/ G ,R» GAMM A, V

3

, DTP

DIMENSION 3 (32) , VP(8)

TYPE DOUBLEG.R, GAMMA

TT = T CALL RH01(P,D,TT)

00=0

00 10 1=1,50

CALL PRESS(PP,OD,TT)

P2 = PP

IFU3S (P-P2)-1.E-7»P)20,20,1

1 CALL DPOD(PP,OD,TT)

DP = PP

C0RR=(P2-P) /OP

0 = 00

IFUBS (CQRR) -l.E-7»0> 20,20,10

10 OD=DO-CORR

20 FINO 0=DO

RETURN

ENO

SUBROUTINE PROPS (PP, DO , T T

)

DIMENSION X (33)

OIMENSION 3( 33) ,G(32)

COMMON /OAT A/ G,R, GAMMA

TYPE DOUBLE B »G , 0 , T ,P , 02 , 03 , 0% , 05,06 , 07 ,08 ,09,01 0 .0 11, 012

If D13,TS,T2,T3,T%,T5,F ,Fl,F2i,F22, F23,F2<t,F2$,F 26, GAMMA,

R

TYPE DOUBLE F 21 2 , F222 ,F 232 , F2%2 , F 252 , F262

TYPE DOUBLE G 1 , G 2 , G3 , , &5 , G6 , X40
EQUIVALENCE <B,X>

OATA(IO-l)

DATA(IZ=1>

C PROPS FOR H2 USING THE STE MART- JACOBS EN EQUATION

1 CONTINUE

IFUZ.LE.O) GO TO 2 IZ«0

2 CONTINUE

C PROPS FOR METHANE USING STEW ART -JACOBS EN EQUATION OF STATE

C PRgjjIHINAR* FIT - DENSE FLUID REGION EMPHASISEO,MCCARTY,<»/26/73

P«PP

T*TT

GM«GAMMA

02*D»D

03*02*0

0<**03*D

D5»Di»»D

Db»05»D

07*06*0

D««O7»0

09«oa»o

D10«O9«O

Dli*O10*O

Di2*01i»D

¥i=6?olT?T)

T2»T»T

T3=T2»T

T<»«T3»T

T5»T*»»T

F=DEXP!GAHHA*D2)

GO TO <1Q0 , 200, 300, "»00 , 500 ,600,700) ,K

ENTRY PRESS

K*l

GO TO 1

100 CONTINUE

B( 1)=02»T

B( 2)=D2»TS

B( 3) =02

B( M = D2/T

B( 5>=D2/T2

3( 6)=03»T

8( 7) =03

B< 8>=D3/T

B< 9>=D3/T2

B<iQ)=0<**T

B(ll)=01t

B(12) = D<*/T

B(13)=D5

8 (m) = D6/T

B(15)=06/T2

B(16)=D7/T

B(17>=D8/T

B(18)=D8/T2

B<19»=D9/T2

B(20)=D3»F/T2

B(21)=03*F/T3

B(22)=05»F/T2

B(23> = D5*F/T<»

8(2<»>=07»F/T2

B(25)=D7*F/T3

B (26>=D9*F/T2

8<27) = D9*F/T->

B(28)=01i»F/T2

B(29)=011»F/T3

B(30>=D13»F/T2

B(31)=013»F/T3

3<32>=D13»F/T<*

IF (ID.GT.O) GO TO 102

8<33)=P-R»3*T

RETURN

132 P=0

M = 3 2 00 101 1 = 1,

M

101 P=P*B(I>»G{I>

P=P+R»D*T

PP = P RETURN

ENTRY DPDD

K = 2 GO TO 1 200 CONTINUE

F1=2.QO*F»SM*0

F21=3.0JG*F*O2 <-Fl*33

F22=5. 000»F*Ch *F1*35

F23=7. OOC*F*06 *F1*D7

41
F2*» = 9. 000*F»O8 -*F1»09

F2 5=11.00»F»Q10*F1*011

F 26=13. 0C*F*D12*F1*D 13

B( 1)=2.00»D*T

B( 2>=2.00*-O + TG

3( 3)=2.00»O

8( <,) = 2.0Q*O/T

8( 5>=2.00»0/T2

8( 6)=3.QQ*02*T

d< 7)=3. 00*02

3( 8)=3.00*D2/T

B( 9)=3.0u»O2/T2

B(10)=<*.00*O3»T

B(11)^.00*D3

B(12)=*t.0Q*D3/T

6 ( 13) = 5.QQ»D<*

B(l<+)=6.00*O5/T

B (15)=6.00*D5/T2

3(16>=7.00»U6/T

B(i7>=8.00»O7/T

8<18)=8.00*0< /T2

3(19)=9.00»Oo/T2

8(20>=F21/T2

8<2l)=F2i/T3

3 (22) =F22/T2

8<23)=F22/TW

8<2<t) =F23/T2

8<25>=F23/T3

8 (26) = F2<*/T2

B<27)=F2'»/T4

8 (28) =F25/'T2

8 (29) =F25/T3

B(30)=F26/T2

B(31)=F26/T3

8(32) =F26/T<»

H = 3 2 IF (ID.GT.O) GO TO 202

B(33)=P-RM

RE TURN

202 P = 0 00 201 1=1,

M

201 P = P<-3< I) &( I

)

P=P+R*T

PP = P RETURN

ENTRY DPDT

K=3

GO TO 1 33J CONTINUE

X( 1)=D2

X( 2) =D2/( 2. Uu*TS)

X ( 3> =0

X( <*)=-D2/T2

X( 5) =-2.0 0 *U2/T3

X( 6)=03

X ( 7) =0

X( 8)=-D3/T2

X( 9) =-2.UQ*'j3/T3

X (1„) = D<* x 1 1 1 ) = o

X ( 12) = -0**/T2

X (13)=D

X (1^)=-D6/T2

X (15)=-2.00»G6/T3

X ( 16) =-07/ T2

X (17)=-U8/T2

X (lo)=-2.00*06/r3

X ( 19) =-2. J 0 * U9/T

3

X(20)=-2.0U*U3-F/T3

X <21)=-3.0Q» 23*F/T<*

X(22)=-2. 00*U3*F/T3

X (23) =-<. I j*J5*F/T5

X (2«t) = -2 . J j * 07*F / T3

X (2'3)=-3.0 0»u/*F/Ti4

X(2o)=-2.00 4 U9*F/T3

X (27) = -<+.00*O9*F/T5

X (28! =-2.00 *Q11*F/T3

X ( 29)=-3. 0 0 + L)ll*F/T4

X (30)=-2.00* J13»F/T3 X(3D = -3.00*U13*F/Tt»

X(32>=-<t.00*D13*F/Ti

1 F ( I . G T . 0 ) G U TO 302

X ( 33) = PP-R* 0 RETURN

3 0 2 P = o DO 301 1=1,32

301 P = P*G( I) *X( I)
PP«P+*»D

RETURN

ENTRY DSDN

GO TO 1 -3D CONTINUE

C S«S|-R»AL0S«D»R»T/P01*COSON«O)-OSONI0l)»i01.32S CPOS(T)

G2*(F*D2-2. 0 0*Gl)/(2 .QQ»&M) G3*(F»Dh-*. 0 0*G2) /(2 . 0 0»SM>

Sk* (F* 06-6.0 :*G3)/ (2 .0 C*GH) G5*(F*O8-3.0Q»Gi»)/C2.Q0*Sl>

G6=(F*D1Q-10 . Q0*G5) / (2 . 3 0*GM>

X( t)=-0

X( 2)=-D/(2.0G*TS)

X( 3)=Q.0u

X( h)=+Q/T2

X{ 5)=2.QQ»D/T3

X( 6)=-O2/2.00

X( 7)=0.C3

X( d)=02/(2. 00*T2)

X( 9>=D2/T3 XUO)=-03/3. go X(il)=0.03

X(12)=03/(3. Q0*T2) X<i3)=Q.D0

X (i<*>*05/<5 . Gu*T2)

X(i5)= 2.00»Q5/(5.03»T3)

X(16)=O6/(6.00+T2)

X(l7)=07/<7. Q0*T2) X(i6l=2.00»D7/(7.00*T3)

X<19)=08/< +. QQ*T3)

X(20)=2.0J»G1/T3

X(21}=3.00»G1/Tit

X (22)=2.00»G2/T3

X (23)=f».00»G2/T5

X (2<»)=2.00»G3/T3

X(25) = 3. 00*G3/Tt»

X (26) = 2.G0*G<*/T3

X (27)=<*.00»Gh/T5

X (28) =2. 0Q»G5/T3

X<29) = 3.03*G5/Tt*

X (30)=2.QQ»G6/T3

X(31)=3.03»G6/Ti»

X (32)=*».00*G6/T5

IFUD.GT.Q) GO TD *Q2

RETURN

h02 P=0

DO hOI 1=1,32

<*01 P=P*G(I)»X(I)

PP=P

RETURN

ENTRY DUDN

K = 5 GO TO 1 33a CONTINUE

C H=H0*(T*OSON(D»-OSOM(an»101.325*IOUDN(O-OUDNI0) )*101.325+CPOHCT1

C (P/0-R*T) 101.325

G1=F/(Z.Q0*6M) G2=(F*02-2.u- *G1 / 1 2 • 0 0 * GUI

G3= (F»Dt-i». 0 U-»G2 ) / (2 .0 0*GH)

G<*=(F»D6-5. 0>j»G3)/(2.0Q»3M> G5=(F*08-3 . 0 G*G* ) / (2 . 0 0* 31)

G6=(F»D10-10.00»G5)/'(2.JO»GMJ

X( 1)=D*T

X( 2)=D*TS

X ( 3) = D X( <*) = D/T

X( 5)=0/T2

X( 6)=D2»T/2.GG

X( 7)=D2/2.0u

X( a)=02/(2

.

U 0 + T)

X( 9) =02/( 2 . 00*T2)

X(lw)=0 3*T/3.J.

X(ll)=O3/3.00

X(12)=D3/(3. 0G*T) <(13)=0Wf . 00

( 1<+) =D5/ ( 5 . uQ*T

)

U5)=D5/(5.0G*T2)

16) =D6/ ( a . 00 + T

)

X (17) =D7/ (7 . 00*T)

X (18)=D7/(7 . 00+T2)

X(19)=D8/(8. G0*T2)

X (20) =Gi/T2

X (21)=G1/T3

X (22) =G2/T2

X (23) =G2/T*

43
X (2<*)*G3/T2

X (25>=G3/T3

X(26)=GWT2

X(27)=G«*/T<»

X (28>=G5/T2

X <29>=G5/T3

X(30)=G6/T2

X (31»=G6/T3

X (32)=66/Tl»

IF(ID.GT.O) GO TO 502

RETURN

502 P=0

00 501 1=1,32

501 P=P+G(I)»X<I>

PP=P

RETURN

ENTRY TOSOT

K-b

GO TO 1 cOO CONTINUE

CV = CVOMTCSDN(/) -TDSDNO) ) •101.325

Gi=F/(2.00»GM)

G2= (F»D2-2. 0 0*G1)/(2.Q0»GM)

G3»(F»Qif-<». 0 0*G2»/(2.0 0*GH)

G*»*(F*D6-6. Q0*G3)/C2.0 0»GM)

G5=(F*D8-8.0 0*G<»>/<2.0 0*GM>

G6»(F*D10-10.00*G5) / (2.0 0»GM)

X( li =0.00

X( 2»=-D/li».00*TS)

X ( 3) =0 .00

X( <*)=2.00*D/T2

X ( 5)=6.00»0/T3

X( 6)=O.O0

X( 71=0. DO X( 8)=D2/T2

X( 9)=3.0Q»D2/T3

X(10)=0.O0

X(11)=0.DQ

X (12)=<2.00»D3)/ (3.0 0»T2>

X(13)=0.D0

X(i<») = (2.0 0*D5>/<5.0 0*T2) X(15)=(6. 00*05) /(5.00»T3)

X <16) = 06/(3 . 00»T2)

X(i7>=(2.0Q*D7)/(7.0 0*T2>

X(18)M6.Q0*D7)/ (7.00*T3)

X(19) = l3.00*O8>/!<».00»T3)

X(20)=6.000*G1/T3

X (21) = 12.QQ»G1/T<*

X (22>=6.000*G2/T3

X(23)=20.00»G2/T5

X (2<*)=6.000*G3/T3

X (251=12. 00*G3/Ti»

X (26)=6.0Q0*GWT3

X <27)=2Q.00»G<+/T5

X (28)=6.CUQ*G5/T3

X(29) = 12.QQ*G5/H»

X(30)=6.000*G6/T3

X <31)=12.0C*Gb/Ti»

X (32>=2Q.00 # G6/T5

IF(IO.GT.O) GO TO 602

RETURN

632 P = 0 DO 601 1=1,32

.01 P = P+G ( I) *X ( I

)

PP = P RETURN

ENTRY 0P2D2

K=7

GO TO 1 703 CONTINUE

Fl=2.»F*GM*0

Fi2=2.»Fl»SM+D+2.»F*SM

F212=3.»F1»D2*3.*2.*D»F*F12*03*F1*3.»D2

F 222 = 5. *F1* Qt* *5 .*<. »D3»F*5. *0«»*F 1 + F12*05

F232=7.*Fl*Dfa*r. *6.»D5»F*7 .06»F1*F12*D7

F2<,2 = 9.*Fl*D8f9.*8.*D7»F*9.*08*FH-F12*09

F252=11.»F1*D10*10. # 11 .»D9*F*11.»010»F1*F12*011

F262=13.*F1*L12*13.*12 .*Dii*F*13. *012»F 1 *F1 2*01

3

B(li=2.*T $B(2)=2.*TS $ B(3)*2.

8(<+)=2./T $ B(5)=2./T2 $ BC6)=6.*D*T

d(7)=6.*0 $ B(8)=6.*D/T $ B(9)=6.*0/T2

8<10>=12.*32*T $ B(lli=12.»D2 $ 8(12)=12.*02/T

B(13)=20.*D3 S B(H») =30.*0<»/T $ B ( 15 ) =30 . *0«»/T2

B(16)=*»2.*D5/T $ B (17) =56.*D6/T $ B ( 18) =56. *D6/T2

B(19)=72.»37/T2 $ B( 20 ) =F212/T2 S B ( 21 » =F21 2 /T3

8 (22)=F222/T2

8(23) =F222/T^ $ B( Zk > =F232/T2 $ B ( 25 ) =F232/ T

3

44
Ms 32

IF(ID.GT.O) GO TO 702

702 P=0

00 701 I=l LH 701 P*P*BII)*GCIT

ppxp

RETURN

END

FUNCTION CP<D,T)

CVEE=CV(D,T)

CALL DPOT(DPT,0, T)

CALL OPOD(DPO,D,T)

CP*CVEE-MT/ (D**2) »<DPT»»2) /DPO) 101. 325

RETURN

FUNCTION C\MO,T)

0ATA(R»6.3l«»3i»)

00=0

TT = T CALL TOSOT(CO,OD,TT>

DO=0

CALL TOSDT(C0.D0,TT)

GV»CPI <TT) (CQ-CD)»101.325

CV=CW-R

RETURN

ENO

FUNCTION ENTROP(O.T)

Rs. 08205615

00 = 0 TT = T CALL DSDN(SD,OO t TT)

OO*0

CALL OSONC30 ,00, TT)

ENTROP»(SO-S0)*101.325-R»ALOG(D»R*T) »101 .325+SI ( T)

RETURN

END

FUNCTION y/ISC{QO,T» COMMON/CRIT/EM

O = OO*EM/100 0 . VISC=0ILV(T> FOCV(D,T) EXCESVCDtT)

RETURN

ENO

FUNCTION FDCV(0,T)

COMMON/ DAT A 1/GV,GT,FV,FT,EV,ET

DIMENSION GV(9) V GT(9) f FV(<») .FT(<») f EV(8) ( ET<<»)

FDCV=(FV<1> Fvm* (FV(3) -ALOG(T/F ) ) »*2) *0

RETURN

ENTRY FDCT

FDCV=(FT (1H-FT(2)MFT{ 3) - ALOG ( T/F T

(

k

)

)

)

**Z) *U

RETURN

END

FUNCTION EXCES\/(0,T>

COMMON/ DAT A 1/G\Z,GT»FV,FT , EV,ET

DIMENSION GV (9) , GT (3 ) , Ft/ Ik) , FT Ik) , EY (8) , ET Ik)

R2=0**(.5) * C (D-EV(8) )/EY (8)

)

R=D**< .1)

X=EV(1) + EV (2)»R2*EV(3) •R+iVik) »R2/(T»T) +E V ( 5

) *R/ J** (1 . 5 ) *EV I 6 ) /

T

l*EV<7) »R2/T

X1=EV( 1) *EV <b)/T

EXCESV=EXP (X)-EXP (XI)

RETURN

ENTRY EXCEST

R=D**(.l)

X=ET(1) +ET (2) »R+ET (3) *R/T»<M1.5) +ET Ik) /T

X1=ET(1) +ET <<*)/T EXCESV=(EXP (X)-EXP (Xl))/10.

RETURN

ENO

ENO

45
FUNCTION THERM(DD,T>

COMMON/CRIT/EM

O=DD*EM/iQQ0

.

THERsDlLTCTt FOCT(0,T) *10 0 .*EXCEST (D, T) +CRITC(O t T> THERM=THER

RETURN

END

FUNCTION CRITCiO.T)

COMMON/CRIT/ EH, EOK, RH, TC» DC, X , PC SIG COIN G/CM3 » T IN K

C 3 JTPUT UNITS ARE HM/H.K

AV«6.0225E*23 $ BK-i .3805<»E-16

C CALCULATE DISTANCE PARAMETER

R= (RN*»2.5>* (0**0.51 * ( AY/EMI **0 .5 R=R*(EQK**Q.5)*X/(T**0.5)

C GENERAL EQUATION

OX=0»100Q. O/EM

C OF»OT IN ATS PER OEG.

CALL OPDT(0PT ? OX,T)

DPT*OPT*i.01325E*6

C 3P0T NOW IN OYNES PRR OEG

C QPDO UN ATS » HOL/L

CALL DPDD(OPO,DX,T)

OPDsDPO*l. 01 325E*6*( EN/10 00. > C0PQ3 NOW IN OYNES, GM/CM3

C VISCOSITY IN GM/CM.S

VISsVISC(D,T)*(1.0E-G6)

COMPRES=1.0/ (O*OPOI**0.5

EX=BK*T**2MQPT**2) *COM?RES

EXB«R*( (BK*T)**Q .5)» (0**0.51 *((AV/ EM) **0.5»

CRIT*EX/(EXB*6.0*3.1<*159*VIS)

C PUT IN DAMPING FACTOR

BDQ= ( C Q-OC) / DC) **k

8TT= ( (T-TC) /TC)**2

FACT=EXP (-18.66*BTT - i».25*BDO>

DELC=CRIT*FACT/1Q0.0

CRITC=DELC

RETURN

END

FUNCTION S3UNO(0,T)

COMMON/CRIT/W

CALL OPOO(DP,D,T)

SOUNO=( (CP(0 ,T) /CV(0,T) ) *OP*101325 . /W) **.5

RETURN

END

FUNCTION PHI (D, T)

CALL DPDT(DT,D,T>

CSUBV=CV(3,T)

PHI=DT/(CSUBV*D)

RETURN

END

FUNCTION THETA(0,T>

CALL 0PDT(DT,D,T)

CALL DPDO(DO,D,T)

CSUBP=CP(3,T>

THETA=D*CSUBP*DQ/DT

RETURN

END

FUNCTION ENTHAL(P,D,T)

R=. 0820561a

00 = 0 TT = T CALL DSON(SO,00,TT)

CALL 0UDN(U0,0D, TT)

00=0

CALL DSDN(SG,OD,TT)

CALL UUON( JO ,00, TT)

ENTHAL=T* (5D-SQ) *13l . 3 25+( UD-UO ) * 101 . 325 *HI (T) (P/D-R*T) *101. 325

RETURN

END

SUBROUTINE RH01 ( PP , 3 D» TT

)

1ST CUT AT RHO FROM P=A«-B*T

REALLY AN ITERATION, BUT IT MAY BE SMALL ENOUGH AND FAST ENOUGH

P IN ATM, T IN K , RHO IN MOLES/LITER

DIMENSION RH0(h3) ,AC*3) ,B(<f3)

C

c
DATA ( ( RHO( I ) ,1=1,43) =1 . , 2 . .3 . ,4 « , 5 . ,6 . , 7 . ,8 t , 9. , 10 , . 11. , 12 , , 13 . ,14

1. 1 15.. 16. ,17. ,18.. 19. 1 20. .21.. 22. ,23. ,24. .25.. 26. ,27. ,26. ,29. ,30.

,

231. ,32. ,33. ,3*., 35., 36., 37., 3ft. .39. ,40. ,41. ,42. ,43.)

DATA ( (A (I) ,1 = 1, 43)«-. 295990 25741, -1.1388466175,-2. 4651376102 ,-4. 30

141059317, -6. 56965539 3. -9. 2 559185405, -12. 336266965.-15. 783961531,-1 29.573370659,-23.6810^1295.-28.086^58255,-32.77255*951,-3^7^

39,-42.942371352,-48.424432541,-54.194465142.-60.287558106,-66.7267

484771,-73.516518375,-80.643765419,-8 8.079905894,-95.781501825,-103

5. 68960 091, -111. 7277 f 842, -119. 79974172, -127. 786579, -135. 94257761.-1

64 3. 83339382, -151. 265 49964, -157. 99608468, -163. 718*857. -167. 27 870589

7, -169. 719223u 8, -169. 76998914, -166. 80 726845, -160. 22976454, -149. 6081

8063, -134. 683 0 0216, -116. 11756401, -94. 290117568, -67 .624247269, -36. 74

92895824,-1.8295576977)

DATA I (B(I) ,1=1, 431 =.0850 8 9208549, .175751585 96*. 2 7172674979, .372976

118724,. 479529114 65,. 59141 276615,. 70864366952, .83124664142,. 9592 644

29114, 1.0928 80265 8,1. 2322367731 ,1.3776455689, 1.5294993198, 1.6663212

3266, 1.654824 9918, 2. 0 30 00 5^26.2. 21510 0610 9, 2. 4112674705, 2. 619506666

W5, 2. 8<»070i»3212, 3.0 75 6<»50 876, 3. 325048554, 3. 5695525006. 3. 669S693541,

54.1658593689,4.4782981944,4.8204022069,5.1840092761,5.569646053,5.

6977252W362, 6. 4G 57175531, 6. 82142<*2299, 7.2760 697825,7. 7404616534. 6.1

7972575502, 8. 6358460 465 , 9. 0 4 732860 76, 9.426095360 7, 9. 6132635045, 10.

2

822281363, 10 . 63012981* ,1 1.0591 3828, 11. 565307194)

10 P*PP

T = TT

IP=1

PLO=0.

C CHECK INPUT DATA FOR RANGE

C INPUT RESTRICTIONS REMOVED

C IFCP.GT.350.) GO TO 15

C IF(T.GT. 300.1) GO TO 15

C GO TO 18

C 15 PRINT 16

C 16 FORMAT (24H INPUT DATA OUT OF RANGE)

C 00=0.0

C RETURN

C CHECK MELTING LINE

18 IF(T.GT.43.600) GO TO 23

PS=.0695-MT-13.8Q3>»30 .3312»EXP ( -5 . 693/T) « ( T-l 3 . 803) *2 • 0 »T/3 .

0

1 +0,000001

IF(P.LT.PS) GO TO 23

PRINT 22

22 FORMAT (26H INPUT CONDITIONS IN SOLID)

00=39.5

RETURN

C PHASEF INDER

23 IF(T.LT.32.938) GO TO 27

PPHASE=-50.6G02+1.92 0 8 8B*T

IFCP.ST.PPHASE) IP=15

GO TO 36

27 PVAP=VPN(T)

IF(P-PVAP) 36,32,35

32 PRINT 33

33 FORMAT 1 53H INPUT PLACES YOU EXACTLY ON THE VAPOR PRESSURE CURVE)

DO=DSATL(T)

RETURN

35 OENL=DSATL(T)

IP=OENL

C START TABLE LOOKUP HERE

3b CONTINUE

DO 40 I=IP,43

PCALC=A(I) *B (I)

M

IF(P.LT.PCALC) GO TO 41

PLO=PCALC

40 CONTINUE

IF(P.LT.PCftLC) 1=43

IF(P.LT.PCALC) GO TO 41

C PRINT 47

C 47 FORMAT < 35H HIGH OEMS IT Y» OUT OF RANGE FOR NOW I DD=44.

RETURN

41 CONTINUE

00 54 J=l,10

RHOF=J

RHOF=RHOF/10 . IF(I.EQ.l) GO TO 50

AA=A(I-1)MA(I)-A(I-1) )*RHOF

8B=8(I-1)*(B(I»-B(I-1) )*RHOF

GO TO 52

53 AA=A(I)»RHOF

BB=B(I) RHOF

5> PCALC=AA+3B»T

IF (P.LT.PCALC) GO TO 55

54 PLO=PCALC

55 FRAC=<P-PL0) / (PCALC-PLO) ^10.

00 =1-1

DD =00 +RHOF-0 .1+FRAC

RETURN

47
ENTRY Ti

FIRST GUESS FOR TEMPERATURE ITERATION OF FINOT

P=PP

D-DO

00 60 J*3 IF (D.LT.RHO CD) GO TO 61

50 CONTINUE

I=<*3

51 FRAC=D»RH0(I-1)

IF(I.EQ.l) GO TO 63

AA=A(I-1)+ ( A (I) -A (I-i) )»FRAC

BB=B(I-l)*(B(I)-B(I-i) )»FRAC

GO TO 62

S3 FRAC=0

AA=A(X)*FRAC

BBs8<Il»FRAC

bZ TT= CP~AA)/B3

FUNCTION ATKINT (X,Y MAT , XMAT* NELMTS, NMAX, NESSY • A CRCY)

C THIS PROGRAM HAS BEEN CHANGED SO THAT THE OSCILLATING NATURE OF

C THE MATRIX TO BE INTERPOLATED EXISTS ONLY AT THE UPPER END OF THE

C TABLE

C THIS ROUTINE WILL TAKE INPUT MATRICES OF UP TO 999 ELEMENTS EACH,

C ARRANGED SO THAT THE X MATRIX(XMAT) IS IN EITHER ASCENDING OR

C DESCENDING ORDER »SELEGT NMAX OF THESE POINTS » CHOSEN §0 THAT

C SUCESSIVE X VALUES OSCILATE ABOUT THE VALUE OF THE ARGUMENT X C UNLESS THE ENDS OF THE XMATRIX INTERFERE (IN THIS CASE THE

C OSCILATORY NATURE IS LOST BUT THE PROGRAM HILL STILL PERFORM AN

C INTERPOLATION) * INTERPOLATE ON THESE NMAX PAIRS OF DATA BY

C AN OSCILATING VARIABLE POINT AITKEN INTERPOLATION ALGORITHM

C EITHER UNTIL THE PERCENTAGE CHANGE IN THE INTERPOLANT IS LESS

C THAN THE ACRCY ARGUMENTUHE ARGUMENT NESSY INDICATES THE

C NUMBER OF THE POINT JUST BEFORE THE LAST ONE CHECKED) OR UNTIL

C THE NMAX POINTS ARE ALL USED. IT IS SUGGESTEO THAT NMAX

C BE LESS THAN 10, AND OF COURSE LESS THAN NELMTS. NELMTS

C INOICATES THE NUMBER OF ELEMENTS IN XMAT OR YMAT

C IF NESSY IS ZERO IT INOICATES THAT THE INTERPOLATION REQUIREMENT

C HAS NOT SEEN SATISFIED* IF NESSY IS 1 IT MEANS THAT THE VALUE OF

C X LIES OUT SIDE THE RANGE OF XMAT.

DIMENSION YMAT(999>, XMAT ( 999) ,A ( 21 , 20

)

103 FORMAT (^^INTERPOLATION REQUIREMENT NOT SATISFI ED( Xe»E16 . 8 ,1H) /33H

1 L AST 2 APPROXIMATIONS OF Y ARE ( Y= , E16 .8 , 1H, , E16 . 8 , 1H)

)

203 FORMAT ( 55HT H IS REPRESENTS AN EXTRAPOLATION OF THE XMAT MATRIX(X*,

1E16.8.1H) /33HN0 CALCULATION HAS BEEN PERFORMEO)

333 FORMAT (2*tH NELMTS IS LESS THAN NMAX)

<*03 FORMAT (22HNMAX IS LARGER THAN 20)

IF(NMAX-20) 71,71,69

59 PRINT t»00 ATKINT=0.0

RETURN

71 IF(NMAX-NELMTS)75,75,73

73 PRINT 300

ATKINT=0.0

RETURN

75 CONTINUE

C FIRST TWO SUCCESSIVE VALUES OF THE XMATRIX THAT STRADDLE THE

C VALUE X WILL BE SOUGHT

JJ1=NELMTS-1

00 20 I=1.JJ1

DIF1=X-XMAT (I)

DIF2=XMAT(I*1)-X

IF(0IF1)16, 15,16

15 ATKINT=YMAT (I)

NESSY =NMAX

RETURN

la I F { OIF 2 ) 18 , 1 7 , 18

17 ATKINT=YMAT (I+l)

NESSY =NMAX

RETURN

13 RATI0=DIF1/ OIF2

IF (RATIO)20 ,20,19

19 IMIO=I

GO TO 32

20 CONTINUE

C AT THIS POINT ONE COULD PRINT THE FOLLOWING STATEMENT

C WRITE OUTPUT TAPE 6,200,X

NESSY=1

ATKINT=0.0

RETURN

32 CONTINUE

C NOTE THAT RATIO IS POSITIVE IF THE TWO POINTS STRADDLE X C REGARDLESS WHICH IS LARGER

JJJ=IMI0

JUP=IMIO

JON=IMID

48
IF< JJJ+NMAX-NELNTS*1>98,98,102

98 DO 201 J=1,NMAX

JJJ=IMIO*J-l

A(1,J)«XMAT ( JJJ)

201 A<2,J)«YMAT (JJJ)

GO TO 203

102 00 kl J»1,NMAX

JJ*J/2

J0E*J-2*JJ

C JOE IS 0 IF J IS EVEN AMD 1 IF J IS 000

IFCJ-1133,%0,33

33 IF(JDN-l)3i»,36,3i»

Zk IFCJUP-NELMTS)35,37,35

35 IF ( JOE) 37,36,37

35 JUP=JUP*1

JJJ=JUP

GO TO kQ

37 J0N«J0N-1

JJJ=JDN

GO TO kO

•*Q A(1,J)«XNAT(JJJ)

A (2,J) «YMAT ( JJJ)

til CONTINUE

203 NNN«NMAX*1

00 6 J*3,NNN

L = J-1

00 5 K»L,NHAX

C J IS THE COLUMN NUMBER

C K IS THE ROW NUMBER

OA (J,K)« (ACJ-l.K) -A (J-l,J-2) )»(X-A(1, J-2) )

/

I A (1 , K ) -A CI, J-2 1

)

1 *A(J-l,J-2)

IF(K-L) 3,2,3

2 IF (ABS (<A(J,l)-A(J-l,L-l) ) / A < J,l ) ) - ACRC Y/10 0 . Q ) 7,7,3

3 CONTINUE

5 CONTINUE

5 CONTINUE

NESSY=0

C AT THIS POINT ONE COULO PRINT OUT THE FOLLONING STATEMENT.

C WRITE OUTPUT TAPE 6, 10 0 , X , A ( NNN,NMAX) , A ( NNN-1 ,NMAX-1)

ATKINT-A (NNN ,NMAX)

RETURN

7 NESSY=J-1

ATKINT = A < J » L)

RETURN

ENO

49
Appendix B. Test Program and Sample Results

PROGRAM VALUES

C A SAMPLE PROGRAM TO CHESK RUNNING AT OTHER INSTALLATIONS.

C CALCULATES T HERMOFUNCT IONS OF PARA-H2 FROM THE 32 TERM MBM*.

C INPUT IS RE AO FROM CAROSf P IN ATM, T IN DEG K. C OUTPUT * HEADING!

C A LIMITED

C REQUIRED

i

C THE PROPERTIES DECK MAS LAST REVISEO ON 75/06/17

t HEADING. UNITS. ANO VALUES ARE PRINTED. HOWEVER. ONLY

ED NUMBER OF THE VARIABLES ARE CHECKED BY THIS SAMPLE OECK.

:0» ALL SUBROUTINES / FUNCTIONS LISTED IN NBSIR 75-014.

VEL/82H

J/MOL-K

RHO

ATM

J/MOL-K

CALL 0ATAPH2

PRINT 20

20 FORMAT I82H1 P

1 C-V C-P

2 J/MOL J/MOL-K

DO 19 I = l,<»

k READ 5,Pjf

5 FORMAT<F6.3,F7.3>

OEN=FINDD(P,T»

H =ENTHAL ( P » OEN *T)

S =ENTROP(OEN,T)

CVE*QV f DENtTI

CPE=CP<DEN,T)

VEL*SOUNO(DEN.T)

12 PRINT 13.P,T,0EN,H,3,CVE»CPE.VEL

13 FORMATIFi0.'*,F12.5,Fi0.6,F12.1 t 3F10.2,F8.0l

19 CONTINUE

END

KELVIN

M/SEC/1H ) MOL/L

P T RHO H S C-V C-P VEL

ATM KELVIN MOL/L J/MOL J/MOL-K J/MOL-K J/MOL-K M/SEC

1.0000 20.00000 35.279160 -521.9 15.84 11.33 19.11 1111

1.0000 30.00000 0.420408 602.5 69.33 12.57 21.79 447

15.0000 34. 00000 17.42'»046 63.1 3<*.59 16.10 301.40 425

70.0000 25.00000 36 .713628 -289.9 17.53 12.30 19.52 1306

Appendix C. Conversion Factors

Temperature 1.8 R = 1 K

Pres sure 14.695949 psia = 1 atm = 1.01325 x 10

5 N/m2

Specific Volume 0.00794590 ft3 /lb =

3

1 cm /mol m 3

= 0.001 liter = 10- 6 m

3

( 1 cm ) Internal Energy, Enthalpy

Entropy, Specific Heat

Thermal Conductivity

Viscosity

Speed of Sound

Molecular Weight

Surface Tension

0. 213405 BTU/lb = 1 J/mol m

0. 118558 BTU/lb R = 1 J/mol-K m

0.0578176 BTU/ft-hr-R = 1 mW/cm-K

0.067196897 lb /ft-s = 1 g/cm-s m

3. 2808 ft/s = 1 m/s

2.01594""

0. 5710147 x 10" 5 lb

f /in = 1 dyn/cm

(1 dyn = 10" 5 N)

12 On the C = 12.000 scale

50
NBS-114A (REV. 7-73)

U.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET

1. PUBLICATION OR REPORT NO. NBSIR 75-814

2. Gov't Accession

No.

3. Recipient's Accession No.

4. TITLE AND SUBTITLE A MODIFIED BENEDICT- WEBB-RUBIN EQUATION OF

STATE FOR PARAHYDROGEN - II

5. Publication Date June 1975

6. Performing Organization Code 275. 02

7. AUTHOR(S)

Hans M. Roder and Robert D. McCarty

8. Performing Organ. Report No,

9. PERFORMING ORGANIZATION NAME AND ADDRESS NATIONAL BUREAU OF STANDARDS DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

10. Project/Task/Work Unit No. 2750423

11. Contract/Grant No.

12. Sponsoring Organization Name and Complete Address (Street, City, State, ZIP)

1 NASA- Johnson Space Center

Houston, Texas 77058

13. Type of Report & Period Covered Final

14. Sponsoring Agency Code

15. SUPPLEMENTARY NOTES

16. ABSTRACT (A 200-word or less (actual summary of most significant information, ((document includes a significant bibliography or literature survey, mention it here.) A 32 term modified Benedict-Webb-Rubin equation of state has been applied to data

for parahydrogen. The adjustable parameters in the equation of state were determined

using 2665 points including very recent measurements at lew temperatures and high

pressures. The new values extend the range of the PVT data sufficiently to warrant

a refitting of the equation of state. Temperatures for the data range from the triple

point to about 700 K with pressures reaching 3000 atmospheres near ambient temperatures. The PVT data were adjusted to the T^g scale. In addition, extensive

modifications have been made to the previously accepted PVT surface in the region near the critical point. These adjustments have been made on the basis of more

recent refractive index data and the application of scaling law equations. Detailed

comparisons between experimental and calculated values are given for density.

Corresponding comparisons are made for enthalpy and the specific heat at constant

pressure.

17. KEY WORDS (six to twelve entries; alphabetical order; capitalize only the first letter of the first key word unless a proper name; separated by semicolons) Critical point; density; enthalpy; equation of state; hydrogen; index of refraction PVT;

saturation properties; scaling laws; specific heat

18. AVAILABILITY [X Unlimited 19. SECURITY CLASS

(THIS REPORT)

21. NO. OF PAGES

1

' For Official Distribution. Do Not Release to NTIS UNCL ASSIFIED 55

1

1 Order From Sup. of Doc, U.S. Government Printing Office Washington. D.C. 20402. SD Cat. No. C13

20. SECURITY CLASS

(THIS PAGE)

22. Price

1 X! Order From National Technical Information Service (NTIS)

Springfield, Virginia 22151 UNCLASSIFIED $4.25