A MODIFIED BENEDICT-WEBB-RUBIN EQUATION

OF STATE FOR GASEOUS AND LIQUID OXYGEN

L A. Weber

'-mi

Cryogenics Division

Institute for Basic Standards

National Bureau of Standards

Boulder,Colorado 80303

April 1978

Prepared for

NASA-Johnson Space Center

Houston, Texas 77058

NBSIR 78-882

A MODIFIED BENEDICT-WEBB-RUBIN EQUATION

OF STATE FOR GASEOUS AND LIQUID OXYGEN

L. A. Weber

Cryogenics Division

Institute for Basic Standards

National Bureau of Standards

Boulder, Colorado 80303

April 1978

Prepared for

NASA-Johnson Space Center

Houston, Texas 77058

U.S. DEPARTMENT OF COMMERCE, Juanita M. Kreps, Secretary

Sidney Harman, Under Secretary

Jordan J. Baruch, Assistant Secretary for Science and Technology

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director

Contents

Page

1. Introduction

.

1

2. The Data 2

3. The Equation of State for Oxygen 4

4. Auxiliary Information 6

4.1. Two-Phase Boundaries 6

4.2. Vapor Pressure Curve 6

4.3. Saturation Densities 7 4.4. Melting Pressures 8

4.5. Ideal Gas Thermodynamic Properties 8

4.6. Summary of Fixed Points for Oxygen 8

5. The Properties Deck 10

6. Quality of the Fit and Accuracy of the Derived Properties .... 10

6.1. Density Deviations 10

6.2. Comparison of Derived Thermodynamic Properties 11

7. Summary 12

8. References 31

Appendix A 13

Appendix B 26

List of Figures

Figure 1. Deviations of the density data from eqn. (1) on several

isotherms 27

Figure 2. Deviations of the saturation densities from the values

calculated via eqn. (1); 0 - liquid, 0 - vapor 28

Figure 3. Deviations of the heat capacity data on three isochores

from values calculated by eqn. (1) 29

Figure 4. Deviations of the velocity of sound data from values

calculated via the equation of state on several isotherms . . 30

I

iii
List of Tables

Table 1. Data Used to Determine the Parameters of the Equation

of State

Table 2. Coefficients for the Equation of State (1)

Table 3. Parameter Values for the Various Equations

al.:!i;;:=.i'':i ''o

, , . >'-!•?.<' - 0 , 51; f t p.; ...

:1^;U'aKV :\j;'.'o

iv
A MODIFIED BENEDICT-WEBB-RUBIN EQUATION OF STATE

FOR GASEOUS AND LIQUID OXYGEN

L. A. Weber

Cryogenics Division

Institute for Basic Standards

National Bureau of Standards

Boulder, Colorado 80302

New PVT data on oxygen at pressures to 800 atmospheres have

been recently reported. These data, along with other NBS data on oxygen have been fitted with a 32 term modified Benedict-Webb-Rubin

equation of state. This equation is valid for gaseous and liquid

oxygen at temperatures from the triple point (54.36 K) to 300 K

and for pressures up to 800 atm. Thermophysical properties

calculated via this equation are compared with the most reliable < values available. The vapor pressure equation has also been

updated to be consistent with more recent experimental work near

the triple point. A computer program is available which calculates

the thermophysical properties of oxygen.

Key words: Density; enthalpy; entropy; equation of state; oxygen; oxygen

property computer program; PVT; saturation properties; specific

heat.

1. Introduction

A modified form of the Benedict-Webb-Rubin (BWR) equation of state [1],

with more terms added, is often the most convenient way of representing the

PVT surface of a fluid and of calculating its thermodynamic properties. Recent

versions, such as the thirty-two term modification of Jacobsen [2], also fit

PVT data with deviations no larger than the experimental uncertainty over most

of the PVT surface.

The fact that many of the derived thermodynamic properties can be calculated

explicitly in closed, analytic form greatly enhances the usefulness of the equation
and the speed of machine calculations. This characteristic of the equation also

allows it to be fit to several properties simultaneously (e.g., density and

specific heat). Such multi-property fitting virtually assures the accuracy of

most derived thermodynamic properties, provided sufficient good data are

available.

In an earlier report, Roder and McCarty [10] described the application of

the 32 term version of the equation to the data on parahydrogen. This report

describes the fitting of the equation to the NBS data on oxygen. Stewart and

Jacobsen [2,6] have successfully applied this equation to oxygen. However, the

availability of significant new PVT data [4] justifies the refitting described

here. ;^.>--ll<:^:J. .^w . ^. - '- -

We have changed slightly the vapor pressure equation used previously, and

the equation with the new parameters is repeated here. Also repeated here are

the representations used previously for the densities of saturated liquid and

vapor and for the melting pressure as functions of temperature.

The Data

All of the data used for fitting the equation came from this laboratory.

These comprise the bulk of all the data available for oxygen, and provide us with a data set having a high degree of internal consistency. The data used

comprise a close-spaced set of 2064 points, summarized in table 1. The

experimental details are given in refs [3,4,5]. With the exception of some critical region measurements from [11], the saturation boundary data are not

independent experimental data but were derived from the intersection of the

single-phase experimental isochores with the vapor pressure curve. The resulting

points were smoothed and interpolated by means of the function given in

section 4.3, and the smoothed data were used in the equation of state fitting.

The Gibbs constraints are not actually constraints but rather calculated data

setting the Gibbs free energy of the coexisting liquid and vapor states equal at

a given temperature.

The earlier PVT data of [3] were not reported on the currently accepted

International Practical Temperature Scale of 1968 (T-n), and thus it was

bo

necessary to convert them to this scale by means of the conversion tables given

in Metrologia [7].

2
Table 1. Data Used to Determine the Parameters of the

Equation of State.

Type of

o n jjata

No. Points Range Source

T(K) P(atm)

PVT 1435 DO - JUL) 0 "3/. n Weber [3]

PVT 348 56 - 300 50 - 800 Weber [4]

Saturation

