NBSIR 81-2253

’ • I A THERMODYNAMIC SURFACE FOR

WATER: THE FORMULATION AND

COMPUTER PROGRAMS

NATIONAL BUREAU of STANDARDS

LIBRARY _ ~ WJG 8 1983

f\o\OJX

GlC vOO

v^3 X

Lester Haar

John S. Gallagher

U S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Measurement Laboratory

Washington, DC 20234

George S. Kell

National Research Council

Ottawa K1A 0R9, Canada

June 1 983

U.S. DEPARTMENT OF COMMERCE, Malcolm Baldrige. Secretary

NATIONAL BUREAU OF STANDARDS. Ernest Ambler. Director

A Thermodynamic Surface for Water:

The Formulation and Computer Programs

by

Lester Haar, John s. Gallagher and George S. Kell

1. Introduction

This report has been prepared to provide computational details and

computer programs for the evaluation of the thermodynamic properties of

water and steam. It is based on an extensive correlation of thermodynamic

properties. 1 Included herein are the basic equations for the description

of the thermodynamic properties of water and detailed listings of computer

programs based on these equations. With this information, thermodynamic

properties can be calculated from the freezing point to '2500 K in

temperature and from the dilute gas to well in excess of 1 GPa (10,000

bars) in pressure for liquid and gaseous states for undissociated water.

In reducing the data the molecular weight of ordinary water has

been taken as 18.0152 g mol" 1

, 2 the universal gas constant as 8.31441

J mol -1

K" 1

, and hence the specific gas constant as 0.461522 J g" 1 K" 1

. 3

Absolute temperatures on the International Practical Temperature Scale

of 1968 have been used as the realization of the thermodynamic scale.

This report contains five sections and 2 appendices: A brief

discussion of the structure of the thermodynamic surface for water is

given in section 2, and in section 3 are the equations, parameters and

constants that define the surface. Section 4 contains the thermodynamic

relations used in the computer program. It also contains a discussion
2

of the general organization and structure of the computer program.

Section 5 contains statements regarding the accuracy and range of

validity of our surface. In Appendix A are the computer listings and a

quide to their use, and, lastly, in appendix B, is a printout of calculations

obtained with the program, so that the user can check his results.

2. Structure of the thermodynamic surface

The thermodynamic surface is composed of three parts: (1) The

first, referred to as the base function, is obtained from a theoretical

equation of state. 4- 12 It properly describes the low-density vapor, the

high temperature behavior at all densities, and the dense fluid at all

temperatures. Except for a large region around the critical point, the

base function yields results that are at least in semi -quantitative

accord with data. (2) The second, referred to as the residual function,

yields corrections to the base function. These corrections are small in

regions where the base function is in good accord with data, and in

regions beyond the range of the data. The contributions to the pressure

from the base function and the residual function are readily integrable

in closed form to yield their respective contributions to the Helmholtz

function. (3) The third includes the thermodynamic properties for the

ideal gas state. These have been reported recently by Woolley 13 as

part of a detailed analysis of the rotation-vibration structure of the

water molecule.

A major part of the work was given to the derivation of the 40

terms that form the residual function. The first 36 terms were used in

a global, least squares fit to data. Each of these terms yields important

contributions over wide ranges of the independent variables. Following

this, small improvements were made by adding 3 terms that contribute
3

only in the immediate neighborhood of the critical point, and a single

term that contributes only in the region of high pressures and low

temperatures. Except in these very limited regions, the residual

function is given by the first 36 terms. A discussion of the thermodynamic

surface obtained with the residual function so restricted (terms 1-36)

is given in reference [12].

The three parts are combined to obtain the expression for the

Helmholtz function for fluid water.

where the independent variables are density (p) and temperature (T).

Eq. (1) is what we mean when we refer to the thermodynamic surface for

water.

3. The Helmholtz function

Presented in this section are the equations that define the Helmholtz

function Eq (1). The units used for the independent variables are

p(g/cm 3

) and T(K), and R (J g _1

K" 1

) for the specific gas constant, so

that A(p,T) is given in joules/g.

The base function:

where a, £, and y are geometric constants 11, 133/3 and 7/2; y is a

dimenionless density, given by y

= bp/4, and P Q

= .101325 MPa = 1 atmosphere.

0 )

• 5Z|±1 + tn £j*I ( 2 ) 0
4

The parameter b is a temperature dependent hard sphere volume, and B is

related to the second virial coefficient.

where T = 647.073 K and the coefficients b. and B. are listed in table 1. o 11

b

i (cm 3 g' 1

)

Table 1

i B.j (cm 3 g" 1

)

.7478629 0 1.1278334

-.3540782 1 - .5944001

0 2 -5.010996

.007159876 3 0

0 4 .63684256

-.00352836 5 0

The residual function:

+ exp * (5)
5

where a = 1 cm 3 /g, and the g.

data. The quantities

6^ and x.

respectively, given by

P - p.

are coefficients determined by fits to are reduced densities and temperatures 5

and Pi , T. are specified density and temperature values. The k(i) and

n(i) are integers. Values for the constants and parameters for the

residual function are listed in table 2.

Table 2

