New forms of state equations for helium

Accurate helium properties must be measured to develop refrigeration systems for superconducting motors and generators. Normally when measuring these properties a double iteration is necessary on the equation o f state to find density and temperature. In this paper new state equations are developed which eliminate iterative procedures. These equations are subsets of three overlapping sets o f reference equations. They are thermodynamically inconsistent to the reference equations by a small degree o f error, but have been found useful for many calculations.

New forms of state equations for helium V. Arp

The development of large-scale helium cooling systems and refrigerators is being stimulated by rapid advances in superconducting technology. While many superconductor designs utilize simple pool boiling in a one-atmosphere helium bath, attention is turning towards more sophisticated cooling systems involving forced helium flow at higher pressures, for example, as required by superconducting power lines. 1,2 Outside the US a number of force-cooled magnets have been built for high energy physics applications. 3-9 Superconducting motors and generators will require forced helium flow to maintain operation. In most of these applications, detailed refrigeration system design requires the use of accurate helium properties, especially, (1) the enthalpy, the change in which is equal to the heat or work input, (2) the entropy, the change in which measures system irreversibilities and refrigeration requirements, (3) the pressure, which is controlled by compressors, valves, and frictional losses, and (4) the temperature, which must be controlled to conform with superconductor requirements. A difficulty which arises in practical calculations derives from the use of density and temperature as the independent parameters in the helium equation of state and transport property functions. This is necessary since most fluid theories are based on these parameters. However, as a consequence, given any pair of the three parameters P, H, and S (calculated from physical requirements), a double iteration on the equation of state must be performed to find density and temperature, and other properties such as compressibility or viscosity. This can be costly with regard to computer time, and requires a carefully chosen initial estimate to be practical. The recent publication by Hands lO reports the development of such iterative routines. In this paper, we report a different approach: the development of new state equations for helium to eliminate the need for iterative procedures in flow and refrigeration calculations.

Selection of equations A study of a variety of flow and refrigeration calculations indicates that the following four equations will allow a variety of problems to be solved without iteration. Pressure and enthalpy are chosen as the basic parameters. To obtain The author is with the Cryogenics Division, Institute for Basic Standards, National Bureau of Standards, Boulder, Colorado 80302, USA. Received 6 August 1974.

CRYOGENICS. NOVEMBER 1974

other properties from the conventional state equations, we need the temperature T(P, H) and the density p (P, H). If the entropy and pressure are known, the enthalpy is calculated using the third equation, H (P, 5). Finally, it is convenient to have the entropy S (P, H), though it can also be calculated more indirectly using the first two equations and the conventional state equations in sequence. In principle, the density and the temperature could also be calculated from the third equation, viz, - V = -OH (P, 5) a n d T = - -OH t' (P,S) p ~P Is ¢3S 1

However, errors as much as an order of magnitude may be introduced by the differentiation, reducing the accuracy of this procedure below an acceptable level.

Development guidelines The equations are based upon the recent helium properties data 11 and state equations 12 developed by McCarty in a thorough analysis of all available experimental data in the range 2 to 1 500 K with pressures to 1 000 atm. To obtain high accuracy over this very wide range, the reference equation of state uses three separate equations covering three overlapping regions which are defined in the density temperature plane. Careful splicing techniques were utilized to obtain smooth transitions in the regions of overlap. For this work we have chosen to express each slate equation by a single mathematical function over a more restricted range, 2.5 to 400 K with pressures to 100 atm, although the resulting equations retain the accuracy of a simple virial approximation considerably above 400 K. For most studies of helium cooling systems, there is little impetus to use pressure above 10 or 20 atm, so that increased weight was given to the fitting at low pressures. A total of 944 data points within this range were generated from the equations and entered into the least squares analysis. The range of data is shown graphically in Fig. 1. All calculations were performed using units o f P (atm), V ( c m 3 gl), H ( J g-l), S ( J gl K-l), T(K), and R = 2.0772258 J grl K-1. Further, the evaluations of enthalpy and entropy use datum values consistent with those in McCarty's 1972 tabulation, 11 which are greater than

593