Boundary 100 56 - 154 0.002 - 49 Weber [JJ

C

V

131 58 - 284 42 - 284 Goodwin & Weber [5]

Gibbs

Constraint 50 56 - 154 0.002 - 49 Weber [3]

3
3. The Equation of State for Oxygen ^ The original Benedict-Webb-Rubin equation of state [1] was modified by

Strobridge [8] in 1962. Strobridge's equation contained 16 terms, and it was

adequate for the data then available for nitrogen. However, when Roder and

Goodwin [9] applied it to a close-spaced, high precision data set for parahydrogen

they were forced to divide the PVT surface into liquid and gaseous regions and

use two sets of parameters in order to fit the experimental data. They were

unable to eliminate discrepancies in the derived properties at the junction of

the two surfaces. Since 1974 McCarty [10] has applied the 32 term version of

Jacobsen [2] to several cryogenic fluids. The resulting "single surface"

representations reproduced the data very well except in the critical region

where no analytic representation can be expected to show proper behavior. This

version was therefore chosen to represent the data for oxygen. The equation of

state is

P = pRT + (A,T + A„T'^''^ + A_ + A, /T + A^/T^)

1 2 3 4 5

+ p^ (A^T + A^ + Ag/T + Ag/T^)

+ p^ (A^qT + A^^ + k^^/l) + p^ (A^3)

+ p^ (A^^/T + A^^/T^) + p^ i^i^l'^)

+ p^ (A^^/T + A^g/T^) + p^ (A^g/T^)

2 p^ (A2q/T^ + A^^/T^) + p^ (A22/T^ + ^21>''^^^

+ p^ (A^^/T^ + A^^/T^)

+ P' (A26/t' + A^/tS

+ p^^ (A28/T^ + A^g/T^)

+ p^^ (A^q/T^ + A^^/T^ + A^^/T^)

4
Table 2. Coefficients for the Equation of State (1).

R = 0.08205616 L atm/mol K y = -0.0056

G(l) -.4308768468E-03

G(2)

G(3)

G(4)

G(5)

-.4143014968E+01

-.1270637452E+05

.1979591095E+00

.1853654396E+03

G(6) .1536388737E-04

G(7)

G(8) -.2199275123E+01

.1326068945E-02

G(9)

G(10)

.4705445127E+04

.4728198017E-06

G(ll) = .2430408198E-02

C(12) = -.1896759615E+00

G(13) = -.6887067207E-05

G(14) = -.6132885180E-03

0(15) = -.1836518694E+00

0(16) = .2575663871E-04

0(17) = -.2415604646E-06

0(18) = .1438680831E-03

0(19) = -.1703915986E-05

0(20) = -.2353705917E+04

0(21) = -.2271707669E+06

0(22) = -.2753815471E+02

0(23) = .9277648729E+05

0(24) = -.4114926856E-01

0(25) = .1982233262E+01

0(26) = -.1239651142E-03

G(27) = -.6322588664E+00

0(28) = -.2443207666E-07

0(29) = .1328704370E-04

0(30) = -.1146313812E-09

0(31) = -.1021169305E-07

0(32) ^ .2334998237E-06

5
The coefficients, given in table 2, were found by least Hquare fitting eq (1) to

the data set in table 1, with the density, f), in mols/ilter, pressure, 1', In

atmospheres, and the temperature, T, in Kelvins.

In the fitting procedure each data point was weighted according to a scheme

which utilized the estimated uncertainties in the experimental quantities. In

addition, three constraints were imposed on the equation. At the critical point,

the density was constrained to be equal to p ; also the derivatives (8P/9p) and

2 2 ^

(9 P/9p )^ were constrained to be equal to zero. The critical point used was the one determined previously [111: p = 13.63 mol/L, T = 154.581 K, and

c c P = 49.768 atm.

4. Auxiliary Information

' - ; •' 4.1. Two-Phase Boundaries

The equation of state representation of the previous section does not

recognize the existence of phase boundaries, and these must be specified independently to complete the description of the PVT surface. The liquid-vapor two

phase region can be inferred from van der Waals type "loops" in the isotherms of

the equation of state surface at subcritical temperatures. A liquid-vapor two

phase boundary (saturation boundary) can be calculated from the equation of state

surface by setting equal the Gibbs free energy of the liquid and vapor states at

a given temperature. The vapor pressure thus derived generally agrees with the

true vapor pressure to within 0.1 - 0.2 percent, liowever, the results are generally more accurate if we specify the vapor pressure curve independently.

. .V 4.2. Vapor Pressure Curve

Although a polynomial representation of the vapor pressure can be made

sufficiently accurate for most purposes, the formulation given by Goodwin [12]

has many advantages

:

Inl^ = Ax + Bx^ + Cx^ + Dxd - X)^ , ^ (2)

where X

= (1 " T^/T)/(l - T^/T^), and P^ and T^ are the triple point pressure

and temperature. This equation has the proper limiting form (£nP a + b/T) at

6
low pressures near the triple point. It has the proper non-analytic form to fit

data near the critical point, and it fits the data with fewer parameters than a polynomial. The parameters have similar values for all cryogenic fluids, and

the exponent, e, is close enough to the value 3/2 that it is sometimes considered

a "universal parameter" for this class of substances. If the triple point

pressure is known, eq (2) is constrained to that value; otherwise it may be used

as a variable parameter. The parameter values used in eq (2) are given in

table 3.

4.3. Saturation Densities

Given the vapor pressure curve, saturation densities can be calculated via

the equation of state (1). It is sometimes convenient, however, to have an explicit formulation of the densities of the saturated liquid and vapor. They

may be expressed by

Pg^^ (liquid) = Pj^ + f (T) (3a)

PSAT ("^P°^) = PRD - ^^^^

where the rectilinear diameter, p , (defined by p = \ [P^at (liquid) +

•^SAT (^^P°^)]' is

PRD = Pc ^l^^c - ^) ^ ^2^\ -

and

(5)

J=l

The parameter values are given in table 3. With the parameters given, eqs

(3-5) should only be used in the temperature range 120 K to T^ = 154.581 K.

7
' 4.4. Melting Pressures

The melting pressure data from [3] and [13] are accurately represented by

tlic equation

P = + [.(T/Tj.)'' - 1] (6)

up to a pressure of 3500 atm. Parameter values are given in table 3.

4.5. Ideal Gas Thermodynamic Properties

By means of the relationships given in [3], equation (1) can be used to

calculate changes in thermodynamic properties along any isothermal path, including

those passing through the liquid-vapor two phase region. In order to relate

property values at different temperatures, we must add to these changes the

values of the thermodynamic properties in the ideal gas state at that tem.perature.

The latter are easily calculated via the relationship used for the ideal gas

specific heat

C°(T) 2 u R

^

J ^

(e"" - 1)^

J=l ^ ^

where u = A^/T. Parameter values are given in table 3. From eqs (1) and (7)

and the values of the enthalpy and entropy in the ideal gas state at one

temperature, all of the thermodynamic properties of the fluid may be calculated

at any temperature and pressure. For the reference state we have chosen the

values at T = 298.15 K, P = 1 atm,

H° = 8682 J/mol

; S° = 205.037 J/mol-K

from CODATA [14].

4.6. Summary of Fixed Points for Oxygen

In this report the following values have been specified as fixed; the

subscripts t, b, and c refer to the triple, boiling, and critical points

respectively.

8
Table 3. Parameter Values for the Various Equations.

Vapor pressure, eq (2)

A = 7.797 7723

B = 4.577 3000

C = -1.928 1264

3.293 8576

1.5

Saturation densities, eqs (4) and (5)

-5

a^ = 6.0402 X 10

= 1.811 3127

A^ = 2.775 1793 x 10

-1

a^ = 9.80 X 10"^

A^ = -7.580 9408 x 10

6 = 0.353

-1

Melting pressure, eq (6):

P = 2637.3 atm

o 0.0014451 atm

1.769

54.359 K

Ideal gas specific heat, eq (7)

A^ = -1.86442361 x 10

A. 2.0784024 x 10

A^ = -3.42642911 x 10

A, = 3.50297163

4

A. = 2.05866482 x 10

-1

-7

A, = -1.11035799 X 10

6

A^ = 2.08612876 x 10

hr, = 1.01894691

o A^ = 2.23918105 x 10"

-11

9
Triple Point : = 54.359 K

• P^ = 0.0014451 atin

Boiling Point : = 90.188 K

P, = 1.0 atm

b

Critical Point: T = 154.581 K

c P = 49.768 atm

c p = 13.63 mol/L

c

5. The Properties Deck

In the appendix is a listing of a set of subroutines useful for utilizing

the foregoing relationships to calculate the thermophysical properties of

oxygen. This properties deck was developed by R. D. McCarty [15] for use on

the Control Data Corporation 6600 computer, and it is available from the NBS

Cryogenic Data Center, Boulder, Colorado. A listing of the deck is given in

Appendix A, and a test program and sample results are in Appendix B.

6. Quality of the Fit and Accuracy of the Derived Properties

Because of the wealth of consistent data available for oxygen, this report

provides a good indication of the ability of the 32 term modified BWR equation

to fit the PVT surface of a classical gas and to provide accurate derived

properties. The following sections provide enough comparisons with experimental

data to show the general quality of the fit and to indicate regions of difficulty.

6.1. Density Deviations

In fig 1 are shown deviations of the density data from the calculated

surface at eight temperatures, 70, 120, 140, 150, 156, 170, 200, and 300 K.

Overall we see that the equation reproduces the data very well. This is especially

true in the compressed liquid. Only in the region of the critical point does the

equation fail to reproduce the PVT surface accurately. The deviation plot for

the 156 K isotherm, and to a lesser extent the 170 K isotherm, are typical of

the Inability of an analytic equation of state to fit PVT data in the critical

region. Neglecting those points, however, the average deviation of the equation

from the data is about 0.04 percent in density. This amount is about twice the

imprecision in the data, but it is less than the overall uncertainty in the data

of about 0.1%, due to possible systematic errors.

10
Deviations of the calculated saturation densities are given In fl;; 2. Here

again it is seen that the deviations are larger in tiie critical region. They

become smaller in the immediate vicinity of the critical point, however, because

the equation is constrained at that point.

6.2. Comparison of Derived Thermodynamic Properties

(a) Specific heat at constant volume : Goodwin and Weber [5] measured

for oxygen at eighteen different densities. One, apparently spurious, set of

data at a density of 5.1 mol/L was discarded, as were several data near critical,

and the rest were incorporated in the data set used for determining the equation

of state. These data had a reproducibility of about 0.25 percent. The average

difference between the calculated and experimental values was 0.67 percent,

approximately equal to the overall uncertainty of the data as expressed in [5].

The maximum deviation for the data used in the fit, 4.2 percent, occurred at 162 K

near the critical density. Calculated values deviated from several other near critical experimental points, not used in the fit, by as much as 20 percent. This

behavior is typical for an analytic equation of state near the critical point.

Deviations of the calculated C 's from the experimental data at three densities

V are shown in fig 3.

(b) Velocity of sound . Extensive measurements of the velocity of sound in

oxygen have been reported by van Itterbeek and van Dael [16] and by Straty and

Younglove [17]. Experimental pressures ranged as high as 900 atm in the

compressed liquid and 350 atm for the high temperature gas. These data have an experimental uncertainty of 0.1 percent or less. Therefore, although they

cannot be readily incorporated into the multiproperty fitting techniques, they

do provide a sensitive test of the quality of the calculated PVT surface.

The data of van Itterbeek and van Dael are on five isotherms from 67.5 to

90.4 K at pressures to 900 atm. The calculated velocities are generally smaller

than the experimental ones by a little over one percent, with the maximum

deviation being 2.0 percent. Comparison with the data of Straty and Younglove

indicates that the calculated values are again smaller, the difference averaging

about 1.4 percent in the compressed liquid and decreasing to about 0.4 percent at

supercritical temperatures. The deviations on several isotherms are shown in

fig 4.

11
(c) Other thermodynamic properties . Other properties such as entropy,

enthalpy, etc. could be used In the multi-property fitting process. However, no experimental measurements exist for these properties of oxygen. Instead we have

compared values calculated via the BWR equation with values calculated in ref [4],

which used a different type of surface representation.

, A spot check of the calculated entropies at 60, 100, 200 and 300 K and

pressures from 50 to 1000 atm showed differences ranging from 0.03 to 0.12 J/mol-K.

This agreement is quite good and is within the estimates of accuracy in [4].

The maximum differences might be expected to occur at the critical point, and

here the two methods disagree by only 0.3 J/mol-K. A similar comparison for

enthalpy showed agreement generally within + 7 J/mol.

'f"^ At the liquid-vapor two phase boundary, the Gibbs "constraint," which sets

G - G„ = AG = AH - TAS = 0 , ^ . , , ,

was satisfied to within 1-3 J/mol at all temperatures from 56 to 154 K.

In the critical region, which for present purposes we may define as + 0.25

at the critical temperature and extending out 20 K above T^ at the critical

density, the isothermal compressibility, (9p/9P)^/p, is subject to large errors.

For this reason values of the specific heat C should be considered unreliable

P

in this region.

7. Summary

A total of 2064 data points for oxygen, including PVT, C^, and Gibbs free

energy data, have been fit with a 32 term modified Benedict-Webb-Rubin equation

of state. The equation represents the data very well in all regions except

near the critical point. The equation is valid over the whole range of the data

from the triple point (54.36 K) to 300 K and to pressures of 1000 atm, but it

should not be extrapolated beyond this range without caution. It allows

calculation of the derived thermodynamic properties with an accuracy equal to

that with which they are known, with the possible exception of the sound

velocities at high densities. Even the latter are calculated with an accuracy

sufficient for most purposes.

12
Appendix A.

SUBROUTINE DATA 02

DIMENSION GI32J , VP(9) ,GI<11»

DIMENSION GV{9»,GT(9),FV<'*),FT('+),EV(8),ET(8) C0MM0N/0ATA1/GV,GT,FV»FT , EV,ET

COMMON/SEN/BETA, XO. DELTA, El, E2, AGAM

COMMON/CRIT/ EM, EOK, RM, TC, DC, X , COMMON/0 AT A/G, R, GAMMA, VP, OTP

COMMON/ISP/N,NW,NWW

COMMON/CPIO/GI

NWW = 0

N=0 $ NW=1

X0=0.183 $ BETA=0.355

EM=31.998S $ EOK=113.0

OC = 0,t»362 5 X=2. 210636

AGAM = 1.190

TRANSPORT PROPERTIES

FROM HANLEY £T. AL.

VOL 3, NO l<^7k

GV<1) =-9.7076378 593E+1

GV(3)=-2.'»668758803E*1 $

GV(5) = 3.78510<»9522E-1 $

GV (7) = 1.113'*<**»130'tE-2 $

GV(9)=i. 0 27 93796<*lE-5

GT <l)=-2. 0395052193E + 5 $

GT (3)=-1.201^175183E + 5 | GT <5) = -5.i»2'»*»239598E + 3 $

GT (7)=-3.285<*821539E + l $

GT <9) =-1.^610986 82 OE-2