1 m) 1(1) S(i) i 3c ( 1 > 1(1) fi ( i ) (J 9’ 1

) (J g’ 1

)

1 i l -.530629685290+03 19 5 4 -.138025771779+07

2 i 2 .227449014244+04 20 5 6 -.251099143690+06

3 i 4 .787793330207+03 21 6 1 .465616261156+07

4 i 6 -.698305273750+02 22 6 2 -.727527732754+07

G.w 2 1 .178638328754+05 23 6 4 .417742461483+06

6 2 2 -.395147315633+05 24 6 6 .140163582446+07

7 2 4 .338038842808+05 25 7 1 -.315552313921+07

8 2 6 - . 138550502027+05 26 7 2 .479296663846+07

9 3 1 -.256374366133+06 27 7 4 .409126647812+06

10 Tw 2 .482125759814+06 28 7 6 -.136263693884+07

11 3 4 -.341830169697+06 39 9 1 .696252208627+06

12 •7w 6 .122231564174+06 30 9 2 -.108349000964+07

13 4 1 .117974336556+07 31 9 4 -.227228274017+06

14 4 2 -.217348101104+07 32 9 6 .383654860007+06

15 4 4 .108299521686+07 33 3 0 .688332579443+04

16 4 6 -.254419980640+06 34 3 3

~ .217572455226+05

17 5 1 -.313777749478+07 35 1 3 -.266279448298+04

18 5 2 .529119107577+07 36 5 3 -.707304180821+05

i MO MO P-j(g/cm- 3

) T-j ( K) a

i 6

i 9-j (J g- 1

) 37

38

2 0

9

.319 640. 34 2000Q -.225

39

C C

r\

. 319 640. 40 20000 -1.68

40

L

A U . 319

*1 r r 641.6 30 40000 .055 *T U 1 . 55 270. 1050 25 -93.0
6

The ideal gas function:

where T

R

= T/100 K.

The coefficients C. are given in table 3.

Table 3

1 .197372710180+02

2 .209662661977+02

3 -.483429455355+06

4 .605743189245+05

5 .225602388500+02

6 -.987532442000+01

7 -.431355365532-01

8 .458155781927-04

9 -.477549017624-07

10 .412384608402-10

11 -.279290527404-13

12 .144816952031-16

13 -.564736589529-20

14 .162004460052-23

15 -.330362277656-27

16 .451916066716-31

17 -.370734122641-35

18 .137546067535-39

The Eqs 1 thru 6 and the values for the coefficients and parameters

included in this section contain the complete thermodynamic description

for water. Since A(p,T) is everywhere analytic, it is straightforward

to evaluate appropriate derivatives and to obtain closed-form expressions
for all thermodynamic properties over the entire fluid range. In the

next section is a discussion of how this is done.

4. Method of calculation and the computer program

The thermodynamic quantities are calculated from A(p,T) using the

following relations:

p s pRTZ = p 2 3A/3p (Pressure) (7)

3P

3p

= - p + p 2

P 3P 2 (8)

3P

3T

= 2

3 2A

p 3p 3l* (9)

S = - 3A

3T (entropy) (10)

U = A + TS (internal energy) (11)

H = U + P/p (enthalpy) (12)

C

v

= - T 3 2 A

1 3T2 (isochcric heat capacity) (13)

G = A + P/p (Gibbs functions) (18)

c

p

= r

. T (3P/3TJ2

L

v p2 3P/3p (isobaric heat capacity) (19)

O) - 1(vc

v fY

/2

(speed of sound) (20)

B n

l h

2 P\

2RT V 3p 2 ] p=0 (2nd virial coef.) (21)

dB n

dl (22)

6

t

- |

I 3h\ = 1 _ T 3P/3T

^Pyj P

p'2

3 P/ 3p

(isothermal J-T coef.) (23)

5 h (

(3T\

= ^

P/H C p

(Joule-Thomson coef.). (24)
8

where except as indicated otherwise all dependent variables are functions

of p and T. Since Eqs 7 thru 18 are linear in A(p,T) each can be represented

as the sum of contributions from its three parts.

The computer print- out in appendix A is the FORTRAN 77 program for

calculating thermodynamic properties based on the equations contained in

this report. All calculations related to the same part of A(p,T) are

included in the same sub- routine. Thus the sub- routine "BASE" includes

Eqs (2-4) for the contribution from the base function to the Helmholtz

function, plus each of the contributions therefrom to Eqs 7 thru 18.

Similarly the sub- routine "QQ" includes all such contributions from the

residual function, and sub-routine "IDEAL" contains the contributions

from the ideal gas. Because the detailed equations are easily "read"

from the respective sub-routines they are not reproduced here.

The contributions to the various thermodynamics properties Eqs 7-24

are summed in sub-routine "THERM". Also included are routines for

properties associated with liquid-vapor coexistence. All calculations

are referred to the liquid at the triple point, for which state the

internal energy and entropy are zero. Though the units for the calculations

are p(g/cm 3 ), T(K) and joules, routine "UNIT" is included which allows

use of other units. Lastly, the iteration sub-routine "DFIND" allows

access with independent variables P and T.

All equations, including parameters and constants listed in sections

2 thru 4 are included in the program listing in Appendix A. Should any

inconsistencies be found between the text and Appendix A, it is likely

that the errors are in the text. Please let us know if you discover

any.
9

5. Discussion

The data selected for the derivation of the thermodynamic surface

were primarily PpT data. These cover the range from 0° <_ t <_ 1000°C and

0 < P < 1000 MPa. Other kinds of data used include values for the

enthalpy of the saturated liquid and for the isothermal compressibility

of the liquid below 100°C. No mathemtaical constraints were used to

impose exact accord with pre-assigned values. The surface yields values

consistent with all input data within reasonable estimates of their

tolerance. It has been further validated by extensive comparison with

high quality thermodynamic data, including other than PpT that had not

been used in the derivation of the surface.

At very high pressure and at very high temperatures thermodynamic

values calculated with the surface are in accord with those given by the

(theoretical) base function, so that the surface should yield useful

estimates well beyond the range of the data used in its derivation.

Thus, except very close to the melting curve, the surface should remain

quantitative at pressures equal to 2000 MPa and at least semi-quantitative

to 4000 MPa, for temperatures in the range 250 ± t <_ 2500 K for undissociated

fluid water. Lastly, because the surface is analytic, it yields some

results at the critical point and its immediate neighborhood which are

at variance with what modern theory predicts. Though we have not found

significant departure from any reliable critical region data, we caution

that the surface may not conform to all theoretical expectations in the

region defined by

T = T + 0.5 K

c —

p = P

c

(1 ± -2) ,

where T

c and p c are values at the critical point.
'

.

.

-
REFERENCES

1. L. Haar, J.G. Gallagher, and G.S. Kell, in preparation.

2. IUPAC Commission on Atomic Weights and Isotopic Aburdances, Pure

and Applied Chem. 52, 2349 (1980).

3. E.R. Cohen and B.N. Taylor, J. Phys., Chem. Ref. Data 2 , 663.

(1973).

4. L. Haar and S.H. Shenker, Proc. Fifth Symp. on Thermophys. Prop.

A.S.M.E. (1970)

5. L. Haar and S.H. Shenker, J. Chem. Phys. ^5 4951 (1971).

6. J.R. Wardell, G. Wilmot, L. Haar and M. Klein, Proc. 12th JANNAF

Combustion Meeting, Chem. Prop. Info. Agency, publ . 273 (1975).

7. L. Haar, G. Powell, M. Klein, J. Wardell, and G. Wilmot, Proc. 13th

JANNAF Combustion Meeting, Chem. Prop. Info. Agency, publ. 281

(1976).

8. G. Powell, G. Wilmot, L. Haar and M. Klein, Interior Ballistics of

Guns, M. Summerfield, Ed., Part IV, P. 307, Am Inst. Aero, and

Astro., New York (1979).

9. M. Klein, L. Haar, G. Powell, and G. Wilmot, 15th JANNAF Combustion

Meeting, Chem. Prop. Info. Agency, Publ. 297 (1978).

10. F. Kohler and L. Haar, accepted for publ. in J. Chem. Phys.

11. L. Haar, J. Gallagher, in preparation

12. L. Haar, J. Gallagher, and G. Kell, Proc., 9th Int. Conf. on Prop.

Steam, Munich (1979), J. Straub and K. Scheffler, Eds., P. 69,

Pergamon Press, Oxford (1980).

13. H.W. Woolley, Proc., 9th Int. Conf. on Prop. Steam, Munich (1979),

J. Straub and K. Scheffler, Eds., P. 166, Pergamon Press, Oxford

(1980).
'

'

It
APPENDIX A

This appendix contains the FORTRAN 77 programs which yield

values for thermodynamic properties for liquid and gaseous states

for water. The independent variables are temperature (T) and

density (D). The programs should yield useful results in the range

250 <_ T(K) < 2500

10“ 8

<_ D(g/cm3

) <_ D(Pmax),

where Pmax is the lesser of the pressure of melting ice or 4000 MPa.

For values of density less than or equal, to 10" 8 g/cm3

, the program

sets the density of 10“ 8 g/cm 3

.
AZ

Outlined below are the procedures for use of the subroutines which

follow. In.this outline, the symbol used t for pressure is "P’ , for

density is "B” , and for temperature is ”T". Calculations are made in the

internal units of the program which are MPa, g/cm3 and deg K respectively.

See III below for conversions to other units. Examples of the use of these

routines are given in three main programs in Appendix B, with samples of

the output which can be used for checking the implementation of the programs.

I. Jor Density and Temperature as independent variables, the sequence

CALL BB(T)

call base (d , t

)

CALL Q0(T ,D)

will place quantities in COMMON from which

P s D^GASCON^T^ZB + Q0

dP/dD = GASCON*.?* (ZB+T*DZB) + Q5

To obtain the other thermodynamic properties:

CALL THERM (D,T)

will olace in COMMON "AB“ (=A/RT), "GB" (=G/RT ) r "SB” (=S/R),

"UB" (=U/RT), "HB“ (=H/RT), ”C7B" (=Cv/R), ‘*CPD' £=Cp/R),

DPET' ( =dP/dT , MPa/K), "DVD?" (=d7/dT, cm3/gl) ,

‘ CJTT"

(isothermal joule-thompson coeff, cm3/g) and ’CJTH" (adiabatic

joule-thompson coeff, K/MPa). These properties in the desired

units can then be obtained by R or RT as appropriate in the

desired units .

II. Pressure and Temperature as independent variables

CALL DJIND( DOUT , F , LGUESS ,T,DPDD)

This will return the density corresponding to the input P and T as DOUT from which the procedures in I can be used for further calculations. This routine requires an initial approximate guess for the

density for use to begin the Newton iteration. dP/dD is also returned

as a byproduct.

III. Use of other systems of units.

A group of three subroutines is included which facilitate the use

of a system of units chosen by the user:

CALL UNIT

will query the user as to which units are desired for temperature,

pressure, density and energy, and set up the necessary parameters for

converting from the desired units to the internal units and baci

again for the output. The names of the units in alphanumeric

characters are also placed in COMMON for use in laoeling.

FUNCTION TTT ( Tex ternal ) , and

JUNCTION TTI ( Tin te rnal

)

will convert from external T to deg £, and from deg K to the external

T respectively.
O

O

O

o

O

A3

BLOCK DATA

implicit double precision (a-h, o-z

)

real p,q

COMMON /ACONST/

COMMON /NCONST/

COMMON /ELLCON/

COMMON /BCONST/

COMMON /ADDCON

/

THIS BLOCKDATA SUBROUTINE

USED IN THE REST OF THE

VM, GASCON, TZ,AA, ZB, DZB ,YB ,UREF ,SREF

G(40) ,11(40 ),JJ( 40) ,NC

G1,G2,GF,B1 ,B2,B1T ,B2T ,B1TT ,B2TT

P ( 10 ) , 0 ( 10

)

ATZ (4 ) , ADZ (4 ) , AAT (

) ,AAD(4)

SUPPLIES MOST OF THE FIXED PARAMETERS

ROUTINES. P IS THE b (i ) of TABLE I,

I , G1,G2 AND GF

G , 1 1 , JJ ARE THE

ARE THE ALPHA,

g ( i ) , k(i) AND

Q IS THE B ( i ) OF TABLE BETA

AND GAMMA OF EC 2, AND

1 ( i ) of EQ 5.

DATA ATZ/2*64.D1,641.6D0,27.D1/,ADZ/3*.319D0,1.55D0/,AAT/2*2.D4

1 ,4.D4,25 . D0/ , AAD/34 . D0 ,4 .Dl ,3 . Dl , 1 .05D3/

DATA WM/18.0152D0/, GASCON/. 461522D0/ , TZ/647 .073D0/ , AA/1 .D0/ , NC /36/

DATA UREF ,S REF/ -43 28 .455039D0 ,7 .6160802D0/

DATA Gl,G2,GF/ll.d0 ,44 .333333333333d0 , 3 . 5d0/

DATA P/.747862S ,-.3540782,2*0. , .007159876 ,0 . ,-.003528426,3*0./

DATA Q/l . 1278334,0 . ,-.5944001,-5 .010996,0. , .63684256 ,4*0 ./

DATA &/- . 52062S68529023D3 , .22744901424408D4 , .78779333020687D3

1,

-. 69830 5273749S4D2, . 17863832875422D5 ,-.39514731563338D5

2, .33803884280753D5 ,- . 13855050202703D5 , - .25637436613260D6

3.

.

4821257 598 1415D6 ,- . 34183016969660D6 , .122231564 17448 D6

4.

.

1 1797433655832D7,-. 217348 101 10373D7, . 1082S952168620D7

5

,

-. 25441998 064 049D6 ,-. 3 1377774 947767 D7 , .52911910757704D7

6

,

- . 13802577 177877D7 , - . 25 10991436900 1D6 , .465618261 15608D7

7, 72752773275387D7, .41774246148 294D6, . 1401 6358244614D7

6 ,- .31 555231392127D7 , . 47929666384584D7 , . 40912664 781209D6

9 ,- . 13626 36S388386D7 , .69625220862664D6 ,- .10 834900096447D7

A

,

- . 2272282740 1688D6 , .38365486000660D6 , .63833257944332D4

B

,

.21757245522644D5 ,-. 2662794482977 0D4,-. 707 30 4180820 74D

5

C

,

-. 225D0 , -1 .68D0 , .055D0 ,-93.0D0/

DATA 11/4*0 ,4*1 ,4*2, 4*3 ,4*4, 4* 5, 4*6, 4*8, 2*2 ,0,4, 3*2, 4/

DATA JJ/ 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5,

1,2, 3, 5, 7, 1,3*4, 0,2, 0,0/

END
o

o

o

A 4

c SUBROUTINE BB (T

)

THIS SUBROUTINE CALCULATES THE B'S OF EQS . 3,4 USING COEFFICIENTS

FROM BLOCXDATA, CALCULATING ALSO THE FIRST AND SECOND DERIVS

W.R.TO TEMP. THE B'S CALCULATED HERE ARE IN CM3/G

.

implicit double precision (a-h, o-z)

real

COMMON /ELLCON/ G1 ,G2 ,GF ,B1 ,E2 ,E1T ,E2T ,B1TT ,B2TT

COMMON /ACONST/ WM , GASCON , TZ , AA , Z ,DZ ,T ,UREF ,SREF

COMMON /BCONST/ P(10),Q(10)

DIMENSION V(10)

V(l)=l.

DO 2 1=2,10

2 V (I )=VII-1)*TZ/T

E 1=P( 1 )+F ( 2)*AL0G ( 1 . /V ( 2 )

)

52=Q(1)

E1T=P (2)*V(2) /TZ

E2T=0

.

B1TT=0 . B2TT=0

.

DO 4 1=3,10

B1=B1+P ( I )*V(I-1)

B2=B2+g(I )*V(I-1)

31T=B IT- ( 1-2 )*P ( I )*V ( 1-1 ) /T

B2T=B2T-( I-2)*0(I )*V(I-1)/T

B1TT=B1TT +P ( I

) * ( I -2 )**2*V ( 1-1 ) /T/T

4 E2TT=B2TT+Q(I )*(I-2 )**2*V ( 1-1 ) /T/T

B1TT=E1TT-E1T/T

B2TT=B2TT-B2T/T

RETURN

END
o

o

o

o

o

A S'

C

SUBROUTINE BASE (D ,T

)

C This function calculates Z [=Pbase/( DRT) ] , C and also Abase ,Gbase , Sbase ,Ubase ,H base ,CVbase , AND 1/(DRT) * DP/DT

C for the base fct (called DPDTB)

C The AB,GB ,SB ,U3,HB and CVB are calculated in dimensionless units:

C AB/RT, GB/RT, SB/R, etc.

IMPLICIT REALMS (A-H.O-Z)

COMMON /ELLCON/ G1 ,G2 f GF ,Bl ,B2 ,E IT ,B2T ,B1TT ,B2TT

G1,G2 AND GF ARE THE ALPHA, BETA AND GAMMA OF EQ 2, WHICH ARE

SUPPLIED BY THE BLOCKDATA ROUTINE. B1 AND B2 ARE THE ’‘EXCLUDED

VOLUME' AND “2ND VIRIAL” (EQS 3 AND 4) SUPPLIED EY THE SUBROUTINE

BB (T )

,

WHICH ALSO SUPPLIES THE 1ST AND 2ND DERIVATIVES WITH

RESPECT TO T ( E1T,B2T ,B1TT ,E2TT ) . COMMON /BASEF/ AB ,GB , SB , UB , HB , CVB , DPDTB

COMMON /ACONST/ WM , GASCON , TZ, A ,Z ,DZ ,Y ,UREF ,SREF

Y= . 25^31*0

X=1 .-Y

Z0=(1 .+G1*Y+G2*Y*Y)/X**3

Z=Z0+4.*Y*(32/B1-GF)

DZ0=(G1+2.*G2*Y)/X**3 + 3.*( l.+G 1*Y+G2*Y*Y) /X**4

DZ=DZ0+4.*(E2/E1-GF)

AB = -DLOG ( X

) - ( G2-1 . )/X+28 . 16666667D0/X/X+4 .^Y* ( B2/B1-GF

)

1 +15 . 166666667D0 + DLOG ( D*T*4 .55483D0

)

G3 = AB + Z BB2TT=T*T*E2TT

UB= -T*B1T*(Z-1.-D*B2)/B1-D*T*B2T

HE=Z+UB

C VB=2 .*UB+ ( Z0-1 .)*( (T*B1T/B1)**2-T*T*B1TT/B1)

1 - D*(BB2TT - GF^BITT^T^T ) -( T*B1T/B1 )**2*Y*DZ0

DPDTB = Z/T + D* ( D Z* B IT /4.+ 323-32/31*3 IT

)

SB = UB - AB

RETURN

END

C
o

o

o

oo

o

o

A 6>

SUBROUTI NS QQ(T ,D

)

THIS ROUTINE CALCULATES , FOR A GIVEN T(K ) AND D(G/CM3), THE RESIDUAL

CONTRIBUTIONS TO: PRESSURE (Q), HELMHOLTZ FCT (AR) , DP/DRHO (Q5),

AND ALSO TO THE GIBBS FUNCTION, ENTROPY , INTERNAL ENERGY, ENTHALPY,

ISOCEORIC HEAT CAFACITY AND DPDT . (EQ 5)

TERMS 37 THRU 39 ARE THE ADDITIONAL TERMS AFFECTING ONLY THE

IMMEDIATE VICINITY OF THE CRITICAL POINT, AND TERM 40 IS THE

ADDITIONAL TERM IMPROVING THE LOV T, HIGH P REGION.

IMPLICIT REAL*8(A-H,0-Z)

COMMON /RESF/ AR,GR ,SR ,UR ,HR, CVR ,DPDTR

COMMON /QQQQ/ Q,Q5

DIMENSION QR( 11) ,QT(10) ,0.ZR(9) ,QZT(9)

EQUIVALENCE (QR(3) ,QZR(l) ) , (QT(2) ,QZT(1)

)

COMMON /N CONST / G (40 ) , 1 1 ( 40 ) , J J( 40 ) ,

COMMON /ACONST/ WM , GASCON , TZ , AA , Z , DZ ,Y ,UREF , SREF

COMMON /ADDCON / ATZ (4 ) ,ADZ (4 ), AAT (4 ) , AAD (4

)

RT=GASCON*T

QR( 1

) =0

.

QS=0.

Q=0 ,D0

AR=0.D0

DADT=0 . C VR=0 . DPDTR=0

.

E=DEXP (-AA^D

)

O10=Da!c D* E

Q20=l .D0-E

QR(2)=Q10

V=TZ/T

QT ( 1

) =T/TZ

DO 4 1-2,10

QR(I+1)=QR(I )*Q20

4 QT(I)=QT(I-1)*V

DO 10 1=1,

K=I I ( I )+l

L=JJ ( I

)

Z z = K

QP=G ( I )* AA*QZ R ( K-l )*OZT (L ) Q=Q+QP

C5 = 05 + AA*(2./D-AA*(1.-E*(K-1)/Q20) )*QP

AR=AR + G( I )*QZR(K)*QZT(L) /Q10/ZZ/RT

DFDT=C20**K* ( 1-L ) *QZT ( L + 1 )/TZ /K

D2F=L*DFDT

DPT=DFDT*Q10*AA*K/Q20

DADT=DADT+G ( I )*DFDT

DPDTR=DPDTR+G ( I )*DPT

10 CVR=CVR-G(I )*D2F/GASC0N

CP=0.

Q2A=0 .
A 7

DO 20 J=37 ,40

If(G(J).EQ.0.D0) GO TO 20

K=I I ( J

)

KM= JJ { J

)

DDZ = AD Z (J-36

)

DEL - D/DDZ ~ 1 IE (DABS (DEL) .LT.l.D-10) DEL=1.D-10

DD = DEL-DEL

EX1 = -AAD( J-36)*DEL**K

DEX=DEXP(EX1 )*DEL**KM

ATT = AAT ( J-36

}

TX = ATZ( J-36 ' TAU = T/TX-1

.

EX2 = -ATT^TAU-TAU

TEX = DEXP ( EX2

)

CIO = TVS’ TTY QM = KM/ DEL - K*AAD (J-36 j^DEL** ( K-l

)

FCT=QM*D**2*Q10/DDZ

Q5T = FCT*(2./D+QM/DDZ) -(D/DDZ )**2*Q10* ( KM/DEL /DEL1 K* (K-l )-AAD ( J-36 j^DEL** ( K-2 ) ) 05 = Q5 + QST^G ( J

)

CP = OP + G ( J )*ECT

DADT = DADT - 2 .*G ( J )*ATT*TAU*O10/TX

DPDTR = DPDTR - 2 .*G ( J)*ATT*TAU*FCT/TX

02A = Q2A + T*G(J)*(4.*ATT*EX2+2.*ATT)*O10/TX/TX

AR = AR + 010-G ( J ) /RT

20 CONTINUE

SR=-DADT/GASCON

UR=AR+SR

C VR=C VR+ Q2 A/GAS CON

Q=0+QP

RETURN

END

C
o

o

o

o

A 8

SUBROUTINE DFIND (DOUT,P , D,T,DPD)

ROUTINE TO FIND DENSITY CORRESPONDING TO INPUT PRESSURE P(MPA), AND

TEMPERATURE T(K), USING INITIAL GUESS DENSITY D(G/CM3). THE OUTPUT

DENSITY IS IN G/CM3 , ALSO, AS A BYPRODUCT, DP/DRHO IS CALCULATED

( "DPD" , MPA CM3/G)

IMPLICIT REAL*8 ( A-H ,0-Z

)

COMMON /QCQQ/ Q0,Q

5

COMMON /ACONST/ WM, GASCON, TZ,AA,Z,DZ,Y, UR EF ,SREF

DD=D

RT=GASCON 5!‘T IF( DD . LE .0 . ) DD=1 .D-8

IF ( DD . GT . 1 . 9 ) DD=1 .

9

L=0

9 L=L+1

11 IF(DD.LE.0.) DD=l.D-8

IF( DD .GT . 1 .9 ) DD=1 .9

CALL BASE (DD , T

)

CALL QQ(T ,DD)

PP = RT*DD*Z + Q0

DPD=RT*(Z+Y*DZ)+Q5

C THE FOLLOWING 3 LINES CHECK FOR NEGATIVE DP/DRHO, AND IF SO ASSUME

C GUESS TO BE IN 2-PHASE REGION, AND CORRECT GUESS ACCORDINGLY.

IF( DPD .GT .0 ,D0 ) GO TO 13

IF(D.GE. .2967D0 ) DD=DD*1.02D0

IF(D.LT. .2967D0) DD=DD*.98D0

IF (L.LE. 10 ) GO TO 9

* 13 DPDX=DPDsic l . 1

IF ( DPDX . LT . . 1 ) DPDX= .

1

DP=DABS ( 1 .-PP/P

)

IF(DP .LT.l.D-e) GO TO 20

IF ( D . GT . .3 .AND. DP.LT.l.D-7) GO TO 20

I F ( D . GT . .7 .AND. DP.LT.l.D-6) GO TO 20

X= ( P-PP ) /DPDX

IF(DABS(X) .GT . .1) X=X* . 1/DABS (X

)

DD=DD+X

IF ( DD . LE . 0 . ) DD=l.D-8

19 IF ( L . LE. 30 ) GO TO 9

20 CONTINUE

DOUT=DD

RETURN

END
o

o

o

A9

c SUBROUTINE THERM(D,T)

IMPLICIT REAL*8 (A-H 0”Z

)

COMMON /ACONST/ WM , GASCON , TZ , AA , ZB ,DZB ,Y ,UREF ,SREF

COMMON /QOQQ/ QP , QDP

COMMON /EASEF/AB ,GB , SB ,UE ,HB , C VB ,DPDTB

COMMON /RESF/AR ,GR , SR, UR , HR , CVR , DPDTR

COMMON /IDF/ AI ,GI ,SI,UI,HI,CVI,CPI

COMMON /FCTS/AD ,GD , SD,UD ,HD , CVD , CPD ,DPDT ,DVDT , C JTT , CJTH

THIS SUBROUTINE CALCULATES THE THERMODYNAMIC FUNCTIONS IN

DIMENSIONLESS UNITS (AD=A/RT , GD=G/RT, SD=S/R, UD=U/RT,

HD=H/RT, CVD=C V/R , AND CPD=CP/R)

CALL IDEAL (T)

RT = GASCON*T

Z = ZB + QP/RT/D

DPDD = RT* ( ZB+Y*DZB ) + QDP

AD = AE + AR + AI - UREF/T + 5REF

GD = AD + Z UB = UB + UR + UI - UREF/T

DPDT = RT*D*DFDTB + DPDTR

CVD = CVB + CVR + CVI

CPD = CVD + T*DPDT**2/(D*D*DPDD*GASC0N)

HD = UD + Z SD = SB + SR + SI - SREF

DVDT=DPDT /DPDD /D/D

CJTT=1./D-T*DVDT

C JTH^-CJTT/CPD/GAS CON

RETURN

END

C

FUNCTION PS ( T

)

C This function calculates an approximation to the vapor pressure, PS,

C as a function of the input temperature. The vapor pressure

C calculated agrees with the vapor pressure predicted by the* surface

C to within .02? to within a degree or so of the critical temperature,

C and can serve as an initial guess for further refinement by

C imposing the condition that Gl=Gv. t IMPLICIT REAL*6 (A-H.O-Z)

DIMENSION A (8 ) /-7 .8889166D0 ,2 . 5514255D0 ,-6 .716169D0

1 ,33 .23949 5D0, -105 .38479D0 , 174 .35319D0 ,-148 .39348D0

2 , 48 .631602D0/

IF ( T . GT . 214 .D2 } GQ TO 2

PL=6 . 3573118D0-8858 .843D0/T+607 . 56335D0*T** (-.6)

PS=.1*D£XP(PL)

RETURN

2 V=T/ 647 . 25D0

W=DABS ( 1 .D0-V

!

B=0 .D0

DO 4 1=1,8

Z = I 4 B=B+A(I)*W**((Z+l.)/2.)

Q=B/V

PS=22 . 093D0*BEXP ( Q

)

RETURN

END
A

C

FUNCTION TSAT(P)

C This function calculates the saturation temperature for a given

C pressure, by an iterative process using PSAT and TDPSDT.

REAL*8 P , PS , TO, TSAT, TDPSDT

TSAT=0

.

IF (P .GT .22.05) RETURN

K=0

PL=2 . 3025£5*DLOG ( P

)

C PL=LOGE(10)*LOGE(P) TO CONVERT EQUATION FROM BARS TO MPa

TG=372 .8 3+?!* ( 27 .7589+PL* (2 .3819+PL* ( .24834 +PL* .0193855 ) )

)

1 IF (TG .LT . 273 . 15 ) TG=273.15

IF( TG .GT .647 . ) TG=647.

I F ( K . LT . 8 ) GO TO 2

WRITE (6, 3 ) K , P ,PP , TG

3 FORMAT ()

GO TO 8

2 K=K + 1

PP=PS ( TG

)

DP= TDPSDT (TG

)

IF(ABS(1.-PP/P) .LT. .00001) GO TO 8

TG = TG* ( 1 . + ( P-PP ) /DP

)

GO TO 1

8 TSAT=TG

RETURN

END

C FUNCTION TDPSDT (T

)

C This function calculates T*(dPs/dT), and is used by the function TSAT

IMPLICIT REAL*8 (A-H,0-Z)

DIMENSION A( 8 ) /-7 .888S166D0 , 2 . 5514255D0 , -6. 716169D0

1 ,33 .239495D0 ,-105 .38479D0 ,174.35319D0 ,-148 . 39348D0

2,48 .631602D0/

V =T/ 647. 25

W=1 .-V

B=0 . C=0 . DO 4 1=1,8

Z=I

Y=A(I }*¥**( (Z+l.J/2.

)

C=C+Y/W*( .5-.5*Z-l./V)

4 B=B+Y

Q=B/V

TDPSDT =2 2 . 093*DEXP (Q )*C

RETURN

END
A n

c SUBROUTINE IDEAL(T)

C THIS SUBROUTINE CALCULATES THE THERMODYNAMIC PROPERTIES FOR

C WATER IN THE IDEAL GAS STATE OF H .W . WOOLLEY

IMPLICIT REALMS (A-H.O-Z)

COMMON /IDF/ Al ,GI ,SI,UI,HI,CVI,CPI

DIMENSION C( 18)/.19720271018D2 f . 209662681977D2 483429455355D0

1 .

605743 189245D1 ,22 .560 2388 5D0 ,-9 .87532442D0 ,- .43135536513D1

2.

458155781D2 , - .47754901883D-1 , . 41238460633D-2 27929052352D-3

3 .

1448 159526 1D -4 , - . 56473658748D-6 , .16200446D-7,-.3303822796D-9

4.

451516067368D-11 . 370734122706D-13 , . 137546068238D-15/

C NOTE THAT THE TEMPERATURES ARE SCALED BY A FACTOR OF 100 HERE SO THAT

C THE EXPONENT OF THE COEFFICIENT OF THE LAST TERM WILL BE WITHIN

C RANGE FOR MOST COMPUTERS.

TT=T/1 .D2

TL=DLOG( TT

GI=-(C(1)/TT+C(2) )*TL

HI=(C(2)+C(l)*(l.D0-TL)/TT)

CPI =C ( 2

) -C ( 1 ) /TT

DO 8 1=3,16

GI=GI-C(I )*TT**(I-6'>

HI = HI+C(I )*(I-6)*TT**(I-6)

8 CPI=CPI+C(I)*(I-6)*(I-5)*TT**(I-6)

AI=GI-1.

UI=HI -1

CVI=CPI-1.

SI=UI-AI

RETURN

END

C

SUBROUTINE SECVIR ( T , VIR ' C THIS SUBROUTINE CALCULATES THE SECOND VIRIAL IN CM3/G

C AT TEMPERATURE T IN K.

IMPLICIT REAL*8 (A-H.O-Z)

COMMON /N CON ST / G (40 } , 1 1 (40 ), J J ( 40 ) , NC

COMMON /ELLCON/ G 1 ,G2 , GF , BB1 , BB2 ,B IT ,B2T ,B1 TT , B2TT

COMMON /OQCQ/ 02,05

COMMON /ACONST/ WM .GASCON , TC ,AA , Z , DZ ,Y ,UREF ,SREF

CALL BB( T

V = T C / T

V IR=BB2

DO 3 J=1 ,NC

IF ( 1 1 ( J ) . ME .0 } GO TO 3

L=JJ ( J

VIR=VIR+G( J)*V**(L-1 '/T /GASCON

3 CONTINUE

RETURN

END
o

o

o

o

c SUBROUTINE CORR (T ,P ,DL,DV,DELG ) C SUBROUTINE CORR WILL CALCULATE, FOR AN INPUT T AND P AT OR

C VAPOR PRESSURE, THE CORRESPONDING LIQUID AND VAPOR DENSIT

C DELG = ( GL-GV ) /RT FOR USE IN CALCULATING THE CORRECTION T

C PRESSURE FOR DELG = 0.

IMPLICIT DOUBLE PRECIS ION (A-H , O-Z

COMMON /QQQQ/ Q00,Q11

COMMON /ACONST/ WM , GASCON , TZ , AA , ZB ,D ZB ,TB ,UREF , SREF

COMMON /FCTS / AD , GD , SD , UD , HD , C VD , CPD , DPDT , D VDT , C JTT , CJ

DLIQ=DL

I F( DL . LE .0 . ) DLIQ = 1 . 11-.0004*T

CALL BB( T

RT=GASCON*T

CALL DFI ND (DL ,P ,DLIQ ,T , DQ)

CALL THERM (DL.T)

GL=GD

DVAP=DV

IF( DV .LE .0 . ) DVAP=P/GAS CON /T

CALL DFI ND(DV,P,DVAP,T,DQ)

IF(DV.LT.5.D-7) DV=5.D-7

CALL THERM (DV ,T)

GY=GD

DELG - GL-GV

RETURN

END

CC

SUBROUTINE PCORR ( T ,P ,DL ,DV

C SUBROUTINE PCORR WILL CALCULATE THE VAPOR PRESSURE P AND T

C VAPOR DENSITIES CORRESPONDING TO THE INPUT T, CORRECTED S

C GL-GV=0 . THE FUNCTION PS IS REQUIRED WHICH WILL GIVE A R

C GOOD APPROXIMATION TO THE VAPOR PRESSURE TO BE USED AS TH

C POINT FOR THE ITERATION.

IMPLICIT DOUBLE PRECISION (A-H, O-Z)

COMMON /ACONST/ WM , GASCON , TZ , AA , ZB ,DZB ,YB ,URFF , SREF

P = PS(T)

2 CALL CORR (T,P,BL,DV,DELG)

DP=DELG*GASCON*T/(l./DV-l./DL)

P = P+DP

IF ( DABS ( DELG ) .LT.l.D-4) RETURN

GOTO 2

END

CC

SUBROUTINE TCORR( T , P , DL , DV

SUBROUTINE TCORR IS SIMILAR TO ''PCORR" EXCEPT THAT THE TEM

CORRESPONDING TO THE INPUT VAPOR PRESSURE IS FOUND. FUNCT

ARE TSAT AND TDPSDT WHICH GIVE AN APPROXIMATION TO T ( SAT

T*DP(SAT)/DT

IMPLICIT DOUBLE PRECISION (A-H, O-Z)

COMMON /ACONST/ WM , GASCON , TZ , AA , ZB ,DZE , IB ,UREF , SREF

T = TSAT(P)

2 CALL CORR(T,P,DL,DV ,DELG)

DP= DELG5* GAS CON^T/tl./DV-l./DL)

T = T* (1 . -DP/TDPSDT (T } ) I F ( DA BS ( DELG ) . LT . 1 .D-4 ) RETURN

GO TO 2

END
cc The following 3 subroutines form a package allowing the us

c operate in a s/stem of units of his choice.

SUBROUTINE UNIT

C THIS SUBROUTINE QUERIES THE USER AS TO HIS CHOICE OF UNITS

C INTERNAL PARAMETERS APPROPRIATELY. THE INTERNAL UNITS OF

C TEMPERATURES IN K , DENSITIES IN G/CM3 , ALL OTHER QUANTITI

C IN DIMENSIONLESS UNITS AND DIMENSIONED AT OUTPUT TIME.

IMPLICIT DOUBLE PRECISION (A-H,0-Z)

DOUBLE PRECISION NT , ND , NP , NH, NNT , NND , NNP , NNH

COMMON /UNITS/ IT , ID , IP , IH , NT , ND . NP , NH ,FT ,FD ,FP , FH

DIMENSION FFD ( 4 ) , FFP ( 5 ) , FFH (6) , N NT ( 4 )

, NND (4 ) , NNP ( 5) r NN

DATA FFD/l.D-3, 1.D0, .0180152D0, .016018D0/

DATA FFF/1.D0 , 10.D0 ,9 .869232667D0 , 145 . 038D0 ,10 .1971D0/

DATA FFH/2*1.D0,18.0152D0, . 23884590D0 ,4 .30285666D0 , .42

DATA NNT/1HK , IH C , 1HR , 1HF/

DATA NND/6Hkg/m3 ,6Hg/cm3 ,6Hmol/L ,6Hlb/ft3/

DATA NNP/6H MPa ,6H Bar , 6H Atm ,6H PS I ,6Hkg/cm2/

DATA NNH/6hkJ/kft ,6H J/g f 6HJ/'mol ,6Hcal/g ,7Hcal/mol

DATA Al, A2 , A3 .A4/8HTEMPERAT ,7HDEN5ITY

1 ,8HPRESSURE , 8HENERGY /

WRITE ( 6,11) Al

30 WRITE (6,12)

READ ( 5 ,10*f END=99} IT

IF( IT . EQ .0 ) STOP

IF( IT .GT . 4 ) GOTO 30

NT=NNT(IT

)

WRI TE ( 6 , 11 ) A2

31 WRITE (6,13)

READ ( 5 ,10 , END=99 ) ID

IF ( ID . GT . 4 .OR. ID.LT.l) GOTO 31

ND=NND( I D

)

FD=FFD ( ID

)

WRITE(6, 11 ) A3

32 WRITE (6, 14)

READ ( 5 ,10 ,END=S9) I?

IF ( IP .GT .5 .OR. IP.LT.l) GOTO 32

N P= N N F ( I P ) FP=FFF (IP)

WRITE (6,11) A

4

33 WRITE(6,15)

READ( 5,10 , END =99) IH

IF ( IH .GT . 6 .OR. IH.LT.l) GOTO 33

NH=NNH (IH

)

FH=FFH ( I H

)

99

10

11

12

13

14

15

RETURN

STOP

FORMAT (

)

FORMAT ( ' FORMAT (

FORMAT (

FORMAT (

FORMAT ( 1RIES/MOL

END

ENTER UNITS

CHOOSE FROM

CHOOSE FROM

CHOOSE FROM

CHOOSE FROM

6=BTU/LE ')

CHOSEN FOR ', A8

)

1=DEG K, 2= DEG C, 3=DEG R, 4=DEG

1 =EG/M3, 2=G/CM3 , 3=M0L/L , 4=LE/F

1 = MPA , 2=BAR , 3=ATM , 4=PSIA, 5=KG

1=KJ/KG, 2=J/G, 3=J /MOL , 4=CALORI
A'*

FUNCTION TTT(T)

C FUNCTION TO CONVERT INPUT TEMPERATURES IN EXTERNAL UNITS T

DOUBLE PRECISION T , TTT ,FT »FD , FP , FH , NT , ND , NP ,NH

COMMON /UNITS/ IT , ID , IP , IH , NT f ND , NP , NH ,FT ,FD ,FP , FH

GO TO (1,2, 3, 4) ,IT

1 TTT=T

FT-1

RETURN

2 TTT=T -*-273 . 15D0

FT=1

RETURN

3 TTT=T/1.8D0

FT= *55555 555 555 56D0

RETURN

4 TTT=( T+459 .67D0 )/l .8D0

FT= . 5555555555556D0

RETURN

END

FUNCTION TTI(T)

C FUNCTION TO CONVERT INTERNAL TEMPERATURES IN DEG K TO EXTE

DOUBLE PRECISION T ,TTI , FT , FD ,FP ,FH , NT ,ND , NP ,NH

COMMON /UNITS/ IT , ID , IP , IH , NT , ND , NP , NH ,FT ,FD , FP , FH

GO TO (5, 6, 7, 8) ,IT

5 TTI =T

RETURN

6 TTI =T-273 . 15D0

RETURN

7 TTI=T*1.8D0

RETURN

8 TTI=T*1.6D0-45S.67D0

RETURN

END
APPENDIX B

This appendix contains sample calculations to help the user test

his program. Also listed are the main programs used to print out

numerical values.
C THIS IS A MAIN PROGRAM FOR THE CALCULATION OF TABLES OF PROPERTIES

C USING A CHOICE OF UNITS, AND AT A CHOICE OF CONSTANT TEMPERATURE,

C PRESSURE OR DENSITY. THE USER IS QUERIED AS TO HIS CHOICES OF

C UNITS AND VARIABLES AT EXECUTION TIME.

IMPLICIT DOUILE PRECISION ( A-H , O-Z)

DOUBLE PRECISION NT ,ND , NP , NH

COMMON /UNITS/ IT , ID , IP , IH , NT , ND , NP , NH ,FT ,FD ,FP ,FH

COMMON /QQQQ/ Q0,Q5

COMMON /FCTS/ AD , GD ,SD , UD ,HD , C VD , CPD ,DPDT

COMMON /ACONST/ WM, GASCON, TC,AA,Z,DZ,Y .UREF ,SREF

COMMON /N CONST / G (40 ) , 1 1 (40 ) , J J ( 40 ) ,NC

DATA NS1 ,NS2/2H m,2Hf t/

CALL UNIT

NS=NS1

IF ( ID .EQ .4 > NS = NS2

15 WRITE (6, 11

READ ( 5 ,* ,£ND=9 ) IOPT.XISO

IF ( IOPT . EQ .0 ) GO TO 9

GO TO (101,201,301) , IOPT

101 WRITE(6, 102 ) 102 FORMAT ( ' ENTER OTHER INDEPENDENT VARIABLE (2 FOR P, 3 FOR DENS'/

1' FOLLOWED BY INITIAL, FINAL AND INCREMENTAL VALUES OF THIS VAR.')

READ(5,*,END=9) JOPT, Y1 , Y2 , YI

IF(JOPT-i) 15,101,103

103 TT=XI SO

T=TTT (TT ) WRITE (6, 444) NT ,NP , ND , NP , NT ,NH ,NT , NH ,NS

IF( JOPT.EC .2) DGSS=Yl/FP/T/.4

IZ=0

CALL 3B( T ) P3S=20000

DVV=0 . IF(T.LT.TC) CALL PCORR( T,PSS, DLL ,DVV

DGSS=DVV

IF( DGSS. EQ .0 . ) DGSS=1.11-.0004*T

PSAT=PSS*FP

I F ( JOPT . EQ .2 .AND. Yl.GT.PSAT) IZ=3

IF(Yl.GT.PSAT) DGSS=DLL

IF ( JOPT. GE .3 ) I Z=3

FI N=Y 1-Y

DIN=PIN

PINC=YI/FP

DINC=YI*FD

22 IF( JOPT.EQ .2) PIN=PIN+YI

II( JOPT. EQ. 3) DIN-DI N+YI

I F( JOPT . EQ .2 .AND. PIN.GT.Y2) GO TO 101

IF( JOPT. EQ .3 .AND. DIN.GT.Y2) GO TO 101

I F (JOPT . EQ .2 ) PRES =P IN/FP

I F ( JOPT. EC. 3) D-DIN/FD

24 CONTINUE
23

26

27

IF( JOPT.EQ .2

IKJOPT.EQ.1

TSAYE=TT-YI

IF ( JOPT . EQ . 1

PSAVE=PI N -YI

IF( JOPT.EQ. 2

IF( JOPT. EQ .1

IZ=IZ+1

.OR. (JOPT.EQ. 2 .AND. PIN .LT .PSAT ) ) GO TO 26

.AND. T.LT.TS) GO TO 26

.AND. I0PT.EQ.2 .AND. IZ.LE.2) TT=TT I (TS

)

.AND. PRES .GT.PSAT/FP .AND. IZ.GE.2) GO TO 26

.AND. T.GT.TS .AND. IZ.GE.2) GO TO 26

IF( JOPT.EQ. 2) PRES =PSAT /FP

IF ( JOPT.EQ. 1) T=TS

IF ( I Z . EQ . 1 .AND. JOPT.EQ.l) DG3S=DLL

IF(IZ.EQ.l .AND. JOPT.EQ. 2) DG3S=DYY

IF( IZ .EQ .2 .AND. JOPT.EQ. 2) DGSS=DLL

IF ( I Z .EQ . 2 .AND. JOPT.EQ.l) DGSS=DYY

C AL L E B ( T ) IF ( IOPT . NS .3 .AND. J0PT.NE.3) CALL DFI ND (D ,PRES , DGSS , T , DQ

)

CALL QQ( T ,D ) CALL BA5E(D ,T

)

RT=GASCON*T

PDUM = D*RT*Z + Q0

IF ( JOPT . EQ .3 .OR. I0PT.EQ.3) PRES=PDUM

I F ( IOPT . EQ .3 .OR. JOPT.EQ. 3) DQ=RT*( Z+Y*D Z

) +Q5

DGSS=D + PINC/DQ

CALL THERM(D ,T

)

U = UD*T*GASCON*FH

C=DSQRT( DABS ( CPD*DQ*1 . D3/CYD )

)

IF( ID .EQ .4 ) C=C*3 . 280833

H = HD^T^GASCON^FH

S = SD*GASCON*FH*FT

dpdti=dpdt :J: f p*f t d'odd=d%*f d

COMP = 1 .E3/D/DQ/FP

DDDTL=1.D3*DPDT/D/DQ

CP=CPD*GASCON*FH*FT

CV=CVD*GASCON*FH*FT

VL=fd/D

DOUT=l./YL

’ POUT=PRES’!‘FP

WRITE (6, 21) TT,POUT,DOUT,DPDTX t DFDD,CY ,CP f S ,H,U,C

21 FORMAT (F9 .3 ,F12 .5 , F12. 7 t FH • 5, FI 1 .3 , 5F12 .4 , Fll .3)

444 FORMAT ( 5X , 1ST , 12X , 1HP , 10X , 1HD ,8X , 5HdP/dT , 5X , 5HdP/dD , 10X , 2HCr

l,10X f 2HCp,llX f 1HS ,9X , 1HH , 11X , 1HU ,7X f 7HYel Snd/3X,4Hdeg ,A1,7X,A8

2 , 4X A6 , 5X ,A3 1H/ , A1 ,4X 7H ,15X,4H ,A6,A1,4H

3 , 13X, 4H A6 ,4H - - ,8X ,A2 , 4H/sec/

)

11 FORMAT ( ' DO YOU WISH TO CALCULATE AN ISOTHERM (ENTER 1), AN ISOBA

1R (ENTER 2) OR AN ISOCHORE? (ENTER 3)V' FOLLOWED BY YALUS OF ISO.

2 (ENTER 0 TO DISCONTINUE)')

IF(IZ.EQ.l) WRI TE ( 6 , 12

)

12 FORMATd

)
IF(IZ.EQ.l) GO TO 22

IF(IZ.EQ.2 .AND. J0PT.EQ.2) PIN=FSAVE

IF(IZ.EQ.2 .AND. JOPT.EQ.l) TT=TSAVE

IE ( IZ . EQ . 2 ) IZ = 3

GO TO (22,204,304 ) ,IOPT

201 JOPT=l

PRES = XISO/F?

202 WR I TE ( 6 203)

203 FORMAT ( ' ENTER FIRST LAST AND INCREMENT OF T '

READ( 5 ,END=S ) T1,T2,YI

IF(T1.EQ.0.) GO TO 15

WRITE (6, 444) NT ,NP ,ND , NP ,NT ,NH , NT , NH ,NS

TT=T1-YI

T=TTT (TT

CALL TCORR( TS ,PRES f DLL ,DVV

D = DLL

IF( (T+YI*FT) .GT .TS ) D=DVV

I Z=3

I F ( (T+YI *FT ) .LT .TS )IZ=0

204 TT=TT+YI

T=TTT (TT

IF(TT.GT.T2) GO TO 202

CALL BB(T)

DGS S=D

GO TO 24

201 JOPT=l

D=XI SO

302 WRITE (6, 203

READ ( 5 , END=9 ) T1,T2,YI

IF ( T1 .LE .0 . ) GO TO 15

WRITE (6, 444) NT ,NP ,ND , NP , NT , NH ,NT , NH ,NS

TT=T1-YI

I Z=3

T=TTT( TT)

204 TT=TT+YI

T=TTT ( TT

IF( TT . GT . T2 ) GO TO 302

CALL BB(T)

GO TO 27

9 STOP

END
o

o

o

o

a

o THE FOLLOWING IS A MAIN PROGRAM FOR CALCULATING THE PROPERTIES

ALONG THE SATURATION CURVE UP TO 646.3 I. THE INPUTS ARE: INITIAL T,

FINAL T AND THE T INCREMENT (ALL IN K). THE VAPOR PRESSURE IS FIRST

CALCULATED USING "PS" AND ADJUSTED FROM THERE SO THAT

GL-GV=0. ALSO IN THE INPUT ARE AN INITIAL GUESS FOR THE LIQUID AND

VAPOR DENSITIES. IF NOT GIVEN, THEY WILL BE INTERNALLY CHOSEN.

IMPLICIT REAL*8 (A-H.O-Z)

COMMON /ACONST/ WM .GASCON, TZ ,AA, ZB ,DZB ,YB

COMMON /FCTS/ AD , GD ,SD ,UD ,HD , CVD ,CPD ,DPDT

3 READ ( 5 , 4, END =99 ) TI ,TF,TD,DL ,DG

4 FORMAT ()

WRITE ( 6,22)

T=TI-TD

6 T=T+TD

IF(T.GT.TF) GO TO 3

IF (T ,GT .646.3D0) GO TO 3

I JK=0

CALL BB (T

)

P=PS ( T) +1 JK*PT

C THE FOLLOWING TWO LINES CALCULATE AN INITIAL GUESS FOR THE DENSITIES

C IF GUESSES NOT SUPPLIED IN THE INPUT.

IF (DL . LE.0 . ) DL=1 . 11- . 0004*T

IF(DG .LE.0 . ) DG=P/GASCON/T

5 P=P+IJK*PT

RT=GASCON*T

CALL DFIND (RL ,P , DL , T, DPDL)

DL=RL

CALL THERM (RL ,T

)

SL=SD*GASCON

HL = HD^GASCON^T

VL=1./DL

GL=GD

CALL DFIND (RG ,P ,DG , T, DPDG)

DG=RG

IF (RG .LT .5 .D-7 ) RG=5.D-7

CALL THERM ( RG ,T ) SG=SD*GASCON

EG = HD*GASCON*T

VG=1./DG

GV=GD

HEAT=HG-HL

DELG=GL-GV

I F ( DABS (DELG ) . LT .2 .D-6 ) GOTO 15

PT = DELG^GASCON^T/ ( VG-VL

)

IF(T .GT .640.D0 ) PT=.1*PT

I JE=I JK+1

IF(IJK.EQ.l .OR. (T.GT.64.D1 .AND. IJI.LE.10)) GO TO 5

15 HEAT=HG-HL

WRITE (6,21) T , P ,RL , RG ,HL ,HG f HEAT ,SL ,SG,VL,VG

21 FORMAT ( F9.3,F13 .6, F9.6 ,F9.7 ,3FS .2 ,2F9 .4,F9 .5,F11 .3)

22 FORMAT ( ' T,K P.MPA DL ,G/CC DV.G/CC HL,J/G

1 LAT HT SL.J/GK 3V,J/G5 VL VG ' /)

GO TO 6

99 STOP

END
Ef>

C THIS IS A SAMPLE MAIN PROGRAM WHICH WILL SERVE AS AN EXAMPLE FOR

C THE USE OF THE SUBROUTINES AND FUNCTIONS GIVEN ABOVE, AND WHICH

C WILL PRINT OUT VALUES OF VARIOUS PROPERTIES CALCULATED FOR A

C GIVEN INPUT POINT TO A LARGE NUMBER OF SIGNIFICANT FIGURES,

C SUITABLE FOR USE AS A CHECK ON THE OPERATION OF THE PROGRAM.

C THE USER IS QUERIED AS TO THE UNITS DESIRED.

IMPLICIT DOUBLE PRECISION (A-H ,0-Z

DOUBLE PRECISION NT ,ND , NP , NH

COMMON /UNITS/ IT , ID , IP , IH , NT , ND , NP , NH ,FT ,FD , FP , FH

COMMON /QQQQ/ Q0,Q

COMMON /FCTS/ AD ,GD ,SD ,UD , HD , C VD , CPD ,DPDT ,DVDT , C JTT , CJTH

COMMON /ACONST/ WM , GASCON , TZ ,AA, Z ,DZ ,Y ,UREF ,SREF

COMMON /NCONST/ G ( 40 ) , 1 1 (40 ) , J J ( 40 )

, NC

DATA NSS1,NSS2/2H m,2Hft/

CALL UNIT

N S=NS S

IF( ID . EC .4 ) N S=NSS

15 WRITE(6, 11)

100 READ( 5,* ,END=9) IOPT,X,TT

T=TTT( TT

RT=GAS CO N*T

CALL EB( T

IF( IOPT.LE .0 ) GOTO 9

GOTO (101, 102, 100, 100, 100), IOPT

101 DD=X

D=DD*FD

CALL CO ( T ,D

CALL BASE (D ,T

PRES = FP* ( RT*D*Z + 00)

DQ=RT* ( Z+Y*DZ ) +05

GOTO 111

102 PRES=X

P=PRES/FP

DGS S=P/T / .4

PSAT=20000.

IF(T.LT.TZ) CALL PCORR (T , PS AT, DLL ,DGSS

I F ( P . GT . PSAT ) DGSS=DLL

IF( P. GE. PSAT .AND. T.LT.523.15) DGSS=1 . 11-. 0004*T

CALL DFI ND(D,P,DGSS,T,DQ)

DD=D/FD

111 CALL THERM ( D , T

U = UD*T*GAS CON*FH

C=DSORT ( DABS ( CPD*D0*1 . D3/CVD )

IF ( ID . EO .4 ) C=C*3. 280833

H = HD*T*GASCON*FH

S = SD*GASCON*FH*FT

CP=C?D*GASCON*FH*FT

C V=C VD*G ASCON*FH*FT

VL=1./D

DPDD = D0’!t FD*FP

DPDT1=DPDT*FP*FT

WRI TE ( 6 , 20 ) TT ,NT ,PRES , NP ,DD, ND j
£7

WRITE(e,21) DPDTl,DPDD,CT,NH,NT,CP,S,H,NH,n,C,NS,CJTT,CJTH,DVDT

20 FORMAT( ' T

= ,F12.4, deg ,A1 ,5X» P = »F13 .6, IX » A6, 5X * D = 1,F14.10,1X,A8)

21 FORMAT ( ' DP/DT = ' ,F16 .9 ,6X

, 'DP/DD =',F16.5/

2

' C V =',F12.6,1X,A6,A1,5X,'CP = ' ,F12 .6 ,6X

, 'S =',F12.6/

3

' H =',F14.6,1X,A6,5X, 'U = ' ,F14 .6 ,6X , 'VEL SND = ' , F14. 6 ,A2

, '/sec '/

4

' JT( T) =' ,F11 .5 ,5X, 'JT(H) = '

, F 11 .5 ,5X

, 'DV/DT =',F12.6/)

11 FORMAT ( ' INTER OPTION, X, AND T, WHERE FOR OPTION=l, X=DENSITY'/

1' AND FOR 0PTI0N=2 , X=PRESSURE (ENTER 0 TO QUIT)') . GO TO 100

9 STOP

END

jXQT POINT

ENTER UNITS CHOSEN FOR TEMPERAT

CHOOSE FROM 1=DEG X, 2=DEG C, 3=D£G R, 4=DEG F

>1

ENTER UNITS CHOSEN FOR DENSITY

CHOOSI FROM 1=KG/M3 , 2=G/CM3, 3=M0L/L, 4=LB/FT3

>2

ENTER UNITS CHOSEN FOR PRESSURE • CHOOSE FROM 1=MPA , 2=BAR , 3=ATM, 4=PSIA, 5=XG/CM2

>2

ENTER UNITS CHOSEN FOR ENERGY

CHOOSE FROM 1=X J/XG , 2=J/G, 3=J/M0L, 4=CALORIES /G , 5=CAL0RI ES/MOL , 6=BTU/L3

>2

ENTER OPTION, X, AND T, WHERE FOR OPTION=l, X=DENSITY

AND FOR OPTI 0N=2 , X=PRESSURE (ENTER 2 TO QUIT)

>1, .9, 873.15

T = 673.1500 deg X P = 7110.605028 Bar D = .9000000000 g/cm3

DP/DT = 14.491600834 DP/DD = 28719.49752

C7 = 2.827220 J/g X CP = 3.615462 5 = 4.064690

H = 2779.151751 J/g U = 1989.062303 VEL SND = 1916.419293 m/sec

JT( T ) = .56718 JT(H) = -.15688 DV/DT = .000623 -

>2 , 225. , 648.15

T = 648.1500 deg K

DP/DT = 3.447886683

CV = 3.743787 J/g X

H = 1965.692198 J/g

P = 225.000000 Bar D = .4103745556 g/cm3

DP/DD = 63.95378

CP = 75.284775 S = 4.221788

U = 1910.664237 VEL SND = 358.617190 m/sec

JT ( T

)

= -205.05536 JT(H) = 2.72373 DV/DT = .320130

>2, .00617, 273.16

T = 273.1600 deg K

DP/DT = -1.576872063

CV = 4.225225 J/g K

H = .000617 J/g

JT ( T ) =

P = .006170 Bar D = .9997782169 g/cm3

DP/DD = 19608 .40085

CP = 4.228690 S = .000000

U = -.000000 VEL SND = 1400.674132 n/sec

1 .02220 JT (H ) = -.24173 DV/DT = - .000080

>
T P D V

deg C MPa g/cm3 cm3/g

50.000 .0100 . 000067 14869.238335

50.000 .0123 .000083 12036.663719

50.000 .0123 .987991 1.012155

50.000 .1000 .988030 1 .012116

50.000 1 .0000 .988422 1.011713

50 .000 10 .0000 .992305 1.007754

50.000 20 . 0000 .996526 1.003486

50 .000 30.0000 1 . 000656 .999344

50.000 40 .0000 1.004700 .995322

50.000 50 .0000 1.008663 .991411

50.000 60 . 0000 1.012550 .987605

50.000 70 .0000 1.016365 .963899

50 .000 80.0000 1.020109 .980287

50.000 90 . 0000 1.023788 .976765

50.000 100 .0000 1.027403 .973328

50.000 200 . 0000 1.060447 .942999

50.000 300.0000 1.088794 .918446

50 . 000 400 . 0000 1.113570 .898013

50 .000 500.0000 1.135810 .880429

50 .000 600 . 0000 1.156346 .864793

50 .000 700 .0000 1.175655 .850590

50.000 800 .0000 1.193766 .837685

50 .000 900 .0000 1.210356 .826203

50 .000 1000 .0000 1.225100 .816260

T P D V

deg C MPa ta / c m3 cm3/g

250.000 .0100 . 000041 24136.161956

250.000 .1000 . 000416 2406.053164

250 .000 1 .0000 . 004298 232.644926

250.000 3.9736 .019956 £0 .110581

250 .000 3 .9736 .799072 1.251452

250 .000 10 .0000 .805899 1.240850

250 . 000 20 .0000 .816284 1.225064

250 .000 30 . 0000 .825733 1.211046

250 .000 40 .0000 .834438 1.198411

250 .000 50 .0000 .842537 1.186892

250 .000 60 .0000 .350128 1.176293

250 .000 70 .0000 .857288 1.166469

250 .000 80 .0000 .864075 1.157307

250.000 90 .0000 .870537 1.148717

250.000 100 .0000 .876711 1.140627

252 . 000 200 .0000 .927623 1 .078024

250 .0 00 300 .0000 .966627 1.034311

250.002 400 .0000 .999498 1 .000502

250.000 500 .0000 1.027884 .972872

250.000 600 .0000 1 .053194 .949493

250.000 700 .0000 1.076161 .929229

250 . 000 800 .0000 1.097271 .911352

250.000 900 .0000 1.116861 .895367

250 . 020 1000 .0000 1 . 135179 .880918
T P D V

deg C MPa g /cm3 cm3/g

375.000 .0100 .000033 29909.055417

375.000 .1000 .000335 2986.855135

375.000 1 .0000 .003395 294.567937

375.000 10 . 0000 .040763 24.532263

375.000 20 . 0000 .130420 7.667535

375.000 22 .3292 .320009 3.124909

375.000 30 . 0000 .558254 1.791298

375.000 40 . 0000 .609556 1.640537

375.000 50.0000 .641323 1.559276

375.000 50 .0000 .665164 1.503388

375.000 70 .0000 .664567 1.460777

375.000 80 . 0000 .701090 1.426350

375.000 90.0000 .715576 1 .397475

375.000 100 . 0000 .728535 1.372617

375 . 000 200 .0000 .816920 1.224110

375.000 300 .0000 .873641 1.144634

375.000 400 .0000 .917244 1.090222

375.000 500 .0000 .953396 1 .048882

375.000 600 . 0000 .984657 1.0155e2

375.000 700 .0000 e 1.012420 .987732

375.000 800 . 0000 1.037537 .963821

375.000 900 .0000 1.060570 .942889

375.000 1000 .0000 1.081910 .924291

T P D V

deg C MPa g /cm3 cm3/g

500.000 .0100 .000028 35679.867768

500.000 .1000 .00 028

3565.549685

500 . 000 1.0000 .002824 354.099874

500 .000 10.0000 . 030503 32.783626

500 .000 20 .0000 . 067711 14.768678

500 .000 30 .0000 .115259 8.676122

500.000 40.0000 . 177974 5.618801

500 . 000 50 .0000 .256947 3.891860

500 .000 60 . 0000 .338443 2.954707

50 0 . 000 70 .0000 .405762 2.464496

500 .000 80 .0000 .456993 2.188216

500 . 000 90 .0000 .496532 2.013967

500.000 100 .0000 .528211 1 .893181

500 .000 200 .0000 .691177 1.446808

500 .000 300 . 0000 .772885 1.293853

500 . 000 400 .0000 .830456 1 .204158

500 .000 500 .0000 .876029 1.141515

500 .000 600 .0000 .914303 1.093730

500 .000 700 .0000 .947615 1 .055281

500 .000 800 . 000

.977308 1.023219

500 .000 900 .0000 1.004227 .995791

50 0 . 000 1000 .0000 1.028941 .971873

\
B 1°

T P D V

deg C MPa S / cm3 cm3 /g

750 .000 .0100 .000021 47219.519942

750 . 000 .1000 .000212 4720 .958893

750.000 1 .0000 .002123 471.103828

750.000 10 . 0000 .021679 46.127641

750.000 20 .0000 . 044390 22.527645

750.000 30 .0000 .068163 14.670731

750 .000 40 .0000 .092985 10.754408

750 .000 50 .0000 .118770 8.419521

750 . 000 60 .0000 .145340 6.880435

750.000 70 .0000 . 172425 5.799609

750 . 000 80.0000 . 199697 5.007588

750.000 90 .0000 .226804 4.409089

750 .000 100 .0000 .253421 3.946003

750.000 200 .0000 .462110 2.163988

750 .000 300 .0000 .584047 1.712190

750 . 000 400 .0000 .665974 1.501561

750.000 500 .0000 .728216 1.373220

750.000 600 .0000 .779010 1.283680

750 . 000 700 .0000 .822326 1.216062

750 .000 800 .0000 .860356 1.162309

750 .000 900 . 0000 .894435 1.118025

750.000 1000 .0000 . .925435 1.080573

T P D V

deg C MPa g

' cm 3

- cm 3/g

1000.000 .0100 .000017 58758.280969

1000.000 .1000 .000170 5875.475023

1000 . 000 1 .0000 .001703 587.198410

1000.000 10 . 0000 .017122 58.404407

1000 . 000 20 . 0000 .034420 29.052786

1000 .000 30 . 0000 .051866 19 .280576

1000 . 000 40 . 0000 .069439 14.401171

1000. 000 50 .0000 .087117 11 .478860

1000.000 60 .0000 .104867 9.535858

1000.000 70.0000 . 122648 6.153425

1000.000 £ 0 . 0000 . 140405 7.122229

1000.000 90 .0000 . 156081 6 .325872

1000.000 100 .0000 . 175613 5.694349

1000.000 200 .0000 .333285 3.000438

1000.000 300 . 0000 .451565 2.214522

1000.000 400 .0000 .539679 1.652267

1000.000 500 . 0000 .609226 1.641425

1000 .000 600 .0000 .666407 1.500564

1000 . 000 700 . 0000 .715336 1 .397941

1000.000 800 . 0000 .756357 1.318641

1000.000 900 .0000 .796939 1.254601

1000 .000 1000 .0000 .832065 1 .201830
B II

Sl <C> pjpjntDGiPj'PPicn.-iair^ Cii to to OJ tv OJ -T 0- OJ 0- S3 *-> (SJ HtOO rH 0- (SJ 0- TH to 'O OJ Ol to CL) OJ *« o- to rH

TS o (0 10 OJ U) ai o- <-t to rO CbJ PJ »i> to CT> fH O) CO H rH CO 'f in pj oj ts <-> S3 iT) 0- pjpicoCTicnncno-in (O 0- CO C\J rH O O* O* ‘O PJ c OJ O) 0- eg *h »-• oj to PJ to • 15 Si .O tO «-» CT) tO PJ CO tO 10 Oj t* rH c a> Th Ol Pj OJ Sl vf rTi CO C- ui UJ ea S3 rH O) O) OJ o- Si Th (g iD CD

CO vO (/j in \ pj to rH r-« pj <5w a’ to to 51 O- to CT) to Si Th oj in an Ol o 0- ai -»« N rH PI to (Jj tO Oj Uj to oj G, in oj OJ o- CL) o- dj S3 in X) Oi 0J Cii pH E & Tp gj e- O) .-1 tO r tO O- CT. rH to rH m in to rH S OJ rH E! OJ OJ •-> in Tjt O- OJ 0- 0- a > «-• 0- CNJ to 0- c* oj ai a) «f OJ K) o

oj r? r}< n if) iT) tTJ tn U J l£j lO CD {Q CD vL» to Sj •—i OJ to Tjt to X s <1) ID to to in ^ oj to 0* CO 01 S3 rH rH Th to a oi ® — OJ PJ to

P* fH f-H rH «-H rH H rH T-H H fH H iH rH rH OJ OJ Cj OJ OJ OJ OJ to 5* rH rH rH rH rH rH rH rH CVJ OJ 0J OJ 0J

rH Th •n r-t 0- oi PJ to to cj —i in m to oi Si <* 01 o- OJ CO CO S) to <* CD tO CO rH CJ <s »-• ui oj tn ® an in SJ tO CD 0- to CJ rP tO

to CaJ rH Sj -u Th PJ to Tji CO C\i f-H vh Sj O) n oj cn rH o- rH Sj C\j eg rH m o rH rH CVJ SJ PJ 0* 0- CO Si o- a> cc to a Pi C\4 tn oi C

S” H}4 <S to to to «-• Oj Th <s C^“ ivl Sl CO 0- IT) Ol tP 05 Ol to Si to in pj 01 0- 01 JMOO) c- in •-) cj o- uu to co cu o- an in o- Ol Pi Pi 3 1 to OJ 01 Ol 01 O to in PJ OJ r-t Sl CD 0- to to cn to CD gi Si h fr- SJ S i CD fs- Ui 01 0- to o- <o cn oi >n oi co cj rH CJ cj to rn O' rH CJ in

"h Th s> S3 Sj s- s* Sj Si s s s Ol Ol O) cu 0- 0- to •-o to lO Th \h 0J CJ to rji OJ C*J 0- to S. CD CO pj c- Pi Si CD vo in tn

1 ir -j* Cg 0/ OJ OJ eg OJ CVJ OJ CJ OJ H rt H rH rH rH rn rH rH rH rH rH 1 01 Ol Ol 0- vf Si 0- to to to ui in in m pj Pi Pi CJ CJ CJ OJ CJ

Cg eg C\J CJ - J CJ OJ CJ rH rH rH rH n rH rH r^ rn rH rH rH rH rn rH rH rH

ac 'a

Kj H O’ PJ Ol to CJ CD tO PJ l"- CD CD n rH Ol CD rH 05 an Pj Pi m Pi rH Si CJ CO Ol PJ Pi CJ Ol CO rH tO 0- rH in oi JV 0) o» in <D m O) OJ S. 0- og S) rp tn pj > CD rH tO •<p tn Ol pj ••P eg in (Si CL) Si T*1 TJ4 Tf< to <71 CJ n Tt< GO rH (Si in Th Si in 0- 0- o- to PJ c- PJ

rH Pi rn Oi n SiO s 'g 4 aj eg to Ol If) U) e- rH n rP rH in rr in 0- O) to rH rH ^Si Oi rH rH CJ u CD O) CO Oi n

rH Oi a. s. O C/ ill Pj CJ s< aD 0- in i-j li J hP CJ Ol CJ n rH (O 32 CD to s. s> tu rH CJ tO Ol CD Ol Si to to 00 o (S) to Pi Si a

Ol O) Si s; rH rH C\2 to ^ in to to o- CD Ol 0- tn Pi S) Ol in CJ Oi n 1 CJ eg rn rH s Ol Ol rH cn c^- 0- P- c- rH to OJ O) in CJ Oi n

m n CJ eg 0J CJ CJ CJ CJ CJ CJ OJ OJ OJ 0J to rP iX) to to > CO O) Ol CJ OJ CJ s gj S) > o- c- to to to to to to 0- 0- CD CD Ol Si <Si rH

CJ CJ Pi Pi Pi to CJ eg rH rH rH rH rH rH rH rn rH rH rH rH rH rH CJ CJ CJ

CJ rH SJ ‘A* to rH a> 0- CD rH Tp ts to s. m U) S5 OJ to 01

f-H n o ig tn s. tP O) rp S) m .—

i

to cn o to o Si rH m Pi CTl CJ c:» ,PJ KJ 0) Sj in (S) cr CJ o- pi or. tC CJ 0J pj "p rH pi CJ Si 0 e- s. s, s. O) O) m CD O) o- o to to in rH c^ PJ O) m rH to rH to

rH Sj > 0- e- to to <o to to to tc to to to n lO in 'p tP -P pj pi CJ co an

i

i

to Pi to to in cn CJ rH CD G) Ol to CJ Pi CJ Ol in rH o- s

to rH eg rH in to CJ cn eg O) PJ <o Ol Pi rH in m Ol o Ol <c Pi Pi

Ol iO CD g* Si 0- pi pi CJ tr. e~ si rH s> Ol If) CJ O) CJ ss OJ

CJ to CD cn CJ CJ pi eg to rn Pi s. CD O) e- S) Uj yf* s to

in •p CJ Si eg p O) ao o- cv to to gj g) CJ CJ rH rH Ol Ol 0) CD

‘ 5 Oi o (O m •p pi Pi CJ pj Pi Pi pj pi CJ Pi pi pj PJ eg CJ OJ CJ

i o

to CJ 0- 0J IJ cd cv in pj CO 0- rH rH ao Pj O Pi ao tn rH S) to o

S3 in m O) pi —• Pi Pi (S 03 S3 Si 0- to IS CD IS to OlHIOfl OJ

rH *rf •-O *^1 •p in m co t-i n CD rH CJ CJ cn l0 rH CU P CJ p in si <S 01 Si CJ rH CD CD in « (S o Ol m 01 CJ eg P to is O) Pi P CJ

Q. S) CTj rH rn Ol 01 CD <D S CJ m s m rH 0- to 0- o- CJ (O rH CD P a. tx CD P CD to cn «s c- m CJ to o- in Pi to Ol <s tn to «-* CO OJ 0J

CJ "v rH rH CD a» o »n PJ «- rtj CD to iO pj CJ SI s Pi CD to CD CD rH eg o \ Th m cj C'- ^ to 0- CJ sj eg eg o- in o- CJ Pi 0- Pi Ol m

O) Oi rH m rH rH rH rH rH s s si S) <s s a. CD > !> ao CJ Pi —3 s si *-< Ol PJ c- to Ol m eg Si CD 0- r-i 01 to o to to in in m

, rH rH g* n* ^ p -p -p -p -p -p tC p p fO Pi Pi Pi Pi PJ p p p 1 CJ CJ CJ pa eg sj rH

CO to m in in in p p p Pi Pi Pi pi Pi pi pi Pi

i

Pi o pi rH H* pi to 0> SJ SI Pi in CJ rH OJ to to p) to CJ SJ Pi Pi