those calculated from his 1973 descriptive text 12 by 58.980 J mole-1 and -0.443 J mole-1 K -1 respectively. Also, our temperatures like those in the reference equations are greater than those calculated from the common Ts8 scale by an amount AT = 0.001 + 0.002 T for any stated pressure along the liquid-vapour equilibrium line.

( Not to scale ) I00 E u

7,

T h e transposed critical line

The transposed critical line also plays an important role in normalizing the equation structure. It is shown in Fig. 1 as the dashed extension of the liquid-vapour equilibrium line, and is the loci of maxima in the specific heat, as seen in Fig.2. More exactly, we define the transposed critical line as the locus of points where

aG I aT p

0.1 2.5 3.5

I0

400

r, K

(/', 7) = o

Fig.1

Different transposed critical lines would result if other equally plausible definitions were used, such as acp/aP = 0 or if either the expansion coefficient a or the compressibility K were substituted for Cp in these above. The differences between transposed critical lines defined in these different ways would be very small near the critical point, but relatively larger far from the critical point where properties maxima become small and spread over wide regions of the P-T plane.

PT plot showing phases& phaseequilibria

SC 2.0

40 30

•

3.0

20 5.0 P= 1.0 atm

Structure of the equations

'v

-

To~

The development of the equations starts with a simple virial expansion, fitted exactly at a selected high temperature and at densities far less than the critical density. The virial expression is then extended to the entire fluid properties range, and a relatively complicated deviation function is added to correct for the'inadequacy of the virial expression. The virial expansions for temperature (TO, entropy (S/), enthalpy (H/), and specific volume (V/), derived from the reference equations at 400 K, are 2H

TI(P, H) = - -

5R

-~ ~

,I,A

IO

I0.0

8 6 5 4

P =20.0 atm

3 2

- 2.8239 + 0.64469 P - 5.615 x 10-6 p2

5

(1)

SI(P, 7) = -~R loge T - R

loge P + 1.79056 + 0.0001077 P

t

[

4 Fig.2

Pitt of specific

I

I

I

6

I

I

8 T, K

I

I0

t

1

12

I

I

14

heat against temperature

(2)

and

V/(P, 7) = 2.8476 + T

( Ra ~ - - 5.63 x lO8 RaP )

(3)

where R a = R/0.101325 = 20.500625. Equations 1 and 2 can be combined to

HI (P' S) = 3"6785757P°'4 exp ( S-00001077P )2.5R + 14.6649 + 0.33479 P - 000002913/~

594

(4)

Here and throughout the work we retain more significant figures than are justified by the basic experimental data, to maintain internal mathematical consistency in the computed results. Equation 3 is a hybrid, in that the second virial coefficient multiplies the pressure divided by the temperature, and the third virial coefficient multiplies just the square of the pressure; the justification for this form is that it results in a better fit than does the two-term expansions of either pressure or pressure divided by temperature. The deviation between real properties and those predicted from the virial expressions become non-zero at enthalpy values of about three times the enthalpy at the transposed critical temperature for the given pressure. This empirical fact is not universal, in that it must change as the zero of enthalpy is changed, but forms a convenient basis for further

C R Y O G E N I C S . NOVEMBER 1974

structuring the equations which are developed here. This is illustrated in Fig.3, which is the deviation between the true temperature and that predicted by the virial expression, (1), for various pressures. The deviations are abrupt at the higher pressures, when plotted on a logarithmic enthalpy scale. This abruptness causes some difficulty in fitting the higher pressure data. Similar plots are found for the other virial equations.

Htc = 16.3713 + 2.21774 P - 0.002606/2

,¢ ,¢

(5)

m

"

Residual uncertainties in these last two equations deserve brief mention. The specific heat and the enthalpy are related to the equation of state through certain derivatives and integrals over the PVT surface. It happens that the region between two separate mathematical functions in the reference equation of state, and properties in this region depend upon a weighted average of the two functions. This is entirely adequate for calculating the usual properties of interest. However, we find that the locus of the transposed critical line, which depends upon third derivatives of the PVT surface, defines distinctly different curves, respectively, below 10 K and above 15 K, with an unrealistic transition in the overlap region. Further, we found that, at higher pressures within this region, the enthalpy calculated by the reference equations differs slightly from f CpdT along an isobar, though the integral over the whole range from 10 to 15 K is correct. The data points within this region were given reduced weight, and (5) and (6) are smooth compromises. The equation

T (P, H)

(H/Htc)

z = n/ncpi

p ;loge

(1 +P)

rr = exp ( -

P/4)

9

20

4060100200400

j g-I

1000

40(30

K-t

Fig.3 Deviation between true temperature and that predicted by equation (1) for various pressures

fori= 5 ~(1 +Ws loge (z/2)) Fs (Ws) =

1 +(z/2) s

for i = 6 and 7

Fi(wi) = rr X exp (-wi X2) and for i = 8 and 9

(logez)Wi p(Wi-2) Fi(wi) =

1 +z 3

Using the definitions

CRYOGENICS. NOVEMBER 1974

H Hte

(9)

t=l°ge I (I + TI(P'H)]20) 1 2 1

2

l )

h-

and

/=0

where, for i = 1 to 4

( l +(1/a)wi 1 + (z/a) wi

The function p (P, H)

4

E Z affFi(wi)pj i=1

Fi (wi) = 1 --

6810

(7)

the temperature as a function of pressure and enthalpy is given by

T (P, It) = TI (P, H) +

4

The constants aii and wi are given in Table 1. The standard deviation between (10) and McCarty's equation is 0.025 K. The equations were so weighted that the accuracy at the liquid-vapour equilibrium line and within a factor of two (in pressure and temperature) of the critical point is higher, generally on the order of + 0.015 K. Along the 100 arm isobar, the accuracy is rtoach less, generally on the order of + 0.09 K.

Using the additional definitions X = loge

2

(6)

This gives the location of the maxima in the curves of Fig.2.

\\

',,',,\ I \

?2

As the temperature drops below the transposed critical temperature, for a given pressure, Cp drops from a maximum to below the ideal gas value, (5]2) R. The enthalpy at the point where Cp = (5/2) R is given approximately by Hcpi = 8.6467 + 1 . 3 2 2 0 P - 0.001823 p2

I

;3

The enthalpy at the transposed critical line is given approximately by

.a=O. 6 an effective temperature can be calculated

595

T a b l e 1.

N u m e r i c a l c o n s t a n t s f o r t h e e q u a t i o n T (P, H )

The

S (P, H)

function

Using the d e f i n i t i o n s given b y ( 7 ) and ( 9 ) , d e f in e al0

=

0

a53 =

0.2116953534

all

=

0

a54 =

0

a12

=

a13

-

a14

=

1.6915507264 1.2909577607

7

a60 =

--11.157137720

a61 =

20.696139062

0.1706960231

a62 =

0

a20

=

a63 =

1.6015258719

=

25.020386928

a64 =

0

a22

= -

28.518685534

a70 =

a23

= =

a30

=

a31

=

a32

=

a33

=

a34

-

a71 =

1.4079442213

a72 =

0

0

a73 =

-

a74 = a80 =

0.1113530701

-

0.0746631692

a81 =

-

0.1033700894

a82 =

a40

=

0

a83 =

3.7959000872

a41

=

a84 =

0

25.813615600 -24.155499152"

a42

-

7.1881925681

a90 =

a43

=

0

a91 --

a44

=

0.1727460143

a92 =

a50

=

5.7871515264

a93 =

a51

-

1.2338412546

a94 --

a52

=

0

Wl

=

2.10

W5 =

W2

=

1.85

W6 =

12

W3

=

1.60

W7 =

3

W4

=

1.35

W8 =

2

W9 =

3

7.0898285089 -

8.1283258235 2.5560991610

-

0.2875898804

Then the entropy is found by using T* in the virial expression:

s(P, H) = s/(p, r*) The coefficients

+(1+h4)-1

~ /=-0

function

s = exp

H (P, S)

(2s)

and using the definition given by (7), we find

8iti +(1 +h2)-I Z

H (P, S) =HI(P, S) ~6+ip°'2i

+IE,

/=I

7

5 bii(s.1)(1-i/3)

p kl

5

Z

i=l

0~](1) '/2p0"2]

]=0

The coefficients

j=l

bq are given in Table 4.

The standard deviation between the calculated enthalpy and that given by the reference programme is 0.050 J gq. If we relate this to a temperature error by

and finally -

Vi (P, T*) 6T-

The coefficients 8, and dij are given in Table 2. The terms 8it' essentially describe the temperature-dependent second virial coefficient. The standard deviation between this equation and the reference equations is 0.28%. In the near critical region, within a factor of 2 of P and T at the critical point, deviations of 0.6% are common. Also, deviations of up to 1.5% are found in the compressed liquid at pressures below 1.5 atm.

596

Cij are given in Table 3.

The standard deviation between the calculated entropy and that given by the reference equations is 0.017 J g-1 Kq" If we relate the error at each point to a temperature error by 6T =(T]Cp)6S, the standard deviation in temperature is 0.014 K. The error tends to be higher than this by a factor of about two at temperatures above 100 K and pressures above 5 atmospheres, and largest in the compressed liquid below 1 atm, where it reaches a maximum of 6S ~ 0.2 J gl K-1 In the near-critical range and along the liquid-vapour equilibrium line, the error does exceed the above-stated standard deviation at only a few percent of the data points.

Defining

0.137115

5

4

O(e, H)

T* = TI ( P, 14")

0.0127568335

6

/=1

#,

and then

The

T*=T[(P,H)+P Z

-(3-:')

0

4.7645900164

10.564283801

C,ih

1=1

,=1

1.3737206422

5.8298990782

34.295641395 -

-

0.2162115516

-33.135574167

x x

-11.089866761

a21

a24

=H+

m

5

6H

the equivalent temperature error is approximately 0.010 K. Over most of the fitted range this error is rarely exceeded, the exceptions being below about 4 K where temperature equivalent errors up to about 0.03 K are common, and at pressures above about 70 atm where temperature equivalent errors up to about 0.05 K are found.

CRYOGENICS.

NOVEMBER

1974

Table 2. Numerical constants for the equation p (P, H) x 101

d23 =

0

4 . 2 8 9 0 6 7 6 0 6 1 x 102

d24 =

2 . 2 3 0 5 0 2 8 8 7 1 x 102

doo

=

do1

=

-8.8532361747

d02

=

do3

=

do4

=

do5

=

dlo

=

dll

=

d12

=

d13

=

-5.1430448958

x 102

d40 =

d14

=

-1.8239511033

x 101

d41 =

0

d15

=

3 . 8 6 8 1 0 3 9 9 9 8 x 101

d42 =

3 . 3 1 7 2 6 2 2 9 2 1 x 101

d2o

= -1.1451954585

d21

=

d22

=

-5.1538182506

x 102

d45 =

-8.5023300915

61

=

--1.6861363337

x 102

67

-4.1151137531

62

=

4 . 3 7 6 1 8 2 1 1 1 1 x 10-3

68 =

63

=

0

69 =

-6.7965868266

x 102

d25 =

4 . 3 2 6 7 3 5 1 5 4 6 x 102

d3o =

-9.0266448352

x 101

0

d32 =

1 . 6 8 2 1 3 2 6 8 4 2 x 102 -7.8689686748

d34 =

1 . 0 9 9 4 3 0 6 5 5 8 x 103

d35 =

x 102

d43 =

4 . 9 4 3 6 9 4 2 0 2 7 x 102

d44 =

=

64

=

- - 6 . 1 5 3 6 0 1 0 1 0 2 x 10 -3

510 =

65

=

- - 3 . 6 8 9 4 2 9 4 0 1 2 x 10-3

611 =

66

= --6,4620070038 x 104

If the input data to any of these four equations corresponds to a point within the l i q u i d - vapour equilibrium region, the calculated property will be meaningless, but this fact will not be obvious by virtue of any unusual numerical magnitude of the result. To test for this error, it is necessary to find whether P and H define a two-phase state point. The enthalpy can either be input data, or be calculated from an input entropy, through H (P, S). Within the two-phase region, the enthalpy is always less than 30.2 J g-l, and the pressure is less than 2.2449 atmospheres. If the point in question falls within this range, a more accurate test must be used. The enthalpy of the saturated liquid as a function of pressure is given by

'7

i

[2 ( J)l

1

2.2449

i=l

and the enthalpy of the saturated vapour, H v, is given by an equation of the same form, but with coefficients designated)~v. The coefficients f//and fly are given in Table 5. If the enthalpy lies between the two values defined by these two equations, for a given pressure, the point in question is in the two-phase region.

CRYOGENICS.

d33 =

x 102

The two-phase region

Hi = 21.36 exp

d31 =

NOVEMBER

1974

-6.8986815155

x 101

3 . 0 7 2 7 2 5 2 6 9 1 x 101 -1.0376918669

x 102

0 1 . 9 9 0 2 7 2 8 4 4 2 x 102 -1.7695524072

x 102

4 . 1 9 6 2 6 0 3 9 8 1 x 101 -2.0829009386

-6.1261451275

x 101

4 . 0 1 7 9 4 0 9 1 0 8 x 101

1 . 5 7 3 3 2 5 4 8 8 6 x 101 -2.1206083607

x 101

1 . 0 3 5 4 7 0 5 9 1 7 x 101 -1.6690511792

The standard deviations between these equations and McCarty's data, stated as a temperature error by the relation 6T= 6H/Cp are about 0.014 K for both equations.

Discussion The reference equations ]2 are a set (that is, three overlapping sets, as described earlier) of equations for pressure, enthalpy, entropy, Cp, Cv, sound velocity, and various derivatives, all as a function of density and temperature. These equations are internally consistent although the accuracy of the fits to the various experimental data is in the range of a few tenths of a percent to a few percent, as explained in reference 12. The equations developed here are a subset c,f the reference equations, set apart by being thermodynamically inconsistent with the reference equations by tile amounts stated as 'error' in this paper. These errors are generally small compared to the fitting accuracy to available experimental data, though they may be important where an unusually high degree of internal consistency is required of the computation. In such cases, the equations could be useful for initial estimates in iterative routines with the reference equations, following Hands. 1° However, we have fi)und the equation to be very usefid for a wide variety of helium refrigeration and flow calculations. The author is grateful to R. D. McCarly for numerous discussions during the course of this work. The work was supported by the Advanced Research Projects Agency.

597

T a b l e 3. N u m e r i c a l constants f o r the e q u a t i o n S (P, H) Cll

=

96.557792404

c44

=

90.187334672

bll

=

b12 =

c12

-

85.342124216

c45

=

c13

=

93.239184634

c51

=

405.81584516

b13 =

c14

=

--32.642437375

c52

=

123.84837652

b14 =

c15

=

c53

=

2.8950104636

C21

=

--471.18800861

c54

=

c22

=

337.90821505

c55

-

c23

=

--195.19240224

c61

=

c24

=

c25

-

c31

=

70.146358374 6.2193650992 884.01164347

c62

-

c63

=

c64

=

8.3939686508

Table 4. N u m e r i c a l constants f o r t h e e q u a t i o n H (P, S)

-321.35063183 89.302684169 8.2508383745 -102.61607953 72.105848121

b21 = b22 = b23 = b24 =

1 . 1 8 2 8 3 4 4 0 3 7 x 103

x 101 b45 = - 9 . 3 2 1 8 9 3 9 1 7 4

x 101

2 . 2 2 9 7 2 6 1 6 7 4 x 101 b51 = - 4 . 4 9 9 8 4 0 6 7 2 4

x 102

-6.7492782747

b52 = - 5 . 4 9 4 5 8 9 2 5 4 0

4 . 0 1 2 1 5 7 3 6 9 1 x 1G 1 b53 = -5.2643461931

4 . 9 5 5 0 3 4 1 8 0 4 x 103

x 10 -2 b54 = - 1 . 3 8 1 9 2 1 2 7 6 3

3 . 6 2 1 6 9 8 8 3 7 6 x 102 b55 = -3.2737113461

x 102 561 =

2 . 9 7 2 6 9 5 6 8 0 9 x 102

9 . 9 1 8 6 5 7 4 0 1 1 x 101 562 =

3 . 2 6 7 8 0 8 3 3 3 8 x 103

-7.0819815035

b63 = - 2 . 8 9 3 3 1 2 4 9 9 2

34.006661384

531 =

-5.2181285811

x 101 564 =

-1.8942905449

x 103 b65 = - 6 . 2 8 7 1 9 2 5 1 3 8

x 101

1 . 7 1 0 7 6 0 2 3 6 5 x 103 b71 = - 6 . 6 2 7 7 4 7 1 1 5 9

x 101

= --359.83067977

c65

=

3.140489957

c33

=

C-/1

=

10.469390204

b33 =

c34

-

6.9236428132

c72

=

11.920614672

534 =

c35

=

0.4327647102

c73

=

16.687143281

535 =

3 . 8 9 2 1 6 7 6 3 7 1 x 101 b73 =

2 . 5 4 2 8 6 3 0 0 5 0 x 102 b74 = - 1 . 7 4 4 7 7 0 8 4 8 7 4 . 5 2 5 5 6 8 2 3 2 6 x 103 b75 =

c41

=

c42

=

c43

=

--821.10597373 51.518553781

c74

=

4.7025518338

b41 =

c75

-

0.4348173542

b42 =

341.82700671

R e f e r e n c e s

1 2 3 4 5 6 7 8 9 10 11 12

598

Meyerhoff, R. W. Cryogenics 11 (1971) 91 Forsythe, E. B. Brookhaven National Laboratory Report No 50325 (1972) Anashkin, O. P. et al Prib i Tekh Eksper No 3 (1968) 218 Morpurgo, M. Particle Accelerators 1 (1970) 255 Lesmond, C., Lottin, J. C., Shimamoto, S. Proceedings Int Conf on Magnet Technology, Hamburg (1970) 925 Bailey, R. L., Colyer, B., Homer, G. J. Rutherford Laboratory Report RHEL/R285 (1972) Berthet, M. Proc ICEC4 Eindhoven IPC, (1972) 209 Stevenson, R. Cryogenics 13 (1973) 524 Vecsey, G. Proc ICEC4 Eindhoven (IPC, 1972) 230 Hands, B. A. Cryogenics 13 (1973) 423 McCarty, R. D. NBS Technical Note 631 (1972) McCarty, R. D. JPhys Chem RefData 2 No 4 (1973)

b43 =

-5.0693398133

x 102 b72 = - 7 . 5 4 7 8 7 6 0 1 7 4

-4.1189023813

x 103

7.8724009925x

c32

-

x 103

1 . 1 0 0 1 8 4 9 7 8 7 x 102

532 =

-- 2 4 . 9 6 1 1 4 1 5 5 7

x 103

b25 =

120.99391576 -

b15 =

6 . 2 0 9 5 6 5 2 6 5 3 x 10-1 b44 = -2.4559275204

102

x 102

6 . 5 4 2 7 8 1 7 2 2 1 x 102 x 102

1 . 3 9 3 9 3 6 1 5 3 6 x 101

x 103

Table 5. N u m e r i c a l constants f o r t h e e q u a t i o n s H (P) at the two-phase boundaries f1£ = f2£ = f3£ = f4£ = f5£ = f6£ = f7£ =

--5.9177497274

flv =

4.5307631243

4 . 7 2 8 0 1 2 4 3 9 3 x 101 f2v = - - 3 . 3 4 9 2 3 6 1 0 7 8 x 101 - - 2 . 2 2 2 2 8 4 6 3 0 9 x 102 f3v =

1 . 3 5 8 2 8 6 7 3 5 3 x 102

5 . 5 9 7 0 5 0 2 6 6 6 x 102 f4v = - - 3 . 0 3 0 3 2 8 9 0 5 9 x 102 - - 7 . 7 3 0 4 2 6 4 0 5 2 x 102 f5v =

3 . 7 2 0 3 0 7 4 2 9 7 x 102

5 . 5 1 0 9 9 4 8 3 0 3 x 102 f6v = - - 2 . 3 5 4 9 2 1 0 5 0 3 x 102 - - 1 . 5 8 9 2 2 6 1 0 1 3 x 102 f7v =

CRYOGENICS.

5 . 9 7 9 1 3 0 5 6 4 8 x 101

NOVEMBER

1974