EV(1) = -1.2152387017E + 1 $ EV (2) =-3. l^»2172899tfE-l

EV (3) = 1. 8201161i*68E + l $ E V ( < ) =2 . 7390'*29525E + 2

EV (5) = -2.7i»989569t»8E+3 $ EV (6) =2. '3'*0689667E + 2

EV(7)=i.l91150ttlO'»E + 2 S EV(8) = .<»35

FVtDs**. 3526515153 $ FV(2)= -2.0361263878

FV(3)=1.'* $ FV('»)=100.

FT(1)=.3060 $FT(2)=.2785

FT(3)=1.12 $ FT(tt)=100.

ET (l) = -1.2310«t00765E + l I ET I 2) = 1 . 679950i*261E +

1

ET(3)=-2.99'»'t87872lE+3 I ET (<) =i*.7350508788E + 2

THERMODYNAMIC PROPERTIES

FROM WEBER - TO BE PUBLISHED

R=8.205616E-2

GAMMA=-. 0056

PC, SIG

$ DELTA = i*. 352 $ El = 2.21 $ £2=0.287

$ RM=3.8896E-08$ TC = 15i*.575

$ PC=i»9. 770 $ SIG = 3.'*37

$ GV(2) =8. 280125i»201E + l

GV('t) = 2.132<*3602^+3

GV(6» =-1. 0«t8 7216090 E-1

GV( 8) =-5.36760 93 757E-tt

GTI2J =2.'*0881ttl709E*5

GT (i*) =3.295if9'*919E + it GT(6) =5.^73^8655^0E2

GT{8)=1. 0753572103

GI(1)= -0.'*981998537119^i3E

GI<2)= 0.2302^7779995218E

GI(3)= -0.3'»5565323510732E

GI(^t)= 0. 352187677367116E

GI(5)= -0,i*35't202l602if'f20E

GI(6)= 0.13«»6353W5013162E-07 GH7)= 0.162059825959105E-10 GI{8)= 0.1031i»6851572565E 01

GI(9)= 0.223918105000000E Qk

GI (10>=GI (11» =0

T0=298.15 $ H0=8682. $

GI (10)=HO-HI(TO)

GI (11) -SO-SI (T0»

G(l)= -.'308768^68E-03

G(2)= .1979591095E*00

G(3)= -.t»li*301'#968E + 01

G(*»)= .185365it396E4-03

GC5)= -.1270637'»52E+05

GI6)= .1536388737E-0^*

G(7)= .13260689'*5E-G2

G(8)= -.2199275123E+01

G(9)= .'705<»**5127E+0**

G(10)= .'728198ai7E-06

G(ll)= .2'*30i»08i98E-Q2

G(12)= -.1896759615E+00

G«13>= -.6887067207E-05

Gllk)- -.6132885180E-03

G(15)= -.183651869'+£*00

G(16)= .2575663871E-0't

G(17)= -•2<»1560<+6it6E-06

GI18)= .lit38680831E-03

G(19)= -.1703915986E-05

G(20)= -.2353705917E + 0^t G(21)= -.2271707669E+06

G<22)= -.2753815'*71E*02

G(23)= .92776i»S729EvC5

G(2«»)= -.<*11'*926856£-01