in co rH U) eO eg p to CD (S ts CD p Oi 03 rH Pi O) m O) PJ ^ CO oto CJ 0- eo eg a; if) Pi Pi Pi Pj m o S. in CJ CJ to rH P Ol p ts

pi p cg eg eg S3 o- tn pj rH Ol 0- m lO eg CD 01 OJ to 0J CD CO SJ Pi p p eg s» Si St 0> Ol Ol o. aj cd CD CD CD to m ai p p PJ Pi p p

h th ^ pj fj pj (O to fj pj f j io i*j to pj n n pj pj to to

CD in cn rH O) in CD CO to Pi to CJ cn m CD to Si pi p- 3 Ol 3

p pj to cu 0

- rH O) Ol CJ O) CD Si 0- pi to in Ol 3 Ol to in > to Oi rH ts pi rH an gj OJ CTl •s m pj to 0- CO CD to tP Ol T|» OJ o ui ao CJ pj o- > !0 eg <D to H* p PJ pj P c^ S) PJ to O) rH to

in m to rH PJ rH Si O) a) UJ CO CD CD CD ao CD Ol 0i 0) ai s. Si rH rH rH CJ Pi to pi CJ CJ CJ CJ CJ CJ OJ CJ CJ eg CJ OJ CJ IO P. PJ

eg CO 1 O) to OJ CD tO i-i ai ai in pj pj pi to p Pi p rH ui to OJ co to

rh eg » UJ rn C'- P; rH p O) PJ 3 pi pi r-» HJ UJ cd in ui HlMlJ c CJ Si 1 c- «.g cj p in P CJ rH tO in to CC 3 tii CJ Th pj eg s rH tt

CD CD 1 e» CD pi OJ in a- Oj CD r^- to eg 3 Ol S3 pj s S3 in rH O) Ui

cu P 1 co Oj Ol p O) TT* O) t

Ol Ol Th Ol 0- U. Uj 3 pi tfj Th to in

TD H r-» 1 tg eg CJ pj PJ Th in eo to c~ CJ L^ OJ 3 pj C» Pi CJ

I

1 CJ CJ CJ eg CJ CJ eg CJ 0J CJ CJ 0J eg Pi PJ g* Th m m Ui to orH cn

t

1 rH Ol in Ol rH in UJ Ol 3 S 3 CD Ol CD C- rH o pi to CD Pi 3

E-» O PT PJ 1 PJ in s rH rn pi eg s CD PJ Th O- s e- 0- o> CD eg rH CO rH GO

T3 "hi Sj Si 1 to s in cn pi rH eg m CJ -H eg gj CJ in m Cl Ui Th rH 3 tO