Oi» 03

01

01

0^*

SO=205.037
G(25)= •1982233262E+01

GI26J= -•12396511<*2E-03

GI27)= -•632258866i*E400

GI28)= -.2'»'»3207666E-07

G{29)= •132870'»370E-0i»

G(30)= -.ll<»63138l2E-09

G(31)= -.1021169305E-07

GI32)= .233if998237E-06

VP (1) = 7. 7977723 $ VP ( 2 ) 57730 $ VP( 3) =-1. 9281261*

VP<<»)=3. 2931232 $ VP<5) = 1.5 $ VPC6) =0. 001«»U51

VP(7)=5U,359 $ VP(8)=15i*.581 VP{9)=0.0

DTP=*f0.8<*

RETURN $ END

FUNCTION VPNITT)

CALCULATES \/APOR PRESSURE IN ATMOSPHERES

FOR AN INPUT TEMPERATURE IN KELVIN

DIMENSION G{32l»VP(9) COMMON/0 AT A/G,R, GAMMA, VP

T = TT

X= (l.-VP (7) /T) /<1,-VP<7) /VP(8)

)

VPN=VPI6) 'EXP IVPU) X+ VPf 2) X'XfVP (3 J X»»3*VP (9) *X»»i* + VP C if ) X*

1 (1,-x) *vp(5n

RETURN

END

FUNCTION OPOTVPITT)

CALCULATES THE DERIVATIVE OF PRESSURE

WITH RESPECT TO TEMPERATURE AT

SATURATION. INPUT IS TEMP. IN K, OUTPUT IS ATM/K.

COMMON/DATA/G.R. GAMMA, VP

DIMENSION G(32)fVP(9)

T = TT

IF(TT.GT. VP(8) IGO TO 1

X= (l.-VP (71 /T) /<1. -VP (7) /VP<8)

)

OXDT= ( VPI75/T»»2)/ C1.-VPC7) /VPI8)

)

DPOT=VP C1»»DXDT4-2.*VP(2) »X*0XDT + VP(3) 3.*X»*2*DXDT + VP(tf) K (l.-X) *»VP<5) )*OXDT + VP(t*)»X* <(1.-X)**(VP<5)-1.))»VP<5)»(-DX0T)

OPOT=OPOT»VPN<T)

DPOTVP=DPDT

RETURN

DPOTVP=0

RETURN

END

FUNCTION FIN0TVIP08S)

ITTERATES THE VAPOR PRESSURE EQUATION

FOR A TEMPERATURE ( IN KELVIN)

GIVEN AN INPUT PRESSURE IN ATMOSPHERES

COMMON/DATA/GiRf GAMMA,VP,OTP

DIMENSION G{32I,VP(9)

T= (VPC7) VP! 8» 1/2.

00 7 1=1,10

P=VPNI T)

IF (A8S « P-POBS J -.000 0 01*POBS} 8,8,6

6 CONTINUE

CORR=(POBS-P) /OPDTVP(T

)

7 T=T+CORR

8 CONTINUE

FINOTV=T

RETURN

END

14
SUBROUTINE PROPS ( PP , DO , TT

)

DIMENSION X(339

DIMENSION B(33) ,G( 32)

EQUIVALENCE (BtX)

COMMON/OATA/G, RtGAMMA

COMMON/l/B

0ATA(ID=1)

DATA (IZ = 1)

1 CONTINUE

IFIIZ.LE.OJGO TO 2

IZ = 0

2 CONTINUE

D = DO

P = PP

T = TT

GM=GAMMA

02=D*D

03=D2»0

0^=03*0

D5=0U»0

06=05»0

07=06*0

08=07»0

09=08*0

010=09*0

Oll=D10*D

012=011*0

013=012*0

TS=SQRT (T)

T2=T*T

T3=T2*T

T<* = T3*T

T5 = Ti»*T

F=EXP CGM*D2)

GO TO ( 10 (,200 , 300 , tfOO ,500 ,600 *700) ,K

ENTRY PRESS

C ENTRY FOR PRESSURE, INPUT IS DENSITY

C AND TEMP. IN MOL/L AND K, OUTPUT IS IN ATM.

K=l

60 TO 1

100 CONTINUE

BC 1)=D2*T

B( 2)=02*TS

B( 3)=02

B( i*)-DZ/l

B( 5)=D2/T2

B( 6)=D3*T

B( 7)=D3

B( 9)=D3/T

B( 9)=D3/T2

B(10)=Oi**T

B(ll)=Oit

B(12>=0*»/T

BC13»=D5

B(1^») = D6/T

B(15)=06/T2

8( 16)=D7/T

B(17)=D8/T

B(18)=08/T2

B(19)=D9/T2

BC20 »=D3*F/T2

BI21»=03*F/T3

B(22)=05*F/T2

B(23)=05*F/Tt»

B(2<»)=07*F/T2

B(25)=07*F/T3

B(26)=D9»F/T2

B(27)=D9*F/Ti»

B(28)=011*F/T2

Bf 29I=D11*F/T3

B(30)=O13*F/T2

B(31) =D13*F/T3

B(32)=D13*F/T<» IFdO.GT.OJ GO TO 102

B(33)=P-R*0*T

RETURN

102 P=0

M=32

DO 101 1 = 1,

M

15
lai P=P*Bf I) GCI)

P=P4-R»0»T

pp=p

RETURN

ENTRY DPOD

C PARTIAL OF PRESSURE WITH RESPECT

C DENSITY - SEE PRESSURE

C ENTRY FOR UNITS

K=2

GO TO 1 200 CONTINUE

Fl=2.00»F»GM»O

F21= 3. 000»F*O2 4-Fl*03

F22 = 5. OOO^F^O'f Fl'DS

F23=7. 000»F»D6 F1»D7

F2<»=9, 000»F»O8 Fl^DS

F25=11.0 0*F»D10*Fl*Oll

F26=l3.0 0»F»O12*Fl»O13

at 1)=2,00»D»T

B< 2)=2.00»D»TS

B( 3)=2.00*O

BC t*)=2.00»O/T

Bl 5)=2.00»D/T2

Bl 6»=3.00»O2»T

B< 71=3.00*02

B< 8)=3.00*D2/T

Bf 9)=3.00»O2/T2

Blia)=«t. 00»D3»T

B(ll)=«». 00*03

8( 12) =^4. 00*D3/T

B<13)=5,00*0*»

B<lt*J=6.00*D5/T

B(15)=e.OO»05/T2 B(16)=7.00»O6/T

Bll7)=6.00*O7/T

B<18)=8.00*D7/T2

B(19)=9.00»O8/T2

8(20)=F21/T2

BI21S=F21/T3

Bf22MF22/T2

BC23)=F22/T^

B(2'») =F23/T2

BI25)=F23/T3

BC26)=F2*»/T2

B(27)=F2'f/Tt*

B(28)=F25/T2

BI29) =F25/T3

B«30» =F26/T2

B(31)=F26/T3

B<32»=F26/Ttf

M~ 3 2

IFaO.GT.0)GO TO 2C2

B(33)=P-R»T

RETURN

202 P=0

00 201 1=1,

M

201 P=P<-B(I) 'GC I)

P = P4-R»T

PP=P

RETURN

ENTRY DPDT

C PARTIAL OF PRESSURE WITH RESPECT

C TO TEMPERATURE - SEE PRESSURE

C ENTRY FOR UNITS

K=3

GO TO 1

300 CONTINUE

X( 1)=02

X( 2)=D2/<2. 00*TS)

X( 3)=0

X( i»J=-02/T2

X( 5)=-2.00*O2/T3

X( 6) =03

X< 7)=0

XI 8I=-D3/T2

X( 9)=-2.00»D3/T3

X(10)=D'»

XI 11)=0

XI12J =-0'*/T2

16
XI13J =0

X(lt»)=-06/T2

X( 15) =-2.00*O6/T3

Xll6>=-07/T2

Xtl7)=-08/T2 X<18)=-2.00»D8/T3

Xll9»=-2. 00*D9/T3

X(20)=-2.0a»O3*F/T3

XC21>=-3.00»D3»F/Tit

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

X(23)=-**. 00»O5*F/T5

Xt2'*)=-2.00»O7»F/T3

XC25>=-3.00»O7»F/Tif

X(26)=-2. 00»O9*F/T3

X(27) 00*D9»F/T5

X(29)=-2.00»D11*F/T3

XI 29)=-3. 00»Dll»F/Tt»

X(30)=-2.00»D13*F/T3

X(31)=-3. OO^OIS^F/Ti*

X(32)=-'».00»D13»F/T5

IF (lO.GT.OJGO TO 302

X(33)=PP-R»0

RETURN

302 P=0

00 301 1=1,32

301 P=P*G(I) 'Xd)

PP=P+R*0

RETURN

ENTRY DSDN

C PARTIAL OF ENTROPY WITH

C RESPECT TO THE G COEFFICIENTS

GO TO 1

kQQ CONTINUE

C S=S0-R»LOGF(O»R»T/PQ) (DSONO) -DSON(O) )*i01.325 +CPOS(T)

Gl=F/(2. 00*GM)

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

G3=:(F»D4-'», 0 0»G2)/ (2.0 0»GM)

G«»=(F»06-6. 0 0*G3)/ (2. OO'GM)

G5=(F*D8-8.0 0»Gi»)/(2.0 0»GM)

G6=(F»D10-10.0 0»G5) /(2 .OO^GM)

X( 1)=-D

X( 2)=-0/ (2. OO'TS)

X( 3)=0.00

X( i»)=+D/T2

Xf 5)=2.00»D/T3

X( 6)=-D2/2,00

X( 7)=0.O0

X( 3)=02/<2. 00»T2)

X( 9J=D2/T3

X(10J=-O3/3.00

X(11)=O.OG

XC 12)=D3/ (3. 00»T2)

Xf 13J =0,00

X(lif)=O5/{5.00»T2) X(15)= 2. 00*05/(5. 00»T3)

X(16)=D6/(6.00»T2)

Xf 17)=D7/I7.00»T2I

X(18)=2.00*D7/(7.00*T3)

X(19) =D8/f<t, 00*T3>

X(20)=2. 00»G1/T3

Xf 21) = 3,00»G1/Tt» X{22)=2.00*G2/T3

X(23)='t. 00»G2/T5

X(2<»)=2. 00»G3/T3

Xf 251=3. 00»G3/T*f

X(26)=2. 00»G'*/T3

X(27) =i*. 00»G«*/T5

X(28)=2. 00»G5/T3

Xf 299 =3. OO'GS/Ti*

X(30)=2.00*G6/T3

X<31)=3. 00»G6/Ti»

Xf 32)=*f. 00»G6/T5

IF<ID.GT.0IGO TO ^^02

RETURN

U02 P=Q

00 tfOl 1 = 1,32

kQl P=P*G(I) Xf I)

PP=P

17
RETURN

ENTRY DUDN

C TERMS NEEDED FOR ENTHALPY CALCULATION

K=5

GO TO 1 500 CONTINUE

C H = HO f (T»DSON<0) -OSON( 0 )

) • 101. 325 ( DUDN ( D-OUON ( 0 ) ) 1 0 1 . 325 + CPOH ( T)

C (P/0-R»T» 101.325

G1=F/ (2.00»GM>

G2=CF»O2-2.0 0»Gl) / (2. 0 0»GM)

G3=(F»O<*-<*.00»G2)/(2. 0 0*GM)

G^fs <F» 06-6, 00»G3)/ (2.0 0*GM)

G5=(F*D8-8.00»G**)/i2.00»GM)

G6= ( F» Dl 0- 10.00* G5)/ (2.0 0 'GM)

X( 1)=D»T

X( 2)=0*TS

XI 33=D

X( ^)=0/T

X< 5)=D/T2

X( 6)=D2»T/2.00

XI 7J=D2/2.00

X( 8)=D2/(2.00*T)

X< 9) =D2/(2.00*T2»

XllO ) = O3»T/3.00

Xlll9=D3/3.00

X(12)=03/(3.00»T)

X(13) =D'»/<».00

X(l*»)=O5/<5.00»T)

Xll5l=O5/C5.00»T2J

X(16)=O6/(6.00*T)

X(17)=D7/<7.00»T)

X(18)xD7/<7. 00»T2»

Xll9}=O8/C8.00»T2l

X(20)=G1/T2

X(21)=G1/T3

X(22)=G2/T2

X( 233=G2/Ti*

X(2'»)sG3/T2

X(25)=G3/T3

X(26)*G<»/T2

XI27)=G^/T%

X(28)=G5/T2

X(29)=G5/T3

X(30)=G6/T2

X(31)=G6/T3

X(32)=G6/T^f

IF(IO.GT.0»GO TO 502

RETURN

502 P=0

00 501 I=l»32

501 P=P*G(I)»X(I)

pp_p

RETURN

ENTRY TDSDT

C TEMP. TIMES THE PARTIAL OF

C ENTROPY WITH RESPECT TO TEMP.

K=6

GO TO 1

600 CONTINUE

C CV=CV0+(TDSDNI/)-TDSON(D) ) '101.325

G1=F/(2.00»GMJ G2=iF»02-2.0 0»Gl)/ 12. 0 0*GM)

G3= (F»O*t-<».0 0*G2)/(2. 0 0*GM)

G«»=(F*O6-6.00»G3)/(2.0 0»GM)

G5=(F»08-8.00»G*t)/ (2.0C»GM)

G6=«F»D10-10.0 0»G5) / (2.0 0»GM)

X(i)=0

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

X(3)=0

XI tfJ=2.00»D/T2

X( 5)=6.00»D/T3

X(6)=0

X(7)=0

XI 8)=02/T2

X( 9)=3.00»O2/T3

X(10)=0

XI11)=0

Xll2»=f 2.00*03}/ 13. Q0*T2)

X(13)=0

18
X(l<» )= (2.00»D5)/(5.00»T2)

X(l5)=(6.0 0»O51/(5,00»T3)

X(16 )=D6/C3. 00»T2)

X ( 17 »= 12.00*078/ 17. 00»T23

X(18>=(6.00»O7}/(7.00»T3)

X(19) = (3.00*O8l/<<*.00»T3)

XC20)=6.000»G1/T3

X (21)=12.00»G1/T«»

X(22)=6.000»G2/T3 X{23)=20.00»G2/T5

XI2<»J=6.000*G3/T3

X(25)=12.00*G3/T*» X(26)=6.000»Gi»/T3

X(27)=20.00»G4/T5

X«28)=6.000»G5/T3

X(29)=12.00*G5/T4

X(30) =6.000»G6/T3

X(31)=12.00»G6/T*»

Xf 32)=20.00»G6/T5

IFCIO.GT.OIGO TO 602

RETURN

602 P=0

00 601 1=1,32

601 P=P*G(I)»XII)

PP=P

RETURN

ENTRY 0P202

C SECOND PARTIAL OF PRESSURE WITH

C RESPECT TO DENSITY SQUARED

K=7

GO TO 1

700 CONTINUE

F1=2.»F»GM»D

Fl 2=2. *F1*GM*D*2.»F*GM

F212=3.»F1»02+3.*2.*0»F+F12»03+F1»3.*02

F222 = 5.»Fl»Di» +5. '.03»F + 5. Di»*Fl + F12*D5

F232 = 7.»F1»06*7.»6. D5»F + 7.*06»F14-F12»07

F2«»2 = 9.*F1»08*9.»8.*07*F*9.»08*F1*F12*D9

F25 2=ll.»Fl*O10*10.*ll.*O9»F*ll.*D10*Fl+F12*Dll

F262=13.»F1»D12*13.»12.»011*F+13. *012*F1+Fi 2*013

BI1)=2.»T $BC23=2.»TS $ BI3J=2.

B(<»)=2./T $ B(5) = 2./T2 $ B(6)=6.»D»T

BC7)=6.*0 $ B<8I=6.»D/T $ B(9)=6.»0/T2

B(10)-12.»O2»T $ B(11)=12.»D2 $ B ( 12) =12. *^02/T

B<13)=20.»03 $ 8(1*»)=30.»D<»/T $ B ( 15) =30 . *D8t/T2

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