Ns a. S3 3 1 C* J pi to cn in eg 01 to pi t-H Ol > m Ol o PJ Th pi PJ 0- to

CL. 21 Cg s 1 J P J PJ m ao rH P) to tn eg e- 3 CO CD UJ Ui PJ 0- P CD

rs Si Si 1 s s s s s rH r rH r eg CJ CJ Pi n CD 3 pi PJ 3 CD ^

1

1 rH rH rH rH r rH rH rH rH rH rH rH rH rH rH eg OJ CJ PJ p p in

0- PJ

1

1 rH (S» eg in to tO 3 Pi 3 in oi co pj e- ^ s 3 to Ui to to 3 Pi to CD 1 CTi PJ eg S3 eg in 3 to in to 3 CD 3 Oi o rH P Ui to UJ 3 £ 3 (S 1 O) s HP pi in <o 0- tu in pj P o Th P 0- a CD PJ tO O Pi «-t o o 3 Si 1 o CD CD CJ to 3 P CD CJ to 3 PT e- s CD Pi ui to n pi 3 in \ 3 3 1 OJ a: CD Ol Ol 3 3 3 f—

»

rn eg cj oj to 'JJ rH pj m o ai rn eg

S. Si s 1 Oi 03 03 Oi Ol S S s S s s s s 3 rH rH rH rH eg eg

1 rH rH rH rH rH rH rH rH rH rH rH rn rH rH rH rn rH

toncoajc^ojtoairpioirQifipitvcn'i'rpiPi^a)

»t<tooic~t>cac-ajvjOjioaic\ipJ'»t<cMio«-*»-t'«t<(OrjCTi

o s w 01 <j> a tp c- in <-j h to rj pj uj ® s in pj w ai 'p «<

TD \ CTi CO O) PJ S3 to eg rH S3 CJ Ol OJ rH S3 rH Th PJ CD CD X) 3) in

(X 0) Ol CL 01 in Pi to C“ to Th PJ r

H

Cl n tt o. 3 Cl to rn n

T3 CJ OJ CJ rH rH CJ PJ Th m to £> UJ -*h 3 in CTl Th o r^ n UU rH CJ CJ CJ PJ Pi Th Th Th

in m CJ rH in S3 to Th to cn to <o 3 ID O) CD PJ C* in rH 3 CJ 3

E-t rH m rH e- Ol OJ m CT) to Ol O) Si U) CJ rH m O) CJ rH Th Ol eg

T3 s rH g^ Th PJ pi eg rH in Pi rH 0- p j S3 UJ 3 Th sh n Th CD eg n N cu s S3 rH tn pj CD Th rH n Si pi c^ to S3 ai o in OJ rH 0- 3 Th Ol

ou z; s Si S3 eg rH <o rH m Th eg CD p > cu PJ eg S3 pj eg S3 u) 3 pj n