B<19)=72.»07/T2 $ B ( 20 » =F212/T2 $ B< 21) =F212/T3

8(22)=F222/T2

B<23) =F222/T^» | B(2^f ) =F232/T2 $ B ( 25 ) =F232/T3

B(26)=F2*»2/T2 § B (27 ) =F2<»2/T«t $ B ( 28 ) =F252/T2

B(29)=F252/T3 $ B ( 30 ) =F262/T2 $ B ( 31 ) =F262/T3

B(32)=F262/T1»

M= 32

IF(IO.GT.0)GO TO 702

B(33)=PP

RETURN

702 P=0

DO 701 1=1,

M

701 P=P+B(I) G(I)

PP=P

RETURN

END

19
FUNCTION FIND OfP,T)

C ITTERATES EQUATION OF STATE

C FOR DENSITY, GIVEN PRESSURE

C ANO TEMP. IN ATM. AND KELVIN. IF

C ITTERATION FAILS TRY USING

C FUNCTION CALLED FIND M

COMMON/OAT A/G,R, GAMMA, VP, OTP

DIMENSION G(32) , VP(9)

TT = T PC = VPN<VPI8H

IF(TT.GT. VP(81 )G0 TO IOC

IF( P.GT. VPN(TT) )G0 TO 101

X=P» (.8-l.J /PC+1.0

00=P/<R»T»X)

GO TO 102

100 X=(1.1/(9.»PC) )»P+.7/9.

00=:P/(R*T»X)

IF (P/PC.GT. 20. .AND.T/VPIS) .LT.2.5) 00= DTP

GO TO 102

101 00=OTP

102 CONTINUE

DO 10 1=1,50

CALL PRESS«PP,OD,TT»

P2 = PP

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

1 CALL DPODIPP,DD,TT)

DP=PP

C0RR=(P2-P) /DP

0=00

IF(ABS CC0RR)-1.E-7»D» 20,20,10

10 OD=OD-CORR

FIND D=0

RETURN ^ 20 FIND 0=00

RETURN

END

FUNCTION FIND M(P,T,DD)

C ALTERNATIVE FOR FIND D, INPUT IS

C PRESSURE IN ATM., T IN KELVIN AND

C DENSITY IN MOL/L. INPUT DENSITY

C IS A STARTING VALUE FOR ITTERATION

C OF EQUATION OF STATE FOR SOLUTION FOR P AND T

TT = T

00 10 1=1,50

CALL PRESS(PP,OD,TT)

P2=PP

IF(AeS (P-P2)-l.E-7»P)20t20,l

1 CALL OPDDfPP,DD,TT)

OP=PP

C0RR=CP2-PI /DP

0 = 00

IFIABS CC0RR»-1.E-7»D) 20,20,10

10 D0=D0-C0RR

FIND M=0

RETURN

20 FIND M=DD - RETURN

END

20
FUNCTION ENTHAHP.OtTJ

C CALCULATES ENTHALPY FOR INPUT OF

C PRESSURE, DENSITY AND TEMP. IN

C ATM.» MOL/L AND K. OUTPUT IS IN

C JOULES/MOL

R=«08205616

D0=0

TT = T

CALL DSDNCSO,OD»TT)

CALL DUDNCUO«DD«TTI

00=0

CALL DSDN(SOf DO,TT)

CALL DUDN(U0«DO«TT)

ENTHAL=T*CSD-SO) 101. 325+ (UD-UO) 101. 3 25+HI (T) CP/D-R»T3»101. 325

RETURN

ENO

FUNCTION ENTROP(D,T)

C CALCULATES ENTROPY

C FOR AN INPUT OF DENSITY AND

C TEMP. IN MOL/L AND K. OUTPUT IS IN

C JOULES/MOL-K

R=,08205616

00=0

TT=T

CALL DSDN(SD,DD.TT)

DD = 0

CALL DSONCS0tOO,TT)

ENTROP= (SO-SO) 101.325-R*ALOG(O*R»T) *101.325+SI<T)

RETURN

ENO

FUNCTION SOUNOCDiTi

C CALCULATES THE SPEED OF SOUNO

C FOR AN INPUT OF DENSITY AND TEMP.

C IN MOL/L AND KELVIN. OUTPUT IS IN

C METERS/SECOND.

COMMON/CRIT/W

CALL DPDD(OP,D,T)

SOUNO=HCP<DtTI/CV(0,T)) DP»101325./W)»».5

RETURN

END

FUNCTION CP(0,T)

C CALCULATES SPECIFIC HEAT CAPACITY

C AT CONSTANT PRESSURE FOR INPUT OF

C DENSITY AND TEMPERATURE IN MOL/L AND K

C CP IS IN JOULES/MOL-K

CVEE=CV(0,T)

CALL DPDT<DPT»D,T»

CALL DPDDIDPCDjT)

CP=CVEE4- (T/CD»»2)»<DPT»*2) /DPD)»101. 32 5

RETURN

ENO

21
FUNCTION CV<OiT)

C CALCULATES SPECIFIC HEAT CAPACITY

C AT CONSTANT VOLUME FOR AN INPUT

C OF DENSITY AND TEMPERATURE IN MOL/L AND K

0ATACR=8.31'»3<»»

00 = 0

TTsT

CALL TDSOTIC0,DO,TT»

00=0

CALL TOSOTIC0»OD,TTJ

CV = CPI(TT)4- (CO-CD) 10 1.325

CV=CV-R

RETURN

ENO

FUNCTION CPUT)

C CALCULATES IDEAL GAS THERMO PROPERTIES

C FOR ALL FLUIDS EXCEPT H2. INPUT IS IN

C KELVIN, OUTPUT IS IN JOULES, MOL/L AND K

C0MM0N/CPID/G(11)

COMMON/ISP/N

K=l

1 U=G(9) /T .

. EU=EXP lUJ - TS=l./T»»i*

GO TO (2>3,^) vK

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

DO 10 1=1,7

TS=TS»T

10 CPI = CPH-GCI)*TS

CPI = CPI»8.31'»3<»

RETURN

ENTRY SI

K=2

GO TO 1

3 CPI=G(8)»(U/<EU-1.)-AL0G«1.-1./EU)) l-Ga»»TS»T/3,-GI2) TS»T*T/2,-GC3l/T4-GC<») ALOGC T) -J-G ( 5) *T+G «6 ) *T»T/

2.+G(7)»T*»3/3.

CPI=CPI*8«31U3^*G (11)

RETURN

ENTRY HI

K=3

GO TO 1 < CPI=G(8) »U»T/IEU-1.) -G<1) /(2.»T»T) -G{2) /T+G«3) »ALOGIT) G{t»)»T

1 + G (5»»T»T/2, Gie) T»»3/3.*GC7I T

CPI=CPI»8.31<»3i» + G(10)

RETURN

ENO

22
FUNCTION VISC{O0,T)

C RETURNS VISCOSITY IN ( G/CM-S)

,

C T IN K, D IN MOL/L

COMMON/CRIT/EM

O=OO*EM/1000.

VISC = DILV(T) FOC\;(D,T) +EXCESV{D,T)

RETURN

END

FUNCTION DILV(T)

C GIVES DILUTE GAS VISCOSITY AND THERMAL

C CONDUCTIVITY FOR AN INPUT TEMP. IN

C KELVIN. OUTPUT UNITS ARE SAME AS

C THOSE IN Vise AND THERM

COMMON/ISP/N,NW,NWW

C0MM0N/0ATA1/GV,GT,FV,FT,EV,ET

DIMENSION GV(9I.GT(9),FV('t),FTJitJ,EVf8),ETC8)

SUM=0

TF=T»»Cl./3.

J

TFF=T*»(-U./3.

I

DO 10 1=1,9

TFF=TFF»TF

10 SUM=SUM+GV(I)*TFF

IF(NWW.EQ.7) GO TO 9

DILV=SUM»1000.

GO TO 11

9 DILV=SUM

11 RETURN

ENTRY DILT

TF=T»»(l./3.)

TFF = T*»(-tt./3. ) SUM=0

DO 20 1=1,9

TFF=TFF»TF

20 SUM=SUM+GT( I) TFF

DILV=SUM

RETURN

END

FUNCTION FDCV(0,T)

C FIRST DENSITY CORRECTION

C FOR VISCOSITY AND THERMAL CONDUCTIVITY

C0MM0N/DATA1/GV,GT,FV,FT,EV,ET

DIMENSION GV(9) iGT(9) , FV{4) ,FT<tt) ,EV(8J ,ET(8) F0CV=(FV(1) FV(2)» (FV(3J-AL0G(T/FV(if) ) )^*2)*D RETURN

ENTRY FOOT

F0CV=(FT(1) +FT(2) {FT (3) -ALOG ( T/FT (^) ) )**2)»D RETURN

END

23
FUNCTION EXCESV(0,T)

CALCULATES EXCESS VISCOSITY

COMMON/0 AT A1/GV,GT,FV, FT tEV,ET

COMMON/ISP/N,NW

DIMENSION GVf9) .GTCg) ,FVf«») ,FT<^» ,EV<8} ,ET(6)

R2=D»»( .5) »<ID-EVI8) ) /EVI8)

)

R=0»» ( •!)

X=EV(1) EVCZ) RZ + EVtS) *R+EV(i*) RZ/ ( T*T ) +EV( 5) R/T^^d.S) 4-EV(6) /T

l+EViZ) *R2/J

X1=E V(l) +EV<6) /T

EXCESV=EXP IX5-EXP (XI)

RETURN

ENTRY EXCEST

CALCULATES EXCESS THERMAL CONDUCTIVITY

IF(NW,EQ,0 I GO TO 3

R=0»»l.l)

X=ET (1) ETC 2) R+ET(3)*R/T*»(1.5) ETC*) /T

X1 = ET(1)+ETU»/T

EXC£SV=CEXP (X)-EXP (Xl))/10,

RETURN

J R2=0**(65)»((D-ET(8) ) /ET(8)

)

R^O^ * ( • 1

)

X=ET<l)*ET«2l»R2 + ET(3)*R+ETi'») R2/(T»T) +ET(5) 'R/T*' il . 5) ET (6) /T

1+ET (7)»R2/T

X1=£T(1) ETCe) /T

EXC£SV=EXP «X)-EXP (XI)

RETURN

END

FUNCTION THERMCOO,T)

RETURNS TC IN MW/M-K, T IN K, 0 IN MOL/L

COMMON/HAN/CR,TCI

COMMON/ISP/N.NM

COMMON/CRIT/EM

O=DD»EM/lO00.

IF(NW.EQ.O ) GO TO 3

CR=CRITC(OiTl

THER=DILT(T) 4.FDCT(D,T) 100 .^-EXCESTCO.T) +CR

TCI=THER-CR

TH£RM=THER

RETURN

CR=CRITC(0,T)

THERM=OILT(T) FDCTtO.T) EXCEST «D , T) CR

TCI=THERM-CR

RETURN

END

FUNCTION SENGI0,T3

SCALED EQUATION OF STATE FOR CRITICAL REGION

COMMON/CRIT/ EM, EOK, RM, TC, DC, X , PC, SIG

COMMON/SEN/BETA, XO, DELTA, El, E2, AGAM

COMMON/CHECK/DELDf OELT.DSTAR.TSTAR

DSTAR= D/OC $ TSTAR=T/TC

BET0=1«/BETA

XX=0ELT/DEL0*»BET0

AG=AGAM-1,0

BET2= 2.0»BETA

AGB=AG/BET2

OEL1=OELTA-1.0

AGBB=«AG-6ET2)/BET2

XX0=fXX* XOI/XO

XXB=XX0*»eET2

BRAK=1«0 E2*XXB

8RAK1=BRAK»»AGB

H=E1»XX0*BRAK1

HPRIM= (E1/X0)*BRAK1 ( AG/XO) E1*E2*< XXB) *( BRAK**AGBB)

RC0M=<DEL0»*0EL1)»<DELTA»H - (XX/BETA ) HPRIM ) RCOMP=1.0/<RCOM»DSTAR»»2)

RCM=RCOMP/«PC»1.01325E+06)

RCM IN DYNES

RCM=RCM»*0,5

SENG=RCM

RETURN

ENO 24
FUNCTION CRITCID.TI

C CALCULATES CRITIC&L ENHANCEMENT

C FOR THERMAL CONDUCTIVITY

COMMON/GRIT/ EM, EOK, RM, TC, DC, X , PC, 5IG

COMMON/CHECK/OELO,OELT .0STAR,TSTAR

COMMON/HJM/EPSI,CPCV,RRR, AKT

COMMON/ISP/N,NW,NWW

C D IN G/CM3 • T IN K

C OUTPUT UNITS ARE MW/M.K

A\/ = 6,0225E + 23 S BK=1. 3805^*E-i6

DEL0=A8S <0-DC)/DC $ DELT=ABS (T-TC)/TC

0 CALCULATE DISTANCE PARAMETER

R=(RM»»2.5) <D»»0.5> *< AV/EM)»»0.5

R = R* <EOK»»0.5) X/(T»»0.5)

C GENERAL EQUATION

DX=0»1COO.O/EM

CALL OPDT<DPT,DX,T)

C DPDT IN ATS PER DEG.

DPT=DPT»1.01325E+6

C OPOT NOW IN DYNES PRR DEG

CALL DPDO(OPO,OX,T)

C DPOD UN ATS, MOL/L

0PO=OPD»<1.01325E*6) •( 10 00.0/EM)

COPDO NOW IN DYNES, GM/CM3

IFC DPD.LT.O.Q) DPD=1,0

9** VIS=VISC(OX,T)*(1.0E-06)

C VISCOSITY IN GM/CM.S

IF(OELD.EQ.0.25. OR. OELD.LT. 0.25) 8,10

8 IF{OELT.EQ.0.025. OR. DELT. LT. 0. 025 ) 9, 10

9 COMPRES=SENG<0,T)

GO TO 12

10 COMPRES=1.0/IO»DPD) 0.5

12 EX=BK»T»*2»I0PT»»2) 'COMPRES

EXB=R»{ (BK*T)»»0.5>»(O»*0.5)»( (AV/EM) O.5)

CRIT = EX/ (EXB»6.0*3.1'tl59*VIS>

C PUT IN DAMPING FACTOR

300=( <0-DC) /DCi**it

BTT= ( (T-TC) /TC3 2

FACT=EXP (-18.66»BTT - ^.25*3D0)

OELC=CRIT»FACT/100.0

CRITC=DELC

AKT=COMPRES*COMPRES

EPSI=R*R»BK»T*«AV»0/EM)*AKT

EPSI=EPSI»»0.5

C CALC CP-CV

CPCV=T» (0PT»»2)*AKT/D

RETURN

END

25
I

Appciul I X W.

PROGRAM VALUES<OUTPUT, INPUT)

CALL DATA 02

PRINT 20

20 F0RMAT(82H1 P T RHO

1 C-V C-P VEL/82H ATM

2 J/MOL J/MOL-K J/MOL-K J/MOL-K

00 19 1=1,5

READ 5,P,T

5 F0RMAT«F9.3,F8.3)

DEN=FINDO(P,T)

H=ENTHAL (Pf 0EN,T)

S=ENTR0PC0EN,T9

CVE=CV(0EN,T)

CPE=CP<OEN,TI

VEL=SOUNDfOEN,T)

PRINT 13,P,T,0EN,H,S,CVE,CPE,\/EL

13 FORMAT(F10««»,F12,5,F10.6,F12.1,3F10.2,F8.0)

19 CONTINUE

END

KELVIN

M/SEC/IH ) MOL/L

P T RHO H S C-V :-p VEL

ATM KELVIN MOL/L J/MOL J/MOL-K J/MOL-K J/MOL-K M/SEC

50,0000 60.00000 40.237637 -5780.0 71,95 35.00 53.16 1127.

50.0000 100.00000 35.252617 -3480 .5 97.93 29.39 52.56 902.

.40.0000 200.00000 13.500484 3013.2 141.86 24o63 71.41 295.