73 S3 3 3 rH CJ X) 0- u3 O) O) 3 3 rH Th gj 0- UJ O) cn 3 3 S3 rn rH rH rH rH 1 1 1 1 2 2 CJ

pi in in pj 3 oi Th tO PJ P 0“ 3 to in 3 rH Th (O o- 3 0-3 S3

PJ pj Pi CTl tO ig Si Ui Ui CJ to tO Ol I> Pi CJ Th *«h (Ji uJ CJ PJ 0- rH E S3 Pi PJ 0- 3 eg Ul PJ P Ui 3) Ui Ui Ol tO OJ Pi to p m n cn a CJ S3 S3 Pi 3 S3 Si co oi p m Th rn in OJ to pj > pi P OJ > rH

V. 3 S3 3 P Pi CJ Ui 3 P to X 3 rn eg p o- rH lf> to p pj 'O CD - S3 Si ^ CS p Pi ui to to to -.0 O- c- CD CD CT) CTi Ol S3 3 3 3

r-» rH rH rH

Si •h i tt S3 Sf Si S3 Si S3 S3 S3 S3 S3 Si (S3 ea S3 S3 ea ea «i Sl ts ts. 3 pj i p J Si CS3 S3 Si 3 3 Ci S3 S3 Si Si S3 S3 Sl S. S3 (S3 S3 G. S3 S3 3 m s> CJ i eg S3 3 3 S3 Ci tS CJ S3 S3 S3 S Si S

;

SI <s S) Si Si Sr S3 S3

CU rH rH i « H S. S3 Si <s Si S3 S3 Si Si Sl Si Si S3 S3 Si S3 Si Sl Sl S3 S3

r. S3 s. i s r—

1

S3 S. S3 S3 Si Si Si Si S3 S3 S3 si Si ss Si Si Cv S3

i t-h Sj g s g S3 Si (S s s S3 si Si S3 S3 S3 s. S) Si

i rn i j pj p a

)

eo O cu fT> S3 S. S3 S. Si Sl S3 Sl (S s.

i

i ig tn ^ lO to X Cl S.rH

S3 S3

i

i 3 Si (53 S3 S3 S- S3 S3 S3 ea (Si Si Si S3 S) ts S3 Si Sl s 5J ®

o S3 3 Si Of S3 Si Si 3 3 SI Si Si <s. SI t. S3 S3 Si S3 CJ S) S3 SJ S. 3 S3 S3 S. S3 S3 S3 Si S3 ca S3 (S3 Si G3 S3 ca Si Si S3 S3 (Si (Sl S3 S3

£-• JD

ai 3 CJ S3 S3 3 3 S3 S3 S3 Cd <S S3 Si S3 Sl S3 S3 S3 ea s cj s s sj

•o in m Ul UJ n n n n n n in uj uj n m Ui in n Ul n ui n n n

.

S3 Sl <SJ S3 S3 S3 <s S3 eg Si s <s S3 S CjJ s S3 s Sl S) Si O

S3 Sl Sl ea S) CJ SJ Si S3 Si S3 S3 Si S. S3 S3 S s Si S S S3 S. CU 75 S3 Sl Si ea s; rn Si S3 S3 (S3 S3 (S ea S, 5; s s s S3 s s S3 (Sl uU rH Si Sl S3 si eg S3 Sl vjl Sl Sl S3 ^a Sl S3 s s s Sl S) c3 S3 Sl 2. Sl rH Si >3 Si Pj eg •g Si Sl G. S3 S3 Si S3 c.a s s Sl >3 s. Sl

S3 (Si eg S3 ea S3 Sl S S3 Si s. Si P; s s s CS Si s s rH T.g ^g p Th uj <0 O vXi * Ti Sl Sr s s s s Cs Sl S3 s

rH 1. s« pj *h U 4 e;> O u. Vji s.

ea Sl S3 S3 S.1 s. S3 Sl Sl <S cj Si S3 S3 S3 S3 Si Si ca SJ S3 S. S3 o Sl S3 S3 S3 51 tS) S3 Si S3 S3 Si Si S3 S3 S3 Si Cv s CJ 3 S3 S3

S3 Si S3 Sl Si G. S3 S3 si eg G> Si g s. s, Si (S3 Sl (Si pg ca s) ca

6-1 OJ in n n uj ui ui m in uj in u j in u‘> in uj UJ in uj in uj in in a)

73 o 0- 0- 0- o- e- 0- 0- C^ o- o- o- c- 0- 0- o- O- 0-

pj pj pj pj pj pi pj PJ Pi PJ PJ PJ Pi cn pj pj PJ PJ PJ Pi PJ Pi PJ
I

I
U.S. DEPT. OF COMM. 1. PUBLICATION OR 2. Performing Organ. Report No. 3. Publ i cation Date BIBLIOGRAPHIC DATA

SHEET (See instructions)

REPORT NO.

NBSIR 81-2253 1981

4. TITLE AND SUBTITLE

A Theimodynamic Surface for Water: The Formulation and Computer Programs

5. AUTHOR(S)

Lester Haar, JohnS. Gallagher and George S. Kell

6. PERFORMING ORGANIZATION (If joint or other than NBS, see instructions) 7. Contract/Grant No. national bureau of standards

DEPARTMENT OF COMMERCE 8. Type of Report & Period Covered

WASHINGTON, D.C. 20234

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City. State. ZIP)

10. SUPPLEMENTARY NOTES

XDocument describes a computer program; SF-185, FIPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual summary of most si gnificant in formatioif. If document includes a si gn i fi cant

bibliography or literature survey, mention it here)

A FORTRAN 77 program is given with which thermodynamic properties for liquid

and gaseous states for water can be calculated for the temperature and pressure

range, 250 <_ TCK) <_ 4000 and 0 <_ P < Pmax , where pmaxis the lesser of the

pressure of melting ice or 4 G Pa (T0,000 ? bar) . The program yields values for

enthalpy, entropy, heat capacity, speed of sound, Jcule-Thomson Coef. , second

virial coef., isothermal compressibility coef., coef. of thermal expansion, and

PpT relations that are everywhere in accord with accurate data for water.

Included is a brief outline of the derivation of the equations

.

12. KEY WORDS (Six to twelve entries; alphabet! cal order; capitalize only proper names; and separate key words by semicolon s) Enthalpy, entropy; equation of state; heat capacity; Hehnholtz function; internal

energy; Joule -Thomson coef., second virial coef., speed of sound; the critical

region: thermodynamic surface: water.

13. AVAILABILITY

[ '/

\

Unlimited

| | For Official Distribution. Do Not Release to NTIS

I | Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402.

Order From National Technical Information Service (NTIS), Springfield, VA. 22161

14. NO. OF

PRINTED PAGES

35

15. Price

USCOMM-OC 6043-P80


»

.