00.0000 200.00000 22.282750 2062.6 132.73 24.40 54.66 497.

00.0000 300.00000 24.548716 6948.2 140. 30 24o72 40. 42 736.

26
0.04

-0.04

0.04

>-t 0

a -0.04

2 0.08

<

S 0

120K

J \ I I I I I I I I \

I Mill

0.08

0.4

0.2

0

-0.2

-0.4

J I I I II

I

J L 140K

0 •--v.v

o Gas

* Liquid

I I I II I II

150K

I i I I 1 1

1

0.2

0

-0.2

0.1

0

-0.1

J I I I LI 170KJ i

I I I 1 1i

J L I I I I I I I 200KJ \

I I I 1 1

1

10 100

PRESSURE, atm

-0.2

1000 10

300K

-i__J L_I.I,.Li..LL ...I J..i_l.l.lJ. 100 1000

PRESSURE, atm

Figure 1. Deviations of the density data from eqn. (1) on several

isotherms

.

27
CO 9

CP

Om O • O • O t OOOOOOOOO

4doooooooooooo i o • o • o

•o

•occ

•oo oo

fiC o

00

CO

CP % '33N3U3dJIQ AllSNBQ

03 1

CO>

0) ,c S-l u 0 o >

fr 1

CO o

0) •H

4-1 T3

•H •H

CO 3 d XT HH d 1 o

•H

O

•u • #t

iH

4-1 ojCO c

(U

(1) 4-1 c« U-l •H o >

COd OJ o 4J •H cC

4-1 rH

nj :3 •H o > H

OJ Q CJi

CM i

<U

00

•H

28
0

^

1 6.1 mol/l -

150

oo

o

I

o

J L J L 200 250 300

5 r

4

3 -

2 -

1 - oo o

<^ o

° o

150

J I Qq 9

200 250

13.1 mol/l

300

O Oo

0 o Q o

° 0°°°* oo°o

° 28.6 mol/l

120 140 160 180 200

TEMPERATURE, K

Figure 3. Deviations of the heat capacity data on three isochores

from values calculated by eqn. (1)

.

29
1

o

o

o

' o

o

o

o

o

o

o

o

o

o

^ CD

r*^ O O O

to m cnvn

CM

I I I I o O <J

O D

43.

to

E

«/3

C/3

€M

CM CM

% '(3|B3 - dx9) jjia Ail3013A QNnOS

<jn \ 4-J <^

1—

1

I-' CJ rH

CJ

UJ

CO p

O

V-i X'-

CO TO •H ^3

rH

t\3 r*\-i I,

•~l fj vU Q K

CO fit

tL_J (0 O

4-J Cy O

O 4J rH CO, > ^\

0

cu u

O i

CO a.

<u: o

•H <U' J-i Xi

cti u

> «d

0) •H>

3

•H

30
8. References

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

Thermodynamic Properties of Light Hydrocarbons and their Mixtures.

I. Methane, Ethane, Propane, and n-Butane, J. Chem. Phys. 8^, 334-45

(1940).

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).

Weber, L. A., PVT, Thermodynamic and Related Properties of Oxygen

2 from the Triple Point to 300 K at Pressures to 33 IW/m , J. Res. Nat.

Bur. Stand. (U.S.) 74A , No. 1, 93-120 (1970).

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

Triple Point to 300 K at Pressures to 1000 Bar, NBSIR 77-865, NASA

RP 1011 (Dec 1977).

Goodwin, R. D. and Weber, L. A., Specific Heats of Fluid Oxygen

from the Triple Point to 300 K at Pressures to 350 Atmospheres, J. Res.

Nat. Bur. Stand. (U.S.) 73A , No. 1, 15 (1969).

Stewart, R. B. and Jacobsen, R. T., A Thermodynamic Property

Formulation for Oxygen from 55 K to 350 K, Proc. of the Seventh

Symposium on Thermophysical Properties sponsored by ASME (May 1977) , to

be published.

The International Practical Temperature Scale of 1968, Metrologia _5, No. 2, 35 (1969).

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 (1962).

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

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

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

Parahydrogen, Nat. Bur. Stand. (U.S.), Internal Rept. 74-357 (1974);

also, Roder, H. M. and McCarty, R. D., A Modified Benedict-Webb-Rubin

Equation of State for Parahydrogen-II , Nat. Bur. Stand. (U.S.), Internal

Rept. 75-814 (1975).

31
[11] Weber, L. A., Density and Compressibility of Oxygen in the Critical

• 1 Region, Phys. Rev. A6, No. 6, 2379 (1970).

[12] Goodwin, R. D., Nonanalytic Vapor Pressure Equation with Data for

Nitrogen and Oxygen, J. Res. Nat. Bur. Stand. (U.S.) 73A , No. 5, 487

^ (1969).

3 k

[13] Mills, R. L. and Grilly, E. R. , Melting Curves of He , He , H^ , D„,

2

2 z

. . Ne, N^, and O2 up to 3500 kg/cm , Phys. Rev. 99, No. 2, 480 (1955).

[14] CODATA Bulletin 10, CODATA Recommended Key Values for Thermodynamics,

1973, Report of the CODATA Task Group on Key Values (Nov 1973).

[15] McCarty, R. D. , private communcation.

[16] van Itterbeek, A. and van Dael, W. , Velocity of Sound in Liquid

Oxygen and Liquid Nitrogen as a Function of Temperature and Pressure,

Physica 28, No. 9, 861 (1962).

[17] Straty, G. C. and Younglove, B. A., Velocity of Sound in Saturated

and Compressed Fluid Oxygen, J. Chem. Therm. 5^, 305 (1973).

32
NBS-n4A 'WEv. 7-73)

U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET

I. 1 t i^l.lv .\liV;i\ VM\ IVi.l IM\1

NBSIR 78-882

2* C'ov I Al l ("'.•.(oil No.

3, l^ci i|>u'i)t .s Av ( f'. sittii No.

4. 1 II 1 1 A N 1 1 n 1 1 1 1 . r A MODIFIED BENEDICT-^sT-BB-RUBIN EQUATION OF STATE FOR

GASEOUS AND LIQUID OXYGEN

5. 1 'uM It .11 u»n n.M

f

April 1978

736

7. .MJrHOR(S)

L. A. Weber

8. Performing Organ. Report No.

9. pi-;ki-()KMIN(i ()K(i an i/.ation .NAM1-; ANO Anni<i-:.s.s NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

10. Projc-ct/Ta.sk/Work Unit No. 7360423

11. Contract/Grant No.

12, Spon.soring Organization Name- and ( tjmplctc Addrcs.s (Street, City, State, ZIP)

NASA-Johnson Space Center

Houston, Texas 77058

13. Type of Report & Period

Covered

T4. Sponsoring Agency Code

15. .slji^pi.i:mi:ntary noils

16. .ABSTRACT (A 200-word or less factual summary of most si^sficant information. If document includes a significant

bibliography or literature survey, mention it here.)

New PVT data on oxygen at pressures to 800 atmospheres have been recently reported.

These data, along with other NBS data on oxygen, have been fitted with a 32 term

modified Benedict-Webb-Rubin equation of state. This equation is valid for gaseous

and liquid oxygen at temperatures from the triple point (54.36 K) to 300 K and for

pressures up to 800 atm. ' Thermophysical properties calculated via this equation are compared with the most reliable values available. The vapor pressure equation has

also been updated to be consistent with more recent experimental work near the triple

point, A computer program is available which calculates the thermophysical properties of oxygen.

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) Density; enthalpy; entropy; equation of state; oxygen;

property computer program; PVT; saturation properties; specific heat.

18. AVAILABILITY

^

Unlimited

! For Official Distribution. Do Not Release to NTIS

I i Order From Sup. of Doc, U.S. Government Printing Office Washington, D.C. 20402, SD Cat. No. C15

rX Ord er From National Technical Information Service (NTIS)

Springfield, Virginia 22151

19. SECURITY CLASS

(THIS REPORT)

UNCL ASSIFIED

20. SECURITY CLASS

(THIS PAGE)

UNCLASSIFIED

21. NO. OF PAGES

38

22. Price

$4.50

USCOMM.DC 29042-P74