A111DS b23773

NATL INST OF STANDARDS & TECH R.I.C.

A1 11 02623773

Morrison, Graham/Application of a hard s QC100 U5753 NO. 1226 1987 V19 C.1 NBS-P

lMfnHiflK3!l!HflVlfflll!OTi*i

NBS

PUBLICATIONS

NBS Technical Note 1226

Application of a Hard Sphere

Equation of State to Refrigerants

and Refrigerant Mixtures

Morrison and McLinden

c NBS NBS NBS NBS NBS NBS NBS NBS NB , . . ?S NBS NBS NBS NBS NBS NBS NBS NBS NBS N

NBS NB BS A BS NBS 'S NBS NBS NBS NBS

IS NBS NBS NBS NBS NBS NBS NBS NBS NBS N

BS \B BS NBS NBS NBS NBS NBS NBS NBS

IS NBS NBS NBS NBS NBS NBS NBS NBS NBS N

NBS NBS NBS NBS NBS NBS NBS NBS NBS NBS

IS NBS NBS NBS NBS NBS NBS NBS NBS NBS N,

NBS NBS NBS NBS NBS NBS NBS NBS NBS NBS

IS NBS NBS NBS NBS NBS NBS NBS NBS NBS N

BS NB National Bureau ofStandards NBl

QC

100

: NBS NBS NBS NBS NBS NBS NBS NBS N

•U5753

#1226

1986

C2

> NBS NBS NBS NBS NBS NBS NBS NBS N

> NBS NBS NBS NBS NBS NBS NBS NBSN
Tm he National Bureau of Standards' was established by an act of Congress on March 3, 1901. The

m Bureau's overall goal is to strengthen and advance the nation's science and technology and facilitate their effective application for public benefit. To this end, the Bureau conducts research and provides: (1) a

basis for the nation's physical measurement system, (2) scientific and technological services for industry and

government, (3) a technical basis for equity in trade, and (4) technical services to promote public safety.

The Bureau's technical work is performed by the National Measurement Laboratory, the National

Engineering Laboratory, the Institute for Computer Sciences and Technology, and the Institute for Materials

Science and Engineering

The National Measurement Laboratory

Provides the national system of physical and chemical measurement;

coordinates the system with measurement systems of other nations and

furnishes essential services leading to accurate and uniform physical and

chemical measurement throughout the Nation's scientific community, in- dustry, and commerce; provides advisory and research services to other Government agencies; conducts physical and chemical research; develops,

produces, and distributes Standard Reference Materials; and provides

calibration services. The Laboratory consists of the following centers:

• Basic Standards2 • Radiation Research

• Chemical Physics

• Analytical Chemistry

The National Engineering Laboratory

Provides technology and technical services to the public and private sectors to address national needs and to solve national problems; conducts research in engineering and applied science in support of these efforts; builds and maintains competence in the necessary disciplines required to carry out this research and technical service; develops engineering data and measurement

capabilities; provides engineering measurement traceability services; develops

test methods and proposes engineering standards and code changes; develops

and proposes new engineering practices; and develops and improves

mechanisms to transfer results of its research to the ultimate user. The

Laboratory consists of the following centers:

Applied Mathematics

Electronics and Electrical Engineering2 Manufacturing Engineering

Building Technology

Fire Research

Chemical Engineering2

The Institute for Computer Sciences and Technology

Conducts research and provides scientific and technical services to aid Federal agencies in the selection, acquisition, application, and use of computer technology to improve effectiveness and economy in Government

operations in accordance with Public Law 89-306 (40 U.S.C. 759), relevant Executive Orders, and other directives; carries out this mission by managing

the Federal Information Processing Standards Program, developing Federal ADP standards guidelines, and managing Federal participation in ADP

voluntary standardization activities; provides scientific and technological advisory services and assistance to Federal agencies; and provides the technical foundation for computer-related policies of the Federal Government. The Institute consists of the following centers:

Programming Science and

Technology

Computer Systems

Engineering

The Institute for Materials Science and Engineering

Conducts research and provides measurements, data, standards, reference

materials, quantitative understanding and other technical information fundamental to the processing, structure, properties and performance of materials; addresses the scientific basis for new advanced materials technologies; plans

research around cross-country scientific themes such as nondestructive

evaluation and phase diagram development; oversees Bureau-wide technical programs in nuclear reactor radiation research and nondestructive evaluation; and broadly disseminates generic technical information resulting from

its programs. The Institute consists of the following Divisions:

Ceramics

Fracture and Deformation 3 Polymers

Metallurgy

Reactor Radiation

'Headquarters and Laboratories at Gaithersburg, MD, unless otherwise noted; mailing address Gaithersburg, MD 20899.

2 Some divisions within the center are located at Boulder, CO 80303.
Research Information Center

National Bureau of Standards ^ ^ Gaithersburg, Maryland 20899

- US-7S3

no, fd-K.

NBS Technical Note 1226

Application of a Hard Sphere

Equation of State to Refrigerants

and Refrigerant Mixtures

Graham Morrison

Thermophysics Division Center for Chemical Engineering

Mark O. McLinden

Building Equipment Division Center for Building Technology

National Bureau of Standards

Gaithersburg, MD 20899

August 1986

Sponsored by:

Electric Power Research Institute Palo Alto, California 94303

W

U.S. Department of Commerce

Malcolm Baldrige, Secretary National Bureau of Standards

Ernest Ambler, Director
National Bureau of Standards U.S. Government Printing Office For sale by the Superintendent

Technical Note 1226 Washington: 1986 of Documents,

Natl. Bur. Stand. (U.S.), U.S. Government Printing Office,

Tech. Note 1226, Washington DC 20402

154 pages (Aug. 1986)

CODEN: NBTNAE
ABSTRACT

This note describes the application of the Carnahan-Starl ing-DeSantis equation

of state to halogenated hydrocarbon refrigerants and their mixtures. A complete and consistent set of thermodynamic functions is derived from the p-V-T

equation of state and the perfect (ideal) gas heat capacities. A thorough

discussion of reference states is included for both pure materials and their

mixtures. Although this model exhibits a critical point, it does not quantitatively represent properties in the critical region. Despite this limitation, this model can represent both liquid and gaseous mixtures away from

their own critical points, even at conditions near to and above the critical

points of their components.

Algorithms and FORTRAN routines for the use of this model are presented along

with the numerical coefficients for 11 pure refrigerants and 7 mixtures.

Routines for evaluating the coefficients from saturation data are included.

Several examples of the application of this equation of state are presented to

demonstrate its versatility. It is shown to predict the properties of pure

materials well and to describe the detailed features of mixtures, both phase

diagrams and thermodynamic properties. The average deviation from the tabulated saturation properties of the 11 pure refrigerants is 0.54 % for pressures, 0.09 % for liquid volumes, and 0.50 % for vapor volumes.

iii
ACKNOWLEDGMENTS

This work has been supported by a variety of organizations and individuals not

only through financial and material support but also advice and criticism,

extended discussion and guidance, and help in the laboratory and at the

workbench. We wish to thank the Electric Power Research Institute, and in

particular project managers, James Calm and Powell Joyner, for sponsoring the

work leading to this report. We thank managers William Levins and Fred

Creswick of Oak Ridge National Laboratory for their support of other aspects

of the refrigerant mixtures program. We also acknowledge the support of the

Building Equipment Division of the Center for Building Technology and the

Thermophysics Division of the Center for Chemical Engineering at NBS. We

particularly want to thank a few of the persons who have helped us in

uncountable ways in this work: David Didion, Peter Domanski, Reinhard

Radermacher, and Howard Ross of the Building Equipment Division*, Gaetane

Francis, Lester Haar, Tanya Neal, and Anneke Sengers of the Thermophysics

Division; and from outside NBS, Helen Connon, Max Klein, Uwe Schultz, and

Eugene Wright. Finally we thank Laurie Watkins, Gail Crum, and Mary Baier who

persevered our handwriting and a recalcitrant word processor in preparing the

many drafts of this Technical Note.

IV
TABLE OF CONTENTS

ABSTRACT iii

ACKNOWLEDGMENTS iv

LIST OF TABLES vii

LIST OF FIGURES viii

NOMENCLATURE x

1. INTRODUCTION 1

2. A COMPARISON OF EQUATIONS OF STATE 4

3. THE EQUATION OF STATE AND MIXTURES 12

4. THERMODYNAMIC FUNCTIONS AND REFERENCE STATES 17

Pore Fluids 17

The Thermodynamic Properties of Mixtures 29

5. CRITICAL BEHAVIOR AND THE EFFECT OF NEARBY CRITICAL POINTS 47

6. NUMERICAL IMPLEMENTATION 55

Determination of Equation of State Parameters 55

Determination of the Interaction Parameter 62

Calculation of Equation of State Parameters from Stored Coefficients. 63

Dew and Bubble Point Pressures 66

Enthalpy, Heat Capacity and Entropy 73

Iteration Technique 75

7. EXAMPLES OF THE APPLICATION OF THE CSD EQUATION OF STATE 78

Pure Refrigerants 78

Equilibrium Properties of the Mixture R13Bl/KL52a 91

Azeotropes and the Equilibrium Properties of RL2/R22 95

8. CONCLUDING REMARKS 103

REFERENCES 105
APPENDIX A: EQUATION OF STATE COEFFICIENTS FOR PURE REFRIGERANTS

AND MIXTURES 108

APPENDIX B: FORTRAN SUBROUTINES Ill

Subroutine FITAB 113

Subroutine FITF 117

Subroutine BCONST 119

Subroutine BUBLT 121

Subroutine ENTROP 127

Subroutine HCVCP 128

Block Data Routine BDESC 129

Subroutine ESPAR 132

Subroutine PLIMIT 134

Subroutine VIT 136

Example Run 13 8

VI
LIST OF TABLES

4.1 Thermodynamic Functions Arising from the Carnahan-Starl ing-DeSantis

Equation of State

6.1 Summary of Thermodynamic Property and Associated Subroutines

6.2 Refrigerants Currently Included in Data Base

7.1 Comparison of ASHRAE and Equation of State Values of Saturation

Pressure, Saturated Liquid and Vapor Volume, and Enthalpy of

Vaporization for 11 Pure Refrigerants

7.2 Comparison of Data Correlations for R113 Based on Saturation Values

of ASHRAE and Saturation Plus Superheated Vapor Data of Mastroianni,

et al.

A.1 Coefficients for Curve-fits of Pure Component a and b Parameters

A. 2 Coefficients for Curve-fit of Pure Component Perfect-gas Heat Capacity

A.3 Interaction Parameters for Refrigerant Mixtures

V 11
LIST OF FIGURES

Figure 4.1. Molar volume of a mixture at a fixed temperature and pressure

showing relationship between pure component and mixture molar

volumes and the partial molar volumes of each component.

Figure 4.2. Molar Gibbs free energy of a mixture at a fixed temperature and

pressure.

Figure 5.1. Molar volume of a mixture at a fixed temperature showing

potential error of estimating mixture properties by a linear

combination of the pure component values.

Figure 6.1. p-V-T behavior of a pure material as represented by the equation

of state; vapor-liquid equilibria criteria represented by

equality of shaded areas.

Figure 6.2. Flow chart for the calculation of dew or bubble point pressure.

Figure 6.3. Graphical representation of the secant and reguli-falsi

iterations to solve forip(6) = 0.

Figure 7.1 Percent difference between ASHRAE values of saturation pressure

and liquid and vapor volumes for R12 and values correlated by

equation of state.

Figure 7.2 Percent difference between ASHRAE values of saturation pressure

and liquid and vapor volumes for R22 and values correlated by

equation of state.

Figure 7.3 A comparison of the molar heat capacity at constant pressure for

saturated liquid R152a given by several sources.

Figure 7.4 Percent difference between values of saturation pressure and

saturated liquid volume correlated by equation of state and

experimental data of Mastroianni et al . (1978) for R113.

Figure 7.5 Percent Difference between values of pressure in the superheated

vapor region as a function of temperature and reduced volume

(V/V_) correlated by the equation of state and experimental data

of Mastroianni, et al . for R113.

Figure 7.6 Saturation pressure for the mixture R13B1/R152a at temperatures

from 15.9 to 55.0°C.

Figure 7.7 Saturated liquid volumes for the mixture B13B1/R152a at temperatures

from 15.9 to 55.0°C.

viii
Figure 7.8 The variations of the azeotropic composition of the R12/R22

mixture as reported by several authors and predicted by the

equation of state.

Figure 7.9 The heat of complete vaporization at 222.00 K for the mixture

R12/R22 as measured by Neilson and White (1959) and predicted by

the equation of state.

Figure 7.10 The heat of mixing for the E12/R22 mixture at 222.00 K as reported

by Neilson and White (1959) and predicted by the equation of state.

IX
NOMENCLATURE

a(a..) equation of state attraction parameter (between species i and

J)

a0' a

l' a2 coefficients for the temperature dependence of a (defined by

Eqn. 6.14)

A molar Helmholtz free energy

b(bj-) equation of state volume parameter (between species i and j)

b/j, bi, b^ coefficients for the temperature dependence of b (defined by

Eqn. 6.15)

B, C second and third virial coefficient

"P' P

C , (C ) molar heat capacity at constant pressure (for the perfect gas

state)

Cy molar heat capacity at constant volume

o Cq, ci, c« coefficients for the temperature dependence of C (defined by

Eqn. 6.16)

D, D* determinates related to mixture critical point (defined by Eqn.

5.8a,b)

E molar internal energy

f .. interaction parameter between species i and j

f/j, fi coefficients for the temperature dependence of f •• (defined by

Eqn. 6.17)

G molar Gibbs free energy

G total Gibbs free energy

H molar enthalpy

H

e molar excess enthalpy of mixing

I ionization potential

Z 1 iquid phase

n number of components in a mixture

n- number of moles of component i

N total number of data values
p pressure

p* arbitrary reference pressure

p° pressure of perfect gas state at specified temperature and

volume

l>low' p lower and upper bound on pressure for a physically meaningful

solution of equation of state

R gas constant

S molar entropy

T temperature

V molar volume

V total volume

xCx^ molar composition (of species i)

y b/4V

z compressibility factor

z

ft , z^ functions of chemical potential (defined by Eqns. 6.24 and

6.25)

Z , Z functions in accentric factor method

a coefficient of thermal expansion

P b/4

ft t isothermal compressibility

r

1 function defined by Eqn. 6.1

6, A change in a quantity

6

' x ' ith guess for iteration variable

0ipos;^ Q(neg) gUesses for 9 giving positive and negative values of function \\)

k^, ^ constants in Eqns. 5.4 and 5.5

p.. chemical potential of component i in phase j

a molecular diameter

y iteration function

(i) accentric factor

xi
t) £* w

t>* wv* wx we i&h tin6 factors for liquid volume, saturation pressure, vapor

volume, and vapor composition data respectively

Subscripts

c value at critical point, calculated value

e experimental value

i,j components i and j

£, 1 iq liquid phase

m mixture

sat saturation conditions

v, vap vapor phase

Superscripts — (overbar), partial molar quantity

perfect gas state

' temperature derivative

'

' second derivative with respect to temperature

( i) iteration index

xii
1. INTRODUCTION

This Technical Note describes the application of an equation of state based on

the hard-sphere fluid to an important class of industrial fluids, the halogenated hydrocarbon refrigerants and their mixtures. These compounds, including

several azeotropic mixtures, have been developed over the past several decades

and are now the dominant working fluids in refrigeration and heat pumping

equipment. Recently, considerable interest has arisen in the use of nonazeotropic mixtures as working fluids because of potential efficiency improvements and operational advantages. The lack of precise information on the

thermophysical behavior of these mixtures has hindered their commercial application. There is a need for a general property formulation that would be

applicable to the analysis of a variety of pure refrigerants and their

mixtures.

The purpose of this Note is to show how an equation of state founded on a

realistic physical model can be used to describe the thermodynamic properties

of a pure fluid or mixture with a minimum amount of experimental information.

Such an equation can be viewed on two levels. First, it represents a correlation of experimental data. But more than merely correlating isolated data, an

equation of state models the thermophysical behavior of a fluid. Such a model

has built into it the virtues of thermodynamic consistency, realistic limiting

behavior, and the power to predict unmeasured properties reasonably when there

is limited experimental information. Unlike simple correlating schemes, an

equation of state founded on a realistic physical model becomes an increasingly powerful tool as ones knowledge about a material increases. The authors
wish to emphasize, however, that there is no substitute for data*, any

predictive scheme, even one founded on a good physical model, can be

considered only a temporary substitute for experimental information.

This Note develops a particular equation of state—a modification of the

Carnahan-Starl ing hard sphere fluid that was first proposed by DeSantis, et

al. (1976):—and applies it to the representation of refrigerant thermodynamic

properties. We will refer to this as the Carnahan-Starl ing-DeSantis (CSD)

equation of state. Much of this development is general, and can be applied to

other similar equations of state. While only refrigerants are discussed here,

the CSD equation of state has been applied to other classes of fluids, such as

hydrocarbons and simple inorganic molecules (DeSantis, et al. 1976).

Section 2 of this Note briefly describes and compares classes of equations of

state. The motivation for choosing the hard sphere reference fluid is

explained. Section 3 presents mixing rules for the extension of the equation

of state to mixtures and introduces an empirical interaction parameter used to

account for deviations from the Lorentz mixing rules. Methods of predicting

the interaction parameter from molecular data are reviewed and the need for

experimental p-V-T-x data for its evaluation is emphasized. A complete set of

thermodynamic functions is developed in Section 4. These formulations use

only the p-V-T information of the equation of state and the heat capacity of a

perfect (ideal) gas. A discussion of reference states is also included.

Section 5 considers the behavior of the CSD and other equations of state near

the critical point. In Section 6, the numerical implementation of the equation of state is discussed, including the determination of the pure component
parameters and the mixture interaction coefficient. A set of computer

subroutines for carrying out these calculations is presented. Section 7

examines the ability of the equation of state to represent the thermodynamic

characteristics of halogenated hydrocarbon refrigerants and their mixtures.

Concluding remarks are presented in Section 8. Finally, the Appendices contain a tabulation of coefficients for use with the hard sphere equation of

state and a listing of the FORTRAN subroutines developed for its

implementation.
2. A COMPARISON OF EQUATIONS OF STATE

An equation of state represents the pressure—volume-temperature behavior of a

pure material or the p-V-T-x behavior of a mixture. When combined with heat

capacity information the complete thermodynamic behavior of a material

(including enthalpy, entropy, Gibbs free energy, etc.) is defined. For the

correlation of experimental data, an equation of state can represent all the

thermodynamic properties with strict consistency, an advantage not necessarily

retained when separate functions are used to correlate individual properties.

An equation of state also permits reasonable evaluation of properties for

which there are no measurements.

There are a number of classes of equations of state which are based on varying

combinations of theory and empiricism and with varying ranges of

applicability. The simplest of all equations of state is the perfect gas law:

£V = i (2.1)

RT

It describes the behavior of an assembly of point masses with no intermol ecular attraction or repulsion. It represents the properties of a gas in the

limit of zero pressure. Depending upon ones requirements, it is useful at

pressures less than 0.1 MPa*, at temperatures well above the critical temperature it can be useful at higher pressures.

The elegant form of the perfect gas law is retained in the virial equation of

state:
£V=1+I+C_+- • • (2.2)

RT V

y

2

where B and C are referred to as the second and third virial coefficients,

respectively. This equation was proposed in 1885 by Thiesen to describe the

behavior of gases and was shown by Kammerl ingh-Onnes (1901) to represent the

p-V-T properties of many gases well. Although originally empirical in origin,

Mayer (1937) showed that this expression arose naturally from a statistical

mechanical model. He was able to relate the experimental values of the virial

coefficients to the molecular properties of the material. Typically, the

virial equation is truncated at the second virial coefficient; such an expression is applicable at subcritical temperatures for pressures below about 1 MPa

and for increasing pressures at higher temperatures, but only for gases. A

major advantage of the virial equation of state is that second virial coefficients are known for a huge collection of materials (Dymond and Smith, 1969)

and can be estimated in the absence of measured values (Hirschf elder, et al. ,

1954).

The reduced equations of state are a second approach. An example is the

Pitzer (1957) accentric factor method:

Jj

= Z°(Tr ,Vr

) + (o Z

1 (T

r ,Vr

) (2.3)

The functions Z° and Z are universal functions of the reduced temperature and

volume derived by correlating the measured properties of many materials; the

accentric factor, w, is an empirical device to compensate for the non-
spherical nature of the molecule. Conformal solution theory (Henderson and

Leonard, 1971) and the TRAPP program recently developed at NBS (Ely and

Eanley, 1981) also fall into this category. A similar idea is the scaled

equation of state around the critical point (Stanley, 1971). These schemes are

based on the idea that by taking molecular size, shape, and attractions into

account, all materials can be compared in the same way. These forms, although

rooted in a theoretical model, depend on data for real materials to define the

functional relationships. Once these relationships have been defined by a set

of well-characterized reference fluids, only a minimal set of properties are

needed to estimate the scaled properties of other fluids. Such equations of

state are appropriate at all conditions, limited only by the knowledge of the

reference fluids and the ability of the scheme to compensate for the various

characteristics of real fluids.

One of the largest classes of equations of state, comprising a large fraction

of the industrially used expressions, are the modified van der Waals equations.

The original equation was proposed in 1873 by van der Waals:

'-.^t:.£

The first term compensates for the excluded volume in a bimolecular collision;

in the classical interpretation, b is four times the molecular volume. The

second term accounts for the attraction between molecules. Other equations in

this class, such as the Relich-Kwong-Soave (1972) equation, and the PengRobinson (1976) equation include modifications of the second term, and use it

to compensate for deficiencies in the first term in addition to accounting for
intermolecular attraction. These equations are cubic in volume, which is the

lowest order of equation able to represent both liquid and vapor behavior with

the same function. Although this class of equations of state has a liquid

branch, the molecular volume has been derived for a gas-like fluid and thus

there is no reason to expect that these equations will accurately produce the

dense fluid (e.g., liquid) conditions. Indeed, Henderson (1979) has shown

that such expressions are fundamentally flawed at high densities because the

packing problem has i.ot been properly addressed. Using an expression

involving a large number of arbitrary parameters can improve the fit in the

dense fluid region but has two disadvantages: first, it gives rise to calculation complexities and second, the physical meaning of the parameters is

lost.

The final approach consists of equations of state based on a theoretical fluid

which is used as a reference to which real fluid behavior is compared. In a

sense this approach is the theoretical analog to the empirical reduced

equation of state. The most important feature of a model fluid is not that it

describe a particular real fluid exactly but that it contain the proper treatment of the high density states. The Leonnard-Jones fluid (Verlet and

Levesque, 1967) is an example; it is composed of molecules that repel one

another at short distances and attract one another at moderate and large

distances. A second class of reference fluids are those characterized by a

hard-body repulsion. The state of development of this class of equations is

demonstrated by the fact that the most recent international steam tables

(Haar, et al. » 1984) are based on a hard convex-body equation of state.
A simple, well-characterized hard-body reference fluid, and the one on which

the equation of state described in this report is based, is the hard-sphere

fluid. Such a fluid exhibits an infinite repulsion force for a bimolecular

collision at some distance of closest approach; there is zero net intermolecular force at greater distances. There have been extensive computer

simulations of this fluid (Barker and Henderson, 1971) and statistical

mechanical treatments to describe it (Reiss, 1965). There are three analytical functions associated with the hard sphere fluid: the Percus-Yevisk

pressure equation (Lebowitz, 1964), the Percus-Yev ick compressibility equation

(Verlet, 1964), and the Carnahan-Starling (1969) pressure equation. The

Carnahan-Starl ing form shall be used in this discussion primarily because of

its familiarity to the engineering community:

£V „ 1 + y ± v

2 - y3

(2 5)

RT (1 - y)

3

where y = b/4V.

Although the Carnahan-Starling equation more closely represents the properties

of the hard sphere fluid at highly compressed, nearly close-packed states, all

the hard sphere fluid representations are equally good at lower, more

physically realistic densities.

This equation was derived as a closed form solution to an infinite geometric

series which was in turn based upon the virial expansion of Ree and Hoover

(1964). Carnahan and Starling demonstrated that eq. (2.5) accurately

describes rigid sphere behavior. In this equation, 'b' is a molecular volume.

8
Although it has a similar interpretation to the b in the van der Waals or

Redlick-Kwong-Soave equations, it will not have the same numerical value. In

the limit of large volumes this expression leads to a van der Waal s-1 ike term.

The converse, however, is not true; at small volumes the van der Waals

equation does not lead to eq. (2.5).

The hard-sphere equation is modified to compensate for long-range attractive

forces by the addition of a second, semi-empirical term to arrive at the

equation of state treated in this work:

pV _ 1 + v + v

2 - v

3 _ RT

' ,, _ v3 (1 -

y)

3

' RT(V + b)

(2.6)

This form, first suggested by DeSantis, et al. (1976), was proposed for the

prediction of multi-component, vapor-liquid equilibria and exhibited good

agreement for simple inorganic molecules, hydrocarbons, and their mixtures.

(The procedures for extending this equation to mixtures are discussed in the

next section.)

The physical significance of the parameters defined by the Carnahan-Starl ingDeSantis equation of state and their temperature dependence can be understood

by comparing the differences between the reference and real fluids. The real

fluid is not composed of molecules with a hard-sphere repulsive core. The

determined value of the 'b' parameter is an effective core that is a measure

of the average closest approach of molecules. In the real fluid, as the
temperature is raised, the average kinetic energy also increases and the

average distance of closest approach becomes smaller*, thus, the molecular size

parameter will decrease with increasing temperature because of the softness of

the short-range repulsion in the intermolecular potential. The parameter 'a'

has the dimensions of energy per mole times volume per mole*, in a spherical,

non-polar molecule it would be connected to the dispersion forces. The

attractive part of the real intermolecular potential can be highly directional

because of dipole moments and higher-order internal charge distributions. As

the temperature is raised, the strong coupling of these directional forces is

'washed out' and the average attractive force drops. The temperature dependence of 'a' is a reflection of the non-spherical character of the force

between molecules. In this work the temperature dependence of 'a' and 'b' are

represented by the following forms:

a = a^xpUjT + a2 T

2

) (2.7)

b = b

Q + b

x T + bjT2

(2.8)

These forms are those suggested by DeSantis, et al. (1976) with the addition

of the T

2 terms. The methods of determining 'a' and 'b' as functions of

temperature from experimental p-V-T data are discussed in Section 6.

This hard-sphere equation of state was selected for the present study of

refrigerant properties because it has the advantages of being a simple expression based on a good physical model. With a few exceptions, the species used

as refrigerants are small, nearly spherical molecules. These attributes in

10
turn allow only two physically significant parameters to represent both the

liquid and vapor phases. The equation defines a fundamental thermodynamic

relationship which, combined with perfect gas heat capacity information,

allows a complete set of consistent thermodynamic properties to be evaluated.

The equation of state can also be extended to mixtures.

11
3. THE EQUATION OF STATE AND MIXTURES

In the previous section, we have discussed the shortcomings of equations of

state with the van der Waals excluded volume term and how the equation of

state for the hard sphere fluid rectifies this serious physical flaw of the

van der Waals model. In this section, we shall discuss how the equation of

state for mixtures is produced.

The first approximation made in the application of the CSD equation of state

to mixtures is the assumption that no substantive change needs to be made to

represent the properties of a mixture. That is to say that there exists an

effective 'a' and 'b' such that an equation of form identical to the one used

for pure materials can be used for the mixture. This section will focus

primarily on the evaluation of such an effective 'a' and 'b'. A different

approach, accounting for the behavior of a mixture of hard spheres of

differing diameters, will be mentioned briefly at the close of this section.

In the simplest mixture models, the effective molecular parameters are defined

as follows:

= X L x

i x

j a

ij (3 * 1)

i=l j=l

and

- E £ *i*jt>ij < 3 - 2 > i=i j=i

12
When i = j, the values of a. . and b^ are those of the pure materials. The

values of &z- and b:: can be obtained if nearly any experimental property of

the mixture is known; examples will be given later in this discussion. The

motivation of a mixing rule scheme is to determine the values of a^. and hzwithout the assistance of measurements. The techniques used for predicting

the mixture parameters fall into at least two classes, those rooted in a

rigorous or at least a semi-rigorous physical model and those rooted in the

correlation of large sets of mixture data (Fredinslund et al. 1977)*, the

latter schemes are typically not applied to species with one carbon.

The origin of the size parameter, b12 > * s tne morc obvious of the two. The

closest distance of approach for two hard spheres is the sum of the two hard

sphere radii which is evaluated as follows:

b

iz = ——r*—

(3 - 3)

In this Note, Eqn. 3.3 is approximated as:

b

12

= (b

x + b

2 )/2 (3.4)

This approximation has been made primarily for mathematical simplicity. When

the values for 'b' are similar for the two components, the two mixing rule

schemes can barely be distinguished. Unless the molecular volumes differ by

more than 42 %, Eqn. 3.3 and Eqn. 3.4 will not differ by more than 1 %*, and,

hence, will not give rise to more than a 1/2 % difference in the value of b.

13
The physical origin of the mixing scheme for a, 2 has its origin in a variety

of molecular and quantum mechanical models. In the equation of state, the value

of 'a' arises from the attraction between the molecules. On a molecular

level, the interaction energy is typically expressed as a correction to the

geometric mean of the pure component energy force constants:

"12-- (1 " f12)(all a22

)1/2 (3 ' 5)

The correction factor, (1 - f- 2

) » attempts to account for differences in

polariz ability of the two species (the ease with which the two electron clouds

are deformed) and the closest distance of approach. Schemes for predicting

f-, are reviewed by Pesuit (1978), who considers them on several levels of

empiricism. What is impressive is the degree to which the interaction

parameter can be estimated for relatively complex molecules. Equally impressive however is the unexpected collapse of these estimation rules. Enobler

(1978) has shown that even for so simple mixtures as He/Ne, He/Ar, He/Kr and

He/Xe that an estimation of the interaction parameter such as the following is

far from adequate:

h + h (<7 U + a22 ) 6

I. is the ionization potential of component i; and a .., the molecular

diameter. The experimental value for (1 - f, 2

) * s 0.611 for He/Xe; its predicted value is 0.824. The conclusion that one would draw from the examples

given by Enobler and Pesuit are that the methods for estimating fj~ are not

14
very satisfactory and that a method appropriate for one class of mixtures may

not be satisfactory for another.

The thesis of these models is that all the information needed to describe the

interaction between unlike species can be derived from information about the

interaction between like molecules. This thesis neglects the possibility that

there can be interactions in mixtures that are either unimportant or do not

occur in the pure fluids; the converse is equally possible. Recently, Gubbins

and his colleagues have attempted to evaluate functions related to fj~ from

the mul t i pole- multi pole interactions between molecules (Gubbins and Twu,

1978). Wallis et al. (1984) have shown that this approach can be used effectively in describing mixtures of carbon dioxide and ethane. In general, one

should approach these estimation schemes with caution and in the spirit that

they can be used as guides when no other information exists. One should be

aware that, at present, the only quantitative estimation of a function such as

fj~ c & n arise from experimental information.

In this Note, the interaction parameter will be evaluated solely from

experimental mixture information. Almost any small set of data on the mixture

is sufficient to evaluate f -. 2

* Because there is but a single parameter to

evaluate for the mixture, one needs in principle only a single measurement;

however, the larger the experimental data set, the greater the certitude in

the evaluation of f.,. In section 6 of this Note, a description of the use of

bubble point information to evaluate f. 2

*- s given. One often finds cross

second virial coefficients, B, 2

» in the literature; such information can be

used to find iyy through the appropriate expansion of the equation of state.

15
One can also use gas phase data. Of the available data, liquid phase densities are perhaps the least reliable because a tiny difference between the

actual density and the equation of state density at the same pressure can lead

to large errors in the pressure.

The above discussion has dealt with finding an effective 'a' and 'b' for a

mixture for use with an expression developed for pure fluids. There are

equations of state that approximate the composition dependence for mixtures of

hard spheres (Mansoori and Leland, 1970) such as the following:

pV = 1 + (1 - 3A)y + (1 - 3B)?2 - (1 - Ov3

RT (1 _

y)

3

(3.7)

In this equation, the various terms in the Carnahan-Starl ing equation have

become functions of the averages of various powers of the hard sphere

diameters. When the two molecular volumes differ by a factor of 1.5, A is no

larger than 0.018, B no larger than .0089, and C no larger than 0.013. In the

limit of a pure material, A, B, and C, become zero and Eqn. 3.7 reduces to the

original Carnahan-Starl ing expression. In this work, we shall assume that the

departure from the Carnahan-Starl ing equation of state caused by the disparity

between the molecular sizes is unimportant and that this disparity can be

compensated for in the interaction parameter f.. «.

16
4. THERMODYNAMIC FUNCTIONS AND REFERENCE STATES

There are two themes found in each of the two parts of this section. The

first is a general description of some important thermodynamic derivatives.

In this Section on pure fluids, the derivation of the thermodynamic properties

from a p-V-T equation of state is discussed. In the part on mixtures, the

ideal mixture model is developed. The second theme confronts the problem of

reference materials and conditions (or reference states). In the section on

pure materials, the perfect gas appears as a necessary part of the derivation

and, through it, other reference states are constructed. In the section on

mixtures, the reference states of the components of the mixture are connected

to the thermodynamic properties of the mixture.

Pure Fluids

In the previous sections, we have discussed the physical origin of the hard

sphere equation of state and the merits of using such a model at high fluid

densities instead of a model where the excluded volume is expressed by the

van der Waals, l/(V-b), term. We begin here a review of the derivation of the

thermodynamic properties of a fluid from its p-V-T equation of state and show

the specific results for the Carnahan-Starling-DeSantis equation. In this

discussion we will find that the pressure equation of state does not include

all the information necessary to evaluate the thermodynamic properties completely. We shall resort to the properties of the perfect gas, spectroscopic

information, and relationships provided by statistical thermodynamics to form

a complete set of thermodynamic functions.

The pressure is connected to the Eelmholtz free energy and, through it, to all

the other thermodynamic properties by the following relation:

17
» - -ft) T

«•"

The Helmholtz free energy can be evaluated from the pressure by an integration

over volume. Several difficulties are encountered in this integration: first

the choice of appropriate limits for the integration*, second, the possibility

of divergences in the integral*, and, finally, the recognition that there will

be integration constants. The problems of integration limits and divergences

are closely connected*, they will be discussed together.

When the Helmholtz free energy is calculated, we would like to choose the

limits of the integration so that one of them is independent of the form of

the equation of state. This can be achieved by having one limit be V = eo

,

where all gases become perfect gases*. At that limit, however, the integral

diverges. This problem can be resolved by evaluating the difference between

the properties of the fluid described by the equation of state and those of

the perfect gas. As the volume approaches infinity, this difference becomes

zero sufficiently strongly that the integral remains finite. Thus,

P"P° - P - f

- -(L(A - A")) (4.2,

*The term 'perfect gas' is used rather than 'ideal gas' in order to preserve

the term 'ideal' for certain properties of mixtures.

18
and

A(V,T) - A° (V,T) = / (p - ^)dV (4.3)

V

Throughout this section, functions applicable to the Carnahan-StarlingDeSantis equation of state (such as those evaluated by Equations 4.2 and 4.3 )

are given in Table 4.1.

Before evaluating the remaining functions, let us review the features of the

Helmhol tz free energy calculation. First, in performing the integration, we

have resorted to a reference material, the perfect gas. It was chosen to

avoid a divergent integral. Although the perfect gas is a construct, it is a

convenient and useful material against which to compare the properties of any

other material—all gases become perfect gases in the V = °° limit and its

properties are known to the last detail. The second feature of the integration is a more general consequence of that operation. Since the integral was

over volume, both A(V,T) and A°(V,T) may be in part functions that are independent of volume. Because the difference between the two free energies is

uniquely defined, the temperature dependent and volume independent parts of

A(V, T) and A°(V, T) must be identical. Unless one is aware of these functions

which are constants with respect to the volume integration, other thermodynamic functions will not be complete.

Once the Helmhol tz free energy has been evaluated, all the other thermodynamic

functions arise by the usual operations

19
5 = - { 2A) (4.4a)

S - S°= -(L (a- A ))

< 4 ' 4b > \ 3T /y

and

E = A + TS (4.5a)

E - E° = (A - A°) + T(S - S°) (4.5b)

and

H - E + pV (4.6a)

H - H° = (E - E°) + (pV - RT) (4.6b)

and

G = A + pV (4.7a)

G - G° = (A - A°) + (pV - RT) (4.7b)

The values of the functions in Equations 4.4b, 4.5b, 4.6b, and 4.7b are given

in Table 4.1. The reference material for the functions enumerated in these

equations has, in all cases, been the perfect gas at the same temperature and

volume as the real fluid. For E and H, the perfect gas values depend only on

temperature, hence, volume is deleted as an appropriate argument for those

20
properties in Table 4.1. The reference condition for the Gibbs free energy is

often chosen terms of a fixed pressure. This choice arises because of the

role of the Gibbs free energy in determining phase equilibrium. Equation 4.7b

can be easily modified to account for a reference pressure, p*, as follows:

G(V,T) - G°(V,T) = G(V,T) - G°/p = ^, T

= G(V,T) - G°(p,T) + G°(p*,T) - G°(p*,T)

= G(V,T) - G°(p*,T) - RT ln^1 p*y

or

G(V,T) - G°(p*,T) = RT £n (%& ) + G(V,T) - G°(V,T) (4 * 8)

\P*V/

This function can be found in Table 4.1.

The heat capacity at constant volume is evaluated as follows:

cv

= (af)v

= T

(fi)v

(4 - 9a >

C - C° = (Z-(B - E°)) = t($-(S - S°)) (4.9b)

Earlier in this section, the need to recognize the existence of volume

independent functions as a consequence of the calculation of the Helmholtz

21
free energy from the pressure was noted. It is in the evaluation of the heat

capacities that the existence of these functions becomes most apparent. The

perfect gas has both internal energy and a heat capacity, yet both are

independent of the volume.

(dE°) = T

(dj>°) _

p

o = RT - RI = o Vav /T VaT A, v v

(4.9a)

/dc°\ _ aV _ aV . T

(d2 v°\ =0 (4.9b)

\3V /T

3V3T 3TdV \ dJ2 )y

The perfect gas values of E and Cy as well as H and C and parts of A, G, and

S can be evaluated from calorimetric measurements*, however, these quantities

are typically calculated through straight-forward statistical mechanical

relationships from infrared and Raman Spectra of low pressure gases.

The heat capacity at constant pressure can be calculated from the enthalpy or

entropy at constant pressure

c> (§?) - T (&-/

p

(4.10a)

however, because the equation of state is a function explicit in temperature

and volume, C i s more readily evaluated from Cy and temperature and volume

derivatives of the pressure:

22
s °v - T

(E)V/(K)T

<"*>

From the beginning of this derivation we have been forced to choose a

reference material. The initial choice, the perfect gas, was made to avoid a

divergence in the value of an integral. There may be reasons for choosing

other reference conditions. Some properties require no reference, properties

such as pressure, volume, temperature, and heat capacity can be measured

absolutely. Indeed, the reason for using the perfect gas reference for heat

capacity is that it can be determined quite accurately. Other properties have

only relative measures; a zero value has a significance only by convention.

For these properties, only changes can have any physical significance. Among

these are internal energy, enthalpy, and entropy. The two free energies

present a very special problem. A change of G when temperature is kept

constant requires knowing the volume. A change of A at constant temperature

requires the pressure, obviously there is no problem in either case. When the

temperature is changed, however, the value of S must be known. Depending upon

ones interpretation of the statement of the third law of thermodynamics, one

can debate whether the entropy has a known absolute value or not, surely its

value relative to a special state at K is known. Thus changes in G and A

cannot be calculated without carrying an unknown quantity T S(0 K) .

Fortunately, one uses G and A only in fixed temperature processes. We have

then at least three kinds of thermodynamic functions: those known absolutely,

those known only relative to a state at the same temperature; and those known

relative to a state at a different temperature.

23
One may have special reasons for choosing a particular reference condition.

One might choose a readily accessible physical state. For example, the reference state in chemical thermodynamics is often the element in its stable form

as found at 101.33 kPa and 25°C. In electrochemistry, all voltages are

typically compared to the standard hydrogen electrode. For gases, one might

choose the perfect gas at 101.33 rPa and 25° C. For refrigerants the reference

is typically the liquid at its saturation pressure at -40°C. There are, of

course, as many reference states as there are applications for thermodynamics',

they are all interconnected through thermodynamic relationships. Some of

these reference states carry the weight of an internationally accepted

definition and are often referred to as STANDARD STATES. The significance of

a STANDARD STATE is its widespread use and precise definition.

The next few paragraphs will be devoted to examples in which the enthalpy,

internal energy, and the entropy are calculated with respect to an arbitrary

reference temperature, T f, and reference pressure, Pre f. Since the pressure

is an explicit function of volume and since several states of a material may

have the same pressure, the description of the reference state is further

defined by indicating the volume parenthetically after the pressure.

The most effective way to go from one thermodynamic state to another is to

break, the path into steps for which a single variable is changed at a time and

for which the changes are readily calculated. For example, for the enthalpy

24
H(p(V),T) = H(p(V),T) - ff>(T)

+ ff»(T) - H»(Tref )

+ H°< Tref> - H<Pref< Vref)' W

+H<Pref< Vref>' Tref> <4;">

The terms in the first and third lines on the right hand side of Equation 4.11

can be evaluated according to the enthalpy expressions in Table 4.1. The

o second line involves the integral of C between the two temperatures, T and

T__«* Th© final term is the reference enthalpy, which may be set to zero by

re i convention or is more properly moved to the left hand side of the equation.

We have then

H(p(V),T) - H(pref <Vref ), Tref ) =

[H(p(V),T) - H»(T)]

T

r o + / C

p

dT

Tref

" [ H( Pref (Vref>' Tref> " H°<Tref )] «.12>

Internal energy can be evaluated in a similar fashion

E(V,T) - E(Vref , Tref ) =

[E(V,T) - E»(T)]

T

r o + 7 Cy dt

Tref

~[E(Vref , Tref) - E«(Tref )] (4.13)

25
As in the case of enthalpy, the first and third lines on the right hand side

of this expression are found in Table 4.1.

The derivation of the entropy is similar, however it has an additional term

because, unlike H° and E°, S° is not independent of volume.

S(V,T) = S(V,T) - S°(V,T)

+ S°(V,T) - S»(Vref ,T)

+ So<Vref ,T) - So(Vref , Tref )

+ S°< Vref Tref> " S< Vref Tref>

+ S(Vref Tref> < 4 ' 14>

In this expression, the first and fourth terms are gr.ven in Table 4.1, the

third term is the integral of c£(T)/T, and Vfe£ is the volume at the reference

state (e.g., for refrigerants, the saturation volume of the liquid at Tre £>

The second term is evaluated as follows:

oo Y /dS°\

S (V,T) - S (Vref ,T) = / [Jy-) dV (4.15)

VT * T ref

The Maxwell relation allows this derivative, (dS/dV)m to be evaluated in a

straightforward manner.

26
(-) = -

3_b\ (4.16)

f n3T V

For the perfect gas:

\3TA, V

The second term then becomes

S»(V,T) - S«(Vref ,T) = - / (^)dV =-Rln J

Vref V ref

The quantity that appears in the refrigerant property tables is as follows:

S(V,T) - S(Vref, Tref) =

[S(V,T) - S»(V,T)] - R ln( ^

Vref

T

+ / (

CV ) dT

Tref V T

" < S < Vref Tref> " S°<Vref> - Tref )] (4.17)

27
Table 4.1 Thermodynamic Functions Arising from the Perturbed CarnahanStarling Equation of State

pV _ 1 + y + y

2 - y3 a RT (1 _ y)

3 RT(V + b)

y = b/4V

A(V.T) - A°(V, T) =-iln (£+-*-) + 4RTP + —*^

b \ V J (V- B> (y _

f p = b/4

G(T,p(V)) - G°(p*,T) = RT m (*£-) - » I n(

W + b

) + —^

(8V2 - 9Vp + 3p2 ) *- Vp*V/ b I V ) (y _

p) 3 H V +

(V - B) 2

lid.pOD.^) = GCT.pm.Xi) + (1 - ii)(|^

Xi

T.p

= ^(,..I., 4

) + RTUW + RTP(4 V - 3 P ) + RI^j(4V2 - 2VP)

+ ab^ /v^ 1 W (V - P) 2 (V - 6)3 b

2 V V

ab i +

2 *i a i + 2 *i a ii

£ / V

b(V + b) b \V + b

Note: the above expression for chemical potential applies only to a binary mixture.

S(V T) - S°(V T) = a

' b ~ ab ' I JW ± b \ +

ab ' - RB(4V - 3 6) _ RTB'UV2 - 2VB)

b

2 V V y b(V + b) (y _ p)

2 (v _

p)

3

E(V T) - E°(T) = a

' bT " ab ' T " ab : n(W + b

) i ab'T _ RT2 B'(4V2 - 2VB)

b

2 V V / b(V + b) (v _

p)

3

H(T p(V)) - H°(T) = ab'T - a'bT - ab „ p/V + b\ , ab'T- ab + RT(4V2 - 2Vp? (p - p'T)

b

2 V b ) b(V + b) (y _ p)

3

Cv(V,T) - c£(T) = gRI^B^CVB - 2V2

) + 2RTV((B"T + 2B')(B - 2V) + B' 2

T) _ Tab'

(V - p) 4 (V - P)

3 b(V + b) 2

+ T(ab"b + 2a'b'b + 2ab' 2 ) _ ( a"b2 T - 2a'b'bT -i- 2ab' 2 T - ab"bT) f ^V ± b

b

2 (V + b)

c

p =cv- t (S)

2

/(Si)

V.

x

T.

x

NOTE: primed quantities indicate a temperature derivative; double primed

quantities are second derivatives with respect to temperature.

28
The Thermodynamic Properties of Mixtures

The above discussion presented in detail for a pure substance: 1) the

calculation of thermodynamic properties from an equation of state, 2) the

perfect gas state, 3) reference states.

In this section, the discussion will be expanded to include the thermodynamic

properties of mixtures. The discussion of mixtures in many ways will parallel

that for pure materials but will differ in two important ways because of the

extra degree(s) of freedom accessible to mixtures, the compositions. The first

topic to be discussed will be reference states for mixtures. Integral to that

discussion w ill be the description of a special class of mixture, the ideal

mixture. The second topic that will be covered is the class of thermodynamic

properties, the partial molar properties, that arise from the derivatives of

extensive properties with respect to changes in the composition. In this

class of properties, the chemical potential will be needed to evaluate conditions of phase equilibrium.

If the components of a mixture did not tend to separate during a phase transition, one could define reference states for mixtures in the same way that such

states are defined for pure materials. To use refrigerant materials as an

example, one could choose the liquid phase at its bubble pressure at -40° C.

The incipient vapor in equilibrium with that liquid typically has a different

composition. If that vapor were condensed and its reference point were

defined as the liquid at its own bubble pressure at -40°C, one would

implicitly have changed the reference point of the parent liquid phase; chaos

would follow! The following paragraphs discuss a solution to this dilemma.

29
The properties of pure materials were initially tied to the properties of the

perfect gas. Through that connection, the properties of the material at any

two arbitrary states are automatically connected. The properties of mixtures

will be connected to the properties of their constituent components through

the ideal mixture. In the discussion of the ideal mixture, we shall find that

the ideal mixture of perfect gases will be particularly useful.

The properties of the ideal mixture can be described in terms of the

following simple experiment: the mixing of the pure constituents all at the

same temperature and pressure, to form a mixture at the same temperature and

pressure. Regardless of the composition, temperature, or pressure, one would

find for an ideal mixture that first, the final volume of the mixture would

equal the sum of the volumes of the unmixed components. That is,

V.CT.p.Cxi)) -2>iVi(T,p) = o (4.18)

i

(Where V^ is the molar volume of component i or the mixture and i- is the mole

fraction of the respective component in the mixture.)

Secondly, the mixture would have to be neither heated nor cooled to bring it

to the same temperature as the unmixed components*, the enthalpy of the mixture

would be the same as the sum of the enthalpies of its constituent pure

components:

yT.p.lxj}) -L»iHi(T,p) = (4.19)

30
By multiplying Eqn. 4.18 by the pressure and then subtracting the result from

Eqn. 4.19 we have a similar relationship for the internal energy, E.

Em(T,Vm

(p),{xi }) -X^x.EidW^p)) = (4.20)

i

When the temperature derivatives at constant pressure are taken of Eqns. 4.18

and 4.19, one finds similar relations for the coefficients of thermal

expansion (weighted by the molar volume) and the heat capacity at constant

pressure.

VT.p.Uil) Op^T.MxiH-I^iViCT.p) opi (T,p) = (4.21)

(Where a is the coefficient of thermal expansion of the respective component

or mixture, I (^) .) V VaTA

C

pm (T'P (V) ' {x i }) -Ex

i Cpi (T' p(V)) - < 4 - 22 > i

(Where C is the molar heat capacity at constant pressure of the respective

component or mixture.)

When the pressure derivative at constant temperature is taken of Eqn. 4.18, a

linear relation for the isothermal compressibilities (weighted by the molar

volumes) can be derived:

31
Vm(T,V.Ui })fiTm(T,p.Ui }) -£x iVi(T,p)PTi (T,p) =0 (4.23)

(Where pj is the isothermal compressibility of the respective component or

mixture, " ff (&) .)

v \ap/T

A similar relationship can be derived for the molar heat capacity at constant

volume. The simple linear relationship is valid, however, only under very

restricted conditions. By using the following relation between C and C : v p

C

y

= C

p

+ TV aJ/PT

(4.24)

the first term on the right-hand- side is linear in composition for ideal

mixtures (Eqn. 4.22). The second term can be made linear with respect to the

composition only if both the numerator and the denominator are insensitive to

the composition. Except for a truly unusual liquid mixture, the only material

that satisfies such a requirement is the perfect gas. The perfect gas thus

plays a central role in connecting the properties of the mixture to those of

its constituent components. For perfect gas mixtures.

CVm(T,V(p),{x i }) "E x i CVi< T' V i<P) > = ° (4 ' 25)

(Where Cy is the molar heat capacity at constant volume for the respective

component or mixture.)

32
There is one further property of the ideal mixture, the entropy of mixing.

Sm

(T,p(V),{xi }) -£*iSi(T,p(V)) =-R2xiln(x i

) (4.26)

i i

By comparing Eqns. 4.19 and 4.26, the Gibbs free energy associated with mixing

is as follows:

Gm(T,p(V),{x i }) = H^LpCY).^}) - TSm

(T,p(V),{x

i })

or

Gm

(T,p(V),{x i }) -^XiG^T^pCV)) - RT^XilnCx.) (4.27)

i i

The definition of the properties of the ideal mixture allow the real mixture

properties and the real pure component properties to be connected, whether it

be the pure components at the same temperature and pressure or the pure

components under some other conditions, in particular, their respective

reference states. The details for determining the mixture equation of state

have been discussed in Section 3 of this Technical Note. Once the form of the

equation of state has been established, all the thermodynamic properties of

the fluid mixture with reference to the respective properties of a perfect

gas mixture with the same composition can be evaluated by using the formalism

described earlier in this Section. Two examples follow: the evaluation of

the enthalpy of the mixture first with respect to the enthalpy of the pure

components at the same temperature and pressure as the mixture and then with

respect to a reference state of the pure materials. For the sake of simpl i33
city, a binary mixture will be discussed; this example may be expanded to any

number of components.

H^T.pW.X!) = H^pCV) ,x

x

) - H^Xj)

+ ^(T^p - XjHjCT) - (1 - x

x)^(T)

+ x^Hjm - H

1 (T,p(V1 ))]

+ (1 - Xl )[^(T) - H

2 (T,p(V2

))]

+ x

1 H1 (T,p(V1

)) + (1 - x1 )B2 (T,p(V2

)) (4.28)

By inspection, one can see that this expression reduces to the identity

Hm

(T,p(V) ,xj) = H^TxpCV) ,xj) . The first, third, and fourth line of the

right-hand- side of Eqn. 4.28 can be evaluated using the functions in Table 4.1

and the appropriate values for the equation of state parameters, a and b. The

second line equals zero because a perfect gas mixture is ideal. The final

line is the weighted sum of the enthalpies of the pure fluids at the same

temperature and pressure as the mixture. The sum of the first, third and

fourth lines in the right-hand- side of Eqn. 4.28 is the enthalpy of mixing of

the components '1' and '2' at fixed T and p. A simple expansion of the

calculation that led to Eqn. 4.28 produces the value of the enthalpy of the

mixture with respect to the reference pressure and temperature of its constituent components, as follows:

34
Hm

(T,p(V), Xl ) = Hm

(T,p(V),x1

) - ^(T.Xj)

+ H^er.xp - x^m - (i - x^hJct)

+ x^Hjtt) - I?(Tlf ,0f)

+ El (Tl,ref'P) " ^^.ref'Pl.ref*

+ H

l< Tl,ref> " %<Tl,ref'Pl,ref<V1 , re£ ))]

+ (1 - x^) [terms for component 2 similar to those for component 1]

+ x

l Hl (Tl,ref'Pl,ref< Vl,ref>>

+ <!" x

l> H2< Tl,ref'P2,ref< V2,ref>> < 4 ' 29>

In this expression, the first and second lines are identical to those in Eqn.

4.28; the third line is the integral of the perfect gas molar heat capacity at

constant pressure between the temperature of interest and the reference

temperature for component 1. The next line equals zero because of the

properties of the perfect gas. The fifth line can be evaluated using Table

4.1. A similar set of terms exists for the second component. The final two

terms on the right hand side of Eqn. 4.29 represent the pure components in

their respective reference states. The enthalpy of the mixture then becomes

the following:

35
H

m

(T,p(V), Xl ) = H^T^V),

x

x

) - H£(T, Xl )

T

.0 + Xl [ / C

pl

(T) dT + Bl°(T1<ref ) - ^(Tltiet ,p1>tet (Vlttef )n

Tl,ref

+ (1 - x^) [terms for component 2 similar to those for component 1]

+ ZlHl^l.ref'Pl.ref^ref^

+ (1 - ^

1 )B2 (T2Tef ,v2fTef (\2ief )) (4.30)

The last two terms cannot be determined absolutely. Setting such terms equal

to zero is an arbitrary assignment of the reference state. More properly,

these terms should be moved to the left-hand side of the equation.

The entropy for the mixture is evaluated as follows:

Sm(T,p(V), Xl ) = Sm

(T,p(V),x1

) - S°(T,V, Xl )

+ S2(T,V, Xl ) - Xl Sf(T,V) - (1 - X2 )S£(T,V)

+ Xl [S?(T,V) - S°(Tlref,V)

+ S

l°< Tl,ref' V> " Sl< Tl,ref'Vi, ref>

(continued on next page)

36
- S

l (Tl,ref' Pl,ref< Vl,ref>> "

S1°< T

1 , ref

Vl , ref >

>

+ (1 - x,) [terms for component 2 similar to those for component 1]

+ i

l S

l^ Tl,ref'Pl,ref< Vl,ref ) >

+ (1 - *2

> S2 (T2,ref' P2,ref< V2,ref>> (4 ' 31)

The path that is followed in the calculation is the same as for the enthalpy.

The first, and fifth lines can be evaluated from Table 4.1 by substituting the

appropriate values for the molecular parameters into the equation of state.

The second line can be evaluated from Eqn. 4.26. The third line is evaluated

o by integrating Cy/T over the temperature range T £ to T. The fourth line is

evaluated by integrating (dS/dV)TM over the volume range and using the Maxwell

relation that OS/3V) T

= OP/3T)y, which equals R/V for a perfect gas.

Equation 4.31 then becomes the following:

S

ffi (T,p(V), Xl ) = S^T.pfV),:^) - S°(T,V,x1

)

- R(xx n(xx

) + (1 - x

1 )£n(l - x

x

)

)

T

+ x^f £v dT + E ln(V/V1 re£ ))

Tl,ref T

" *<*..£ Pl,ref < Vl,ref>> "

S°< Tref Vref>l

+ (1 - x1

) [terms for component 2 similar to those for component 1]

+

*1 h (Tl,ref Pl,ref < Vl,ref>>

+ (1 - x

x

) S

2 (T2>ref , P2jiet (V2jref)) (4.32)

37
The evaluation of any mixture property proceeds by a similar method. First

evaluate the difference between the real mixture and the ideal-perfect gas

mixture. When the components are separated, the derived nature of the ideal

mixture prescribes the associated change in the property. After the

separation has been made, each component can be treated as if it were a pure

material.

The additional degree(s) of freedom with mixtures, the compositions, gives

rise to a set of properties that typically have a trivial meaning in pure

materials, the partial molar properties. Earlier in this section, the

experiment that led to the definition of the ideal mixture was discussed. Let

us discuss an experiment performed under similar conditions for which the

outcome is different. Suppose one had a mixture of two materials, A and B, at

a particular temperature and pressure. Let us further suppose that a small

amount of B were added to the mixture while keeping at T and p constant. If

one measured the ensuing change in volume and attributed that change to the

'effective' volume of the added component in the mixture, that effective molar

volume would be calculated as follows:

VB

= [V(T,p,nA

,nB

+ 6nB

) - V(T, p, nA, n

fi

) ] /6nB

(4.33)

where the bold face V indicates a total (rather than a molar) volume, and nrefers to the total number of moles of the respective components. In the

limit of adding an infinitesimal amount of B, the effective molar volume

described in Eqn. 4.33 becomes the value of the following partial derivative:

3S
VB

= (6V(T,p,nA

,i!B)/anB

) (4.34)

T, p, n^

This partial derivative is the definition of the 'partial molar volume'.

(Partial molar quantities are indicated by an overbar.) One could have

evaluated the partial molar volume of component A, V», in terms of the partial

derivative of the volume with respect to the amount of A. One can define

similar partial molar quantities for any extensive property of the system. In

this example, one can be comforted by the explanation that the partial molar

volumes are effective volumes. Partial molar quantities often do not retain

this simple interpretation; thus one should keep in mind the formal mathematical definition of these quantities.

In the previous paragraph we have discussed an operational definition for the

partial molar volume. Typically, it is experimentally determined in a different way, by noting the variation of the molar volume with composition. In the

following paragraphs, the connection between the molar volume and the partial

molar volume will be discussed and an example of an important class of

thermodynamic relationships will be encountered.

Because volume is an extensive property, more properly a first order

homogeneous function of the amount, one can write the following relation

between the total volume and the molar volume of a material:

va,?,^,^) = nv(T, P, J^ir)

= nV(T,p, Xl ) (4.35)

39
where n = n-.+ iu. The partial molar volume of the two components can be

calculated as follows:

1 l,p, n2

= V (T.p.iij) + n( — > ' - 3xl/T,p \ an

p \"»1/T,p,ll2 l/T,p,

= V + (1-x

l

) (dx1

) T)

(where x* = n^/ (n^+v^))

3V

V2 (T,p, Xl ) =(^"

2

/ T,p,n

1

= V (T.p.nj) + n f — l '-^

Sxi/T.pVani/T,?,^

= V

/av\

(4.36a)

(4.36b)

One could expand this definition to as many components as one wished. The

resulting partial molar volume would be

i i x 'P' xk^i

40
where, 8 — is the Kronecker delta (8"=! when i=j and 8,.=0 when i^j) , and xis the set of independent mole fractions. Recall that the number with

independent mole fractions is one fewer than the number of components since

the sum of the mole fractions is unity.

Figure 4.1 shows the molar volume plotted versus mole fraction for a binary

mixture. Eqns 4.36a and b have an immediate geometric interpretation; V, and

V2 are respectively the Xi =l and x*=0 intercepts of the tangent to V at X...

It immediately follows that

V = x^ + (1 - x

± )y2 (4.37)

which is consistent with the notion that the partial molar volumes are, in a

sense, effective volumes for the components in the solution. If one moves to

mixtures with more than two components, the partial molar volume becomes the

intersection of a many dimensional plane tangent to the V-composition surface

at the composition of the mixture with the respective pure component axes.

By evaluating the composition derivative of the molar volume, another important

relationship can be derived.

Wt,p

1 2

* Wt,p

* Vs Vt, (4.38)

41
LU

O

>

<

x B,i 1

COMPOSITION (mole fraction component B)

Figure 4.1. Molar volume of a mixture at a fixed temperature and pressure

showing relationship between pure component and mixture molar

volumes and the partial molar volumes of each component.

42
\ The difference V< - V2 can be evaluated by subtracting Eqn. 4.36b from Eqn.

4.36a; however, that difference equals ( —

\3Xl

7 T,p

One concludes then that

avi

x, "1 + (1-x,) 2^2

a Wt, p

^

9x1

dV, = (4.39)

This relation, one of a set of relations known as 'Gibbs-Duhem' relations, can

be justified geometrically by noting that when the tangent is rolled along the

V - x* curve that changes in the right and left intercepts are related by

construction of similar triangles as follows:

^2 = SVl

1-x-,

For the purposes of phase equilibrium in mixtures, the most important partial

molar quantity is the chemical potential, which can be defined in the

following ways:

aA

= OE<S.V,{n i »/8nA

) 8-Vfni ^ A

OH(S.pW.{n1 ))/8nA

) gfp A

OA(T,V,{n.})/anA

) T)V#ni ^ A

(8«T.p(V),{ni })/a.A

) T|P(ai<lA

(4.40)

43
The chemical potential for the Carnahan-Starl ing-DeSantis equation of state is

most easily evaluated by differentiating the Hemholtz free energy; the

resulting expression is given in Table 4.1 for a binary mixture.

Whenever two or more phases are in equilibrium, not only must the temperature

and pressures of all phases be the same, but also the chemical potential of

each component must be equal in all phases. This result is a consequence of

the thermodynamic stability requirement that, at a fixed temperature and

pressure, the Gibbs free energy of a system will seek a minimum. How the

phase equilibrium arises can be seen in Figure 4.2, where the molar Gibbs free

energies for the two phase forms intersect at XqJ however, between ii and iy,

the mixture can have yet a lower Gibbs free energy by forming two phases. The

free energy of the two phase system lies on the common tangent to the free

energy curves for each phase form. In an argument similar to the one for the

volumes, the xR

= and xB

= 1 intercepts represent the respective partial

molar Gibbs free energies (or chemical potentials). With the common tangent

construction, one can readily see that the respective chemical potentials of

each of the components will be the same in both phases. Were one dealing with

a three component system, there would be a common tangent plane construction^

for systems with more than three components, there would be a common tangent

hyperplane, a problem that is more easily contemplated mathematically than

visual ized.

In this section, we have reviewed several important ideas. First, we found

that the properties of the perfect gas state for both pure materials and for

mixtures was an essential device for connecting one state of a material with

44
>-o

DC

HIz

LU

111 LU

CC

LLCO

00

COO

o:<

Chemical

potential of

more volatile

component

Two-phase

region

Vapor

composition

*1

Liquid

composition

X2

Chemical

potential of

less volatile

component

COMPOSITION, xb (mol frac less volatile component)

Figure 4.2 Molar Gibbs free energy of a mixture at a fixed temperature and

pressure. The stable phase is the one having the lower Gibbs

free energy; for compositions lying between a common tangent to

the liquid and vapor branches a two-phase mixture posses a lower free energy than either phase.

45
another. Second, we found that the p-V-T equation of state alone was not

adequate to describe the temperature dependence of a material; knowledge of

the perfect gas heat capacities was essential. Finally, we found that the

properties of a mixture could be connected to the properties of its constituent components through the ideal-perfect gas mixture and that, by using this

device, the reference states of the constituent components of a mixture define

the values of the thermodymamic properties of the mixture itself.

46
5. CRITICAL BEHAVIOR AND THE EFFECT OF NEARBY CRITICAL POINTS

Through the past few sections, we have developed a fluid model based upon a

hard sphere reference fluid. By using this model, the thermodynamic properties of fluids, both pure components and mixtures, can be described by having

only a modest set of experimental information. In this development we have

avoided critical points, one phenomenon that cannot be described quantitatively by this or any of the commonly used industrial equations of state.

First, let us note that this model has a critical point. A single component

fluid is at its critical point when the first and second derivatives of

pressure with respect to volume are zero:

cLe\ = (&A =

(5 ' 1}

av/T \av2

/T

For the Carnahan-Starl ing-DeSantis equation of state this yields:

(ta\m o = SI T-v4 - bv3 - ^£ + iiv - blJ

{dW{ v2/v b\4 L

4 " 256. riv-j

+ ai2Vjt_bl_ (5 2)

V

2 (V + b) 2

e 2 p\_ ft ... RT Lv5 . IbV* . Sb2 V

3 _ 5b3 V

2 . 5b4 V _ b£. |= = ^% v

3

( v - j;

.^5 + TbVl + 5t£v

L 2 2 4

' «/ v\5 2 4 16 128 512

+ a(-6V2 - 6Vb - 2b2

) (5#3)

V

3 (V + b) 3

where all quantities are evaluated at the critical point.

47
Equations (5.2) and (5.3) could in principle be solved to yield the critical

temperature and volume in terms of the a and b parameters*, the critical

pressure would then be given directly by the equation of state. The

complexity of these expressions, however, makes an analytical solution

impossible. Thus we postulate:

Vc = X1 b

c (5.4)

where X, is a constant. Such a relationship holds for the van der Waalslike equations of state. Substitution of Eqn. (5.4) into Eqn. (5.2) yields

for the critical temperature:

_

a c^2 .. a T

c^2 /* o

« " m7 " E£b7 (5 - 5)

where X n is a complex function of X, and thus also a constant. Equation (5.3)

is satisfied using the above expressions for V. and T_, thus confirming the

assumed relationship for V (Eqn. 5.4). A numerical solution of the above

c system of equations yields:

V„ = 3.006818 b„ (5.6)

0.227329a c T„ = e (5.7)

Rbc

The critical pressure is given by

48
0.315714 RTC

Pc - y (5 * 8>

The constant 0.316 in Eqn. (5.8) is the value of the compressibility factor at

the critical point. This compares to a critical compressibility of 0.375 for

the van der Waals equation and 1/3 for the Redl ich-Kwong-Soave equation of

state. A wide variety of organic fluids have a Z_ of approximately 0.27-0.29

(Reid, et al., 1977) , thus none of these equations of state can quantitatively

predict near critical behavior. Were we to use data far from the critical

point (Sengers et al., 1981) to evaluate the parameters in the equation, we

would find the following discrepencies between the measured and predicted

critical properties: first, the predicted and measured critical points would

not coincide within experimental uncertainty; second, the predicted values of

all the extensive properties—volumes, enthalpies, entropies, etc.—would

differ from the measured values in fundamental ways; finally, the values of

the thermodynamic response functions—-heat capacities at both constant pressure and vol Time, the isothermal compressibility, and the thermal expansion

coefficient—would all diverge more strongly near the critical point than

those respective properties predicted by the equation of state.

One is tempted to resolve these differences by forcing the equation of state

to match the critical behavior. Such coercion does not represent an

acceptable solution to what is a fundamental physical problem. Forcing the

equation of state to match the critical behavior affects the temperature

dependence of the molecular parameters a and b; thus, states far from the

critical density—and critical point—but near the critical temperature would

'sense' the critical point in a physically unrealistic way. The alteration of

49
an equation of state to produce both proper near-critical and far-fromcritical behavior is a major task (Woolley, 1983) (Fox, 1983) and will not be

discussed here. As long as one is not operating too near the critical point,

a somewhat subjective criteria that depends upon the property being

considered, 'classical' equations of state, such as the one discussed in this

paper, can describe the properties of a fluid quite accurately.

The criteria for critical points in binary mixtures are the following:

D =

3^A

3 V

2

T,x

3jfA_

3x3V

(Hi)

\3x3V/

d±k

3x2

V,T

- (5.9a)

and

D* =

3D

3V

3D

3x

T,x

T,V

3_2a\

3x3V/

3 2 A

3x2

V,T

= (5.9b)

or, by straightforward thermodynamic transformation

3 2 G

3x

= ®M [ ^

N

) - T.P 3x-

(5.10)

T,p

50
Although the criteria for mixtures appear to be different than for pure

materials, the physical and thermodynamic reasons for them are the same. The

differences are only the consequences of the additional degrees of freedom, the

compositions, available in mixtures and not in pure materials. Indeed, the

pure material criteria expressed in Eqn. 5.1, are a special case of the

criteria in Eqn. 5.9a and b. The essence of Eqns. 5.1, 5.9a and 5.9b is the

thermodynamic stability requirement arising from the second law; for a system

to be stable any perturbation of the system at fixed volume and temperature

must cause the Helmholtz free energy to rise. The first of the criteria in

Eqn. 5.1 and that in Eqn. 5.9a define the boundary between a locally stable

and unstable thermodynamic state. A critical point is a special place on this

boundary that requires that any perturbation will lead into a stable state; at

other points on the boundary there will be perturbations that lead into

unstable states. The important fact to note however is that the critical

point for a mixture is not given by the same function of a and b as the pure

materials.

Even though the equation of state cannot quantitatively predict the near

critical region, applying such an expression to mixtures is superior to the

usual mixing rules. The simplest mixture approximation, and one that is often

applied, is that of an ideal mixture where the mixture property (e.g., molar

volume) is given by the mole fraction weighted average of the pure component

v al ue s

:

Vm = *AVA

+ *BVB

< 5 ' n>

51
Such mixing rules can work well when both, pure components are well below or

above their critical points and if no strong intermolecular interactions (such

as hydrogen bonding) are present. These rules can collapse, however, when one

of the components is near its critical point, even though the mixture itself

may be an ordinary fluid (Morrison, 1985).

Figure 5.1 shows how errors can arise with simple mixing rules. In this

illustration, which plots molar liquid volume against composition, the

temperature is near, but below, the critical temperature of component A and

well below the critical temperature of component B in the representative

mixture A/B. The curve ab is the locus of saturated liquid volumes. The

curve af is the isobar at the saturation pressure of component A, The strong

curvature of ab and af near the pure A side of the figure are indicative of

the near critical state of component A The ideal mixture approximation would

give a straight line between a and b or, more correctly between a and f. In

both cases, the straight line deviates from the corresponding curves by about

15 percent in the mid composition range for this example. Of course, by the

nature of the mixing rule, the curve and line must coincide at pure A and B. A

well-chosen equation of state, however, would follow the proper curved path,

although perhaps deviating from the actual volume as the mixture approaches

its critical point at high concentrations of A.

In summary, so long as one stays away from the critical point of the mixture,

classical equations of state can be made to work quite adequately for the

quantitative evaluation of the properties of the mixture. As noted before,

one's only recourse should one of the components of this mixture be near to or

52
120

110

o

E

CO

E

o

LU

O

>

DC

<

100

90

80

Mixture volume

given by linear

combination

estimate

Error in

volume

\

predicted

by linear

combination

Mixture \

volume at v saturation

Mixture volume at

constant pressure

^

(P = P, Component A) . \

COMPOSITION (mol frac)

Figure 5.1. Molar volume of a mixture at a fixed temperature showing

error of estimating mixture properties by a linear

combination of the pure component values.

53
above its critical point is the use of an equation of state; all schemes used

for ordinary liquids break down completely under such conditions. As indicated previously, the problem is more fundamental than the absence of liquid

data for the respective component above its critical point; the presence of a

nearby critical point causes properties that otherwise change slowly and

smoothly with changes in composition, temperature, or pressure, to change

rapidly. As a consequence any linear interpolation scheme ceases to serve as

a quantitative way of estimating a property.

54
6. NUMERICAL IMPLEMENTATION

The algorithms which express the p-V-T and thermodynamic relationships

presented in the preceeding sections have been implemented in a set of

computer subroutines written in standard ASCII FORTRAN (FORTRAN 77). These

are summarized in Table 6.1 and fall into four general categories: 1) subroutines to determine the a and b parameters and the mixing coefficient, f, from

saturation data, 2) routines which store coefficients for the a, b and C°

expressions for pure components and calculate these parameters as a function

of temperature and composition for a mixture, 3) property routines for the

calculation of enthalpy, heat capacity, entropy, specific volume, and saturation pressure and 4) auxiliary routines which are referenced by the other

property routines. The subroutines are interdependent with one calling

another. This section details the methods of solution employed by these

routines. The required inputs and resulting outputs for each routine are

summarized in Table 6.1 and discussed in detail in Appendix B.

Determination of Equation of State Parameters

The a and b parameters in the equation of state must be determined as a

function of temperature for each of the pure materials of interest. Data are

required for this determination. Laboratory measurements or original literature values are preferablej tabular values (e.g., from handbooks) may also be

used but are less desirable because considerable data smoothing has taken

place to generate the tables. The most commonly available and generally most

reliable vapor-liquid equilibria data are the saturated liquid and vapor

specific volumes and the equilibrium vapor pressure at a given temperature.

The present algorithm is developed for these data although the overall

methodology is general and can be modified to accomodate other data.

55
Table 6.1. Summary of Thermodynamic Property and Associated Subroutines

Name De sc r ipt ion/ Input s/ Outjguts

Data Fitting Routine s

FITAB

outputs: a,b

Other Routines Called

determine pure component a,b parameters

inputs: T, p satf6 , V l>e , Vye , <o p

, « £f »y

FITF determine interaction parameters from

mixture data

inputs: T, x. , some combination of:

p sat, e' V £, e' Vv, e' x v,e

output : f 22

Routine to Access Stored Parameters

BCONST access eqn state parameters from data base

and calculate reference states for H,S

inputs: code nos. for pure comp. , fQ., f,

outputs: (contained in common blocks)

note: must be referenced once for each

mixture before calling any of

following routines.

Property Routines

PL I MIT, VIT

BCONST, BUBLT, ENTROP

ESPAR, HCVCP, PLIMIT,

VIT, ZXLSF

BDESC (via common

blocks), BUBLT, ENTROP,

ESPAR, HCVCP, PLIMIT,

VIT

BUBLT calculate bubble or dew point pressure ESPAR, PLIMIT, VIT

inputs: T, x of parent phase

outputs: P, x of incipient phase, V« , V

ENTROP compute specific entropy

inputs: T, V, x output: S (returned as function value)

HCVCP compute enthalpy and/or heat capacity

inputs: T, V, x, IQ - output qualifier

outputs: H and/or C^ and/or C (as

specified by IQ)

Auxiliary Routines (transparent to user but required ic executable element):

ESPAR

ESPAR

BDESC

ESPAR

PLIMIT

VIT

ZXLSF

Block data routine containing pure component data

compute a, b, C , f^ as functions of T, x from stored coefficients

find upper and lower bounds on pressure given T

specific volume for liquid or vapor as function of T, p

or equivalent one-dimensional minimization routine from

math/ statistics library (required only with FITF)

56
Pressure and liquid and vapor volume are more than sufficient to determine 'a'

and 'b' and thus the parameters are chosen to minimize the following sum of

squares:

ru.,x> . ^(^f+ »v(^-)

2

* »pf-^7-f «.i>

where the e and c subscripts refer to experimental data and calculated value s

respectively and the ta terms are weighting factors which can be adjusted to

match the reliability of the various measured quantities. The calculated

pressure and volumes are evaluated using the equation of state subject to the

constraint that the Gibbs free energy of the liquid and vapor phases be equal.

The function T is at a minimum where the partial derivatives with respect to

'a' and 'b' are zero:

*£ = o = J"*,' - \.°yjl>£ + . /Vv,.-Vy,c\!Vv^

a A v

2

J

da ~

y

\ y

2 / -C» e v, e

+ w /p sat,e Psat,c\ ?P_s_at^c (6#2)

P

\ «2 / 3a

N Psat.e '

*L m o . JH.* ' Vl.c) tlltc + w /Vv.e-Vv^c^jlv^c

3b \ v\ e I ab v

\ v* ) *b

+ u feat.e Psat,c\ a Psat,c (6#3) P\ 2 / 3b

X Psat,e '

Because the equation of state is fifth order in volume, Equations 6.3 and 6.4

cannot be solved explicitly. To set up an iteration in a and b the calculated

57
quantities resulting from small changes in a and b (6a and 8b) are expanded in

terms of partial derivatives; for example:

3V/ c (a,b,T) dVp c (a,b,T)

V£,c (a + 8a' b + 6b ' T) = v&c(a * b ' T) + aa 8a +

Qh

8b (6 ' 4)

Hie combination of equations 6.2 and 6.3 with equation 6.4 and similar

expressions for V„ „ and P 00+. „ yields:

(i) , !L-(*?l*o\

2 + _^L_ /^M + *P /

8p sat t c\

2

l 6a + v

2 \ da / y

2 V 6a J 2 \ da / J £, e v, e F sat, c J

Lvi

i£_ dYLc dYl.c + _^v_ aVv,c avv,c + S

da db v

2 da db

e v, e sa p_ a Psat.c a Psat.c

da db

t, e J

6b

da

Y^c +

^Psat.e ~ Psat.c ) a Psat,c

'P\ 2 /

' da

K sat, e

(6.5)

±L- dV^c 6V£,c + _^y_ aVv,c aVv,c

u2 db da ,,2 db da

L v £,e v, e % a Psat.c

,2 3b

'sat, e J

8a +

v

2 \ db / v

2 v db y L x., e

¥ v,e

fa>p /3P

2 \ db

p sa t,

sat, c 6b

(Continued on next page)

58
= J*L,'-*l,c) m^ + (*,,. - \,<)

"^ , \ v 2 /3b v

X \ v 2 /3b V£,e X Vv,e

/Psat.e Psat.c] a Psat,c

(6 6) P V p 2 / 9a

F sat, e

where all quantities are evaluated at (a,b). The partial derivatives are

approximated numerically, for example:

*V/.c ~ Vfc<* + Aa ' b > - Vfc<*.*>>

(6 7)

da Aa

3V/ c „ V/(a,b + Ab) - V/(a,b) —^5 = JL L ( 6 . 8 ) 3b Ab

Equations 6.5 and 6.6 form a linear system in 6a and 6b. Given starting

values of a and b the system can be solved to give improved guesses:

,(1+1) « a (D + 6a (6.9)

b

(i+l) = b (D + 6b (6.10)

where the subscript is the iteration index.

In the implementation of the above method, the routine FITAB generates

starting guesses of a and b using the saturated volume data alone. Starting

at b = 1/2 V» the value of the a parameter is expressed in terms of b by

equating the pressure of the liquid and vapor:

59
RT(1 + yi + ji

2 - Ji

3

) RT(1 + yv

+ yv

2 - yv3

)

a= (6.11)

VV£

+ b) Vv

(Vv

+ b)

where

b . b

Yl " 4V ; 7yr ~ 4V

I, e v, e

A Newton's method iteration is carried out to find the value of b which

satisfies

Y(b) = = G£(a,b,T,V^ e

) - Gv(a,b,T,Vv>e ) (6.12)

If only specific volume data were available this iteration could be used to

determine a and b.

With the above starting values for a and b the iteration outlined in eqns.

6.5 to 6.10 is then carried out by FITAB. For each guess of (a,b) an inner

iterative loop is used to calculate the pressure and liquid and vapor volume.

This iteration is based on the equality of the liquid and vapor Gibbs free

energy and is identical to the one used in the BUBLT routine discussed below.

The iteration is repeated at (a + Aa,b) and (a,b + Ab) in order to evaluate

the partial derivatives (where Aa = 0.001 a and Ab = 0.001 b). The system of

equations 6.5 and 6.6 is then evaluated and solved for 5a and 6b to generate a new (a,b) by eqns. 6.9, 6.10. The iteration is repeated until 'a' and 'b' are

determined to a specified accuracy.

60
The a and b determined for individual data sets are then fit to simple

functions of temperature (such as eqns. 6.14, 6.15) to arrive at the final

expressions required by the property routines. Users are warned that because

the FITAB routine assumes saturation data it cannot be used at or above the

critical point. Furthermore, attempts to force the equation of state to fit

liquid-vapor equilibrium data near the critical point (i.e., at a reduced

temperature above about 0.95) will cause a degradation in accuracy for points

near the critical temperature but removed from the critical region. This

effect occurs with any equation of state that arises from the so-called mean

field approximation. This class includes all industrial equations of state.

The subroutine FITAB computes the a and b parameters at a single temperature

from saturation data at that temperature. If the various saturation

quantities are not all measured at the same temperature or if superheated

vapor data is to be used in computing a and b, the subroutine FITAB cannot be

used, rather it is necessary to carry out a non- linear least squares

regression of the entire experimental data set. Non-linear regression

routines are commonly available in statistics libraries but because they, as

well as the data to be fit, vary widely in form only a general discussion of

using such a routine is given here. A regression routine would find the set

of parameters that minimizes the sum of squares of residuals between the

experimental data points and corresponding calculated values. Typically, an

external subroutine computes the residuals as a function of the parameters to

be optimized (e.g., aQ, a^, aj , bQ, b^, b2 in eqns. 6.14 and 6.15). The process

is iterative and generally requires reasonably good initial guesses for the

parameters. For each set of parameter estimates generated by the iteration in

the regression program the external model subroutine is called to compute a

61
revised set of residuals. This external subroutine must therefore update the

common block containing the a. and b i coefficients and call the routine BCONST

(which computes reference state values as described below) before computing

the residuals with BEBLT or similar routines.

Determination of Int er action Parameter

The mixing rule for the 'a' parameter for a mixture involves the interaction

parameter, fij, which must be determined from experimental data. The approach

is again to find the value of fj» which minimizes the sum of squares of

relative deviations between measured and calculated quantities:

+ "x < xv,e" Vc>

2 (6 - 13)

Because the vapor composition must be between zero and unity, the last term in

eqn. 6.13 is expressed as an absolute error. Only a single measured quantity

in addition to T and x» is necessary to determine fj^* If all of the data

indicated in eqn. 6.13 is not available, the corresponding weighting factors

in the expression for T are set to zero.

The subroutine FITF will accept any combination of Psat> \ > V

y

and x^ data

but requires the a and b parameters for the pure components (i.e., pure

component data in the BDESC routine described below). For given values of T

and %n and the current iterative value of f 12 ' tne BUBLT subroutine (discussed

below) is used to provide the calculated pressure, volumes and vapor

62
composition needed to compute eqn. 6.13. The actual minimization is carried

out by the subroutine ZXLSF contained in the proprietary IMSL library accessed

through the MBS central computer. It is expected that users will have access

to a similar one-dimensional minimization routine and can modify the

appropriate lines of FITF.

The values of f,

2 calculated for different individual data sets should be

constant, or at most a function of temperature, for a given pair of pure

components. A wide variation in f .« with composition indicates either poor

mixture data or an inability of the equation of state to represent the entire

composition range of the mixture. In the latter case, the user must choose a

value of f|« consistent with the composition range of interest. The values of

f]9 determined near the pure component compositions typically will have greater uncertainty than values determined in the middle of the composition range.

Calculation of Equation of State Parameters from Stored Coefficients

o The coefficients of the curve fits to the parameters a, b and C for 11 pure

refrigerants are stored in common blocks initialized in the BLOCK DATA element

BDESC This element also contains mol ecul ar weights, critical properties,

reference state conditions for enthalpy and entropy and character variable

representations of the refrigerant names. The refrigerants currently in the

data base are listed in Table 6.2. The data arrays are dimensioned for 20

components so that users can add data for other pure materials.

The information in BDESC for the particular refrigerant pair of interest is

accessed by the subroutine BCONST and stored in subsidiary common blocks which

are referenced as needed by other subroutines. The enthalpy and entropy

63
Table 6.2. Refrigerants Currently Included in Data Base

Code No. ASHRAE Designation

1 Rll

2 R12

3 P13

4 R13B1

5 R14

6 R22

7 R23

8 R113

9 R114

10 R142b

11 R152a

Chemical Formula

trichlorof luorome thane

di chl orodif1 uor ome thane

chlorotr if luorome thane

br omo t r i f1 uor ome tha ne

tetraf luorome thane

chl or odif1 uor ome thane

tr if 1 uor ome thane

1,

1

,2- tri chl oro tr if1 ur oe thane

1 , 2-di chl or ote tr af 1 uor oe thane

1-chloro-l, 1-dif1 uor oe thane

1, 1-dif luoroe thane

64
routines are called at the reference temperatures of each component to

establish reference values of E and S. The inputs to BCONST consist of the

code numbers for the pure components (as listed in Table 6.2) and the interaction coefficient for the mixture (expressed in terms of f^ and f, as given in

eqn. 6.17). For a pure component the code numbers will be identical and f«

and fi should be specified as zero. It is necessary to call BCONST only once

(before any other property routines are referenced) for each pair of pure

components considered.

The thermodynamic property routines reference the subroutine ESPAR when it is

necessary to calculate the a and b parameters as a function of temperature and

composition. When called by the enthalpy, entropy or heat capacity routines

o the C for the pure components and the temperature derivatives of a and b are

also calculated. This arrangement completely isolates in ESPAR the

temperature and composition dependence of the equation of state parameters

except for the Gibbs free energy and chemical potential statements in BUBLT.

Thus, alternative temperature dependencies or mixing rules can be accommodated

by changing only two subroutines (plus of course the corresponding data in

BDESC) .

As presented in Appendix B the ESPAR routine calculates the a and b

parameters as:

a = aQ exp(a

x T + ajT2

) (6.14)

b = b + bjT + b

2 T

2 (6.15)

65
The ideal gas heat capacity is represented as:

C* = cQ + c

x T + CjT2 (6.16)

The interaction parameter is expressed as:

f 12 " f + f

l T (6 ' 17 >

For mixtures with an interaction parameter independent of temperature, fi is

set to z ero.

Values for the a

A

, b^ and c^ are tabulated in Appendix A for 11 refrigerants.

The value of the interaction parameter, fio» f°r several mixtures for which

p-V-T-x data are available are also given in Appendix A.

Dew and Bubble Point Pressures

The hard sphere equation of state represents the p-V-T behav ior as an explicit

expression for pressure in terms of specific volume and temperature. The

application of the equation of state for a pure component, however, often

requires saturation pressure as a function of temperature alone. For a mixture the bubble or dew point pressure is a function of temperature and the

composition of the appropriate phase as well. This calculation is carried out

by the subroutine BUBLT: given the temperature of a mixture (or pure

component) and the composition of one phase this routine calculates the saturation pressure, the composition of the other phase, and the specific volume

of each phase.

66
Vapor-liquid equilibria becomes meaningless, of course, above the critical

temperature and BUBLT returns a warning message without carrying out any

calculations when the temperature approaches the critical temperature. For

pure components, a warning is returned at temperatures greater than 0.99 of

the critical temperature predicted by eqn, 5.7. For mixtures the limit is

0.99 of the critical temperature of a hypothetical pure material having the

same a(T) and b(T) as the mixture. The critical temperature of a mixture is

always higher than that of the corresponding pseudo-pure material and thus,

BUBLT will return a message of 'critical point exceeded' for temperatures

which may be significantly below the mixture critical temperature. This limit

does not arise from the equation of state itself, but results from a failure

of the solution logic above the pseudo-pure material critical point. An

alternate version of BUBLT which allows computation to the mixture critical

point is under development.

All physically meaningful solutions to the equation of state must have a

negative slope of pressure with respect to volume. Furthermore, for there to

be both liquid and vapor solutions for a pure component the pressure will lie

between the limits where (6p/dV) T

= 0. These limits are a function of

temperature and are indicated as Pi ow and p in Figure 6.1. The slope of

6p/3V was given by eqn. 5.2 and is repeated here:

(*J>) = RT /y* _ bV3 _ b¥ + bjv __bl\+ a(2V + b)

V9V/ / . x 4 \ 4 16 256/ v2, v + ^

ylL_ b\4

V

4 16 256/ V

2 (V + b) 2

(6.18)

In the limit as V approaches b/4 (i.e. , as the volume approaches the molecular

volume) Op/8V)j approaches negative infinity. At a volume of V = 3.007b the

67
slope is positive for all temperatures below the critical temperature at least

for values of a and b which are reasonable for refrigerants. (The constant

3.007 is the ratio of the critical volume to b(T J.) At the other extreme, as

c V approaches infinity Op/dV)TM approaches zero but remains negative. Thus

between the volumes of b/4 and infinity the slope dp/dV must twice pass

through zero, corresponding to Pi ow and ?„„• la the critical region of a

mixture the pressure is not bounded by Vi cm and p and the solution logic

fails; this is the reason that BUBLT is limited by the critical temperature of

the corresponding pseudo-pure material.

The pressure limits Pi^ and p are found for a given a,b and T by the

auxiliary subroutine PLIMIT. Starting at V = 3.007b the volume is decreased

and eqn. 6.18 evaluated until a negative slope is found. Using these limits,

the bisection technique is applied to volume to find the point at which

Op/dV)y = 0. The value of pj is then evaluated from the volume by eqn.

2.6. If pi is negative, a small positive value is used. To find the upper

bound on pressure, the volume is increased from 3.007b until a negative slope

is found*, the bisection method is then used once more.

For a pure component, the equil ibrium criteria is the equal ity of Gibbs free

energy between the phases:

G£(T,V£) = G

y

(T,V v

) (6.19)

The specific volumes Vp and V are in turn functions of T and p. Graphically,

eqn. 6.19 is equivalent to the requirement that the areas between a line of

constant pressure and the p vs V line be equal as indicated by the shaded

68
areas in Figure 6.1. The general method of solution is to iterate on pressure

starting at the limits of p^ and p . For each guess of pressure, the

liquid and vapor volumes are found as a function of T and p using the

subroutine VIT (discussed below). The G for each phase is then computed as a

function of T and V by eqn, 4.8 (except that the G terms will always cancel

and are thus omitted). The iteration is continued by a combination of the

secant and reguli-falsi methods until the computed liquid and vapor Gibbs free

energies are equal to within a convergence tolerance.

The computation of saturated volumes as a function of pressure requires an

iterative calculation which is contained in the subroutine VIT. Newton's

method is used to iterate on volume until the pressure computed by the

equation of state agrees with the input pressure. The previous converged

value of Vp or V is used as the initial guess. To speed convergence the

iteration is carried out in terms of ln(V)*, for the vapor volume iteration the

pressure is transformed into logarithmic coordinates as wel 1. Convergence is

also aided by constraining guesses of vapor volume to values greater than the

critical volume*, liquid volumes are constrained to lie between b/4 and V_.

The need for finding the upper and lower bounds on the pressure iteration

stems from a potential problem with solving for V(T, p). For example, if a

pressure greater than p is input to the iteration for vapor volume no

solution exists on the vapor branch of the p-V curve and the iteration wil 1

either fail to converge or will converge to a solution on the liquid branch.

This latter case would result in apparently equal liquid and vapor volumes and

consequently equal Gibbs free energies for the two phases. This would be a

trivial and physically meaningless solution to eqn. 6.19. This is a problem

69
UJ

DC

CO

CO

LU

DC

Q.

—

1 r

Critical point Liquid-vapor

coexistence curve

Liquid branch

Vapor branch

T=TC

"

Vliquid

V

vapor

MOLAR VOLUME

Figure 6.1. p-V-T behavior of a pure material as represented by the equation

of state; vapor-liquid equilibria criteria represented by

equality of shaded areas.

70
of the solution logic and not the equation of state itself. By bounding the

pressure iteration a false solution is avoided.

For a mixture the criteria for liquid-vapor equilibrium (in addition to

temperature ard pressure equality) are the equality of the chemical potential

of each component in the two phases:

|i al (T,x1 ,V1

) = ua2 (T,x2 ,V2

) (6.20)

^

lbl (T' xl' Vl

) = »1b2 (T' x2' V2

) (6.21)

where the subscripts refer to components a and b in phases 1 and 2. The

temperature of the mixture and the composition of the parent phase are input

quantities. (The phases with the known and unknown compositions are referred

to as the parent and incipient phases respectively. For a bubble point

calculation the parent phase is liquid* a parent phase of vapor gives the dew

point pressure.)

Two concentric iteration loops are employed in the BUBLT routine to solve for

the saturation pressure and the composition of the incipient phase. For each

guess of pressure the volume and chemical potentials of the parent phase are

computed. The inner iteration loop for composition is entered and the

incipient phase a, b, V and u's are computed for the current guess of x2

.

The chemical potentials are combined with the current guess of x

2 to yield the

function used with secant/regul i-fal si iteration:

¥(x,) = —^

x, (6.23) 1 Z„ + Zr L

71
where z

fl = Xj exp (|ji al - |i a2 ) (6.24)

zb = (1 - x

x

) exp (nbl - jib2 ) (6.25)

For each guess of pressure the composition iteration continues until ¥(x9

) is

within a convergence tolerance of zero. Successive guesses of pressure in the

outer iterative loop are generated by the secant/regul i-fal si method and the

function:

<Mp) = 1 - (z

a + zb

) (6.26)

This function goes to zero when the chemical potentials are equal between the

phases.

The routine BUBLT applies the pure component iteration for Gibbs free energy

to the mixture to calculate the saturation pressure of a pure component having

the same properties as the parent phase of the mixture. This mixture is

unstable with respect to two phases at this pseudo-pure component pressure and

thus it represents a lower bound on pressure for a bubble point calculation

and an upper bound for a dew point. The pseudo-pure component pressure is

used as the starting value in the mixture iteration. The second guess for

pressure (necessary to start the secant method iteration) is given by:

P

= < z

a + zb)P (for tub^le point)

(6.27)

p

(2) = p

< 1 )/( z + z ) (for dew point)

To start the composition iteration an initial guess of equal parent and

incipient phase compositions is made. For the second and subsequent guesses

72
of pressure, the Xj iteration is started with the previous converged value of

Xa* In each case the second guess value for composition is given by:

A2) m _ia (6.28) 2 z a + 2b

A flow chart of BUBLT is given as Figure 6.2. This routine imposes numerous

checks and bounds on the calculations. Although these slow the execution time

of the routine they were found to be necessary to insure reliable convergence,

particularly near the critical region. BUBLT has been checked for numerous

refrigerant mixtures and was found to converge quickly and reliably for a wide

range of temperatures. Occasional convergence problems were encountered at

very low reduced temperatures (< 0.2 T ) but this is generally well below the

v freezing point where the equation of state is no longer valid anyway.

Enthalpy, Heat Capacity and Entropy

The calculation of enthalpy, heat capacity and entropy are straightforward,

non-iterative implementations of the expressions given in Section 4. Molar

specific enthalpy and heat capacities at constant volume and pressure are

computed by the subroutine HCVCP. Depending on the value of an input

parameter one, two or all three of these quantities are determined as a

function of temperature, composition and specific volume. The function

subroutine ENTROP computes molar enthalpy by eqn. 4.49, also for a given T, x

and V. The calculations for E and S employ the reference state calculated by

BCONST for each pure component.

73
cInput T, Xi )

i

a 1, b^fCT, *1)

Print

warning

C Return J

New

P

Compute:

v^, v v

, G^, G v

Return

Initial guess:

P = Ppseudo-pure

1

Compute:

v-|> Mai- Mb!

Initial guess:

*2 = *1

a2, b2 = f(T, xg)

I

Compute:

V2. Ma2. Pb2

converged

New

x 2

^^

us \

converged

No New

P

|Yes

Return )

Figure 6.2. Flow chart for the calculation of dew or bubble point

pressure

74
Iteration Technique

Several iterative calculations required to compute thermodynamic quantities

employ a combination of the secant and reguli-falsi methods to converge to a

solution. The equations to be solved are arranged to give a function ¥(0)

which is zero at the converged value of the iteration variable, 6. In the

secant method an improved guess for is generated by assuming a linear

function ¥(0) between two successive guesses and solving for ¥(©) = 0:

6

<i+ » =0<i>

- 4<(0 (i)

) f.f

- 6(A

T" (6.29)

¥(0 (l)) - ¥(0 (l 1})

where the superscript is the iteration index. This is represented in Figure

6.3 as the extrapolation from points (2) and (3) to find o' 4

. (If the

derivative 8^/30 can be expressed analytically the bracketed term in Eqn.

6.29 is replaced by 1/ (9^/80) to generate Newton's method.) The most recent

values of and ^(0) are retained and the iteration continues. The regulifalsi method is similar except that the new guess for is based on points

having a positive and negative value of ¥ :

e

<i+i> = e

<p°s) _ Y ( e<Pos >) e

(pos) - a

(neg)

(6.30)

¥(6 (pos)

) _ ^(neg))

This method is shown in Figure 63 as the interpolation between points (1) and

(3); note that the oldest guess for is not necessarily discarded.

Given starting values for which result in positive and negative values of

f(0) (such as the pressure iteration starting at p, and p„_) the reguli-falsi J. OW ^]? method bounds the solution and will always converge for a continuous function ^(0)

75
Function +(0)

<X>

neg

Secant method

iteration

V^pos Reguli-falsi

/ N method iteration

neg

Solution

ITERATION VARIABLE 6

Figure 6.3. Graphical representation of the secant and reguli-falsi

iterations to solve for ^(6) = 0.

76
The secant method, because it is always based on the two most recent guesses

of 6, often converges faster but can also fail to converge or even diverge

from the solution, In combining these two methods, a new guess for 6 is

generated by the secant method. However, the information necessary for the

reguli-falsi method is retained so that if &' 1+ ' given by the secant method

lies outside the bounds set by reguli-falsi (i.e., ©'Pos ^ and 0^ ne 6') a new

guess for Q^ 1

' is generated. For example in Figure 6.3, the secant method

extrapolation from (2) and (3) to (r ' is diverging from the solution in this

case the reguli-falsi interpolation between (1) and (3) is used to compute a

(A l\ new guess, &

y '.

77
7. EXAMPLES OF THE APPLICATION OF THE CSD EQUATION OF STATE

In this section, the perturbed hard-sphere equation of state developed in the

preceeding sections is applied to a particular class of fluids—the

halogenated hydrocarbon refrigerants—and their mixtures. First, the

thermodynamic characteristics of 11 pure refrigerants as represented by the

equation of state are compared with published values. For most of these

fluids, tabular values taken from ASHRAE (1981) are used, however, a detailed

discussion for R113 using the original literature data is also included.

Finally, the ability of the equation of state to represent the behavior of

mixtures is considered for two cases: the nonazeotropic mixture E13E1/R152a

and the system R22/R12 which has an azeotrope.

Pure Refrigerants

The saturation data of ASHRAE (1981)* were used to generate the a and b

parameters of the equation of state for 11 common one- and two-carbon

halogenated hydrocarbon refrigerants. The saturation pressure and saturated

liquid and vapor volumes over a reduced temperature range of approximately 0.6

to 0.9 were supplied (with equal weighting) to the programs described in

Section 6 to obtain the a and b which best fit the data at a particular

temperature. These individual values were then fit to the following functions

of temperature:

a = aQexpUjT + a2 T

2

) (7.1)

b = b

Q + bjT + b2 T

2 (7.2)

The data of the 'SI Unit Formulations' were used; these data are improved and

revised from earlier ASHRAE compilations presented in English Units.

78
These temperature dependencies are equivalent to those suggested by DeSantis,

et al. (1976) except for the addition of the Tr terms. The perfect gas heat

o capacities, C . were obtained from literature values derived from

spectroscopic measurements and were fit over a limited temperature range to

the function:

C* = cQ + CjT + CjT2 (7.3)

The values of the a^, hi and c- for each refrigerant are given in Appendix A.

The success of the equation of state in representing the ASHRAE data is

indicated in Table 7.1. The differences between the tabulated values and

those predicted by the equation of state (using the a and b given by Eqn. 7.1

and 7.2) are listed as RMS differences over the temperature range indicated.

The RMS differences for the saturation pressure ranged from 0.10 to 1.34

percent with an average for all 11 refrigerants of 0.54 percent. The liquid

volume representations ranged from 0.01 percent to 0.34 percent with an

average of 0.09 percent. The RMS error for the vapor volume ranged from 0.14

to 1.16 percent with an average of 0.50 percent. These differences are on the

same order as the accuracy to which the corresponding experimental quantities

are known (for example, see Mears et al. (1955)).

The temperature dependence of the differences in saturation pressures and

volumes for R12 and R22 are shown in figures 7.1 and 7.2. The smooth curves

are a result of comparing the equation of state with tabular data generated by

another correlation rather than with experimental measurements. The greater

differences at higher temperatures are primarily a result of the approaching

79
Table 7-1. Comparison of ASHRAE and Equation of State Values of Saturation

Pressure, Saturated Liquid and Vapor Volume, and Enthalpy of

Vaporization for 11 Pure Refrigerants

Temp. Range of RMS Difference (%)*

Refrigerant Cor relation (K) Psat V*

Vv ABva

Rll 240-420 0.41 0.04 0.45 0.87

R12 200-340 0.26 0.05 0.29 0.74

R13 180-272 0.10 0.02 0.14 0.17

KL3B1 200-300 0.64 0.07 0.57 1.46

R14 140-200 0.18 0.01 0.18 0.35

R22 220-330 0.50 0.06 0.45 1.13

R23 180-272 1.34 0.17 1.16 2.88

R113 240-43 0.22 0.03 0.23 0.61

R114 210-370 0.72 0.34 0.70 2.33

R142b 220-360 0.92 0.10 0.83 2.41

R152a 200-340 0.62 0.09 0.53 1.88

Average 0.54 0.09 0.50 1.35

RMS Difference = 1N

~ N / Y (ASHRAE - Aj V EOS

EOS

r 100%)

1/2

80
LUO

z

LU

C£LU

LUO

CCLU

Q.

OO

CO

O

LU

I LU<

DCI

CO

<

CO

o

LU

1.0

, r

,

0.5 Vapor volume

\^^~~~^\ /

^^ ""* L—Liquid volume '

-0.5

^-Saturation pressure

-1.0 i i 200 250 300

TEMPERATURE (K)

350

Figure 7.1 Percent difference between ASI7P.AE values of saturation pressure

and liquid and vapor volumes for EL2 and values correlated by

equation of state.

LU CO

LL o Li. LU Q 1 LU h- <

Z CC

LU I O CO

cc < LU

0.

0.5

-

-0.5

-1.0

1.5

Liquid volume

Saturation pressure

Vapor volume

i 200 250 300

TEMPERATURE (K)

350

Figure 7.2 Percent difference between ASHRAE values of saturation pressure

and liquid and vapor volumes for R22 and values correlated by

equation of state.

81
critical point; there is also an increased error in Eqn. 7.1 and 7.2 in

representing a and b at the extremes of temperature used in the data

correlations. The equation of state consistently overpredicts the pressure

and volumes for R22 (fig. 7.2). Although one might expect the differences to

be centered about zero, the observed behavior is the result of simultaneously

fitting three quantities with two parameters which interact in the equation of

state.

A more stringent test for the equation of state is its ability to predict

quantities not used in the correlation of the parameters. The enthalpy of

o vaporization is such a quantity (and one that is independent of C ) . The RMS

differences between the predicted and tabulated values (given in Table 7.1)

vary between 0.17 and 2.88 percent with an average of 1.35 percent. Although

these differences are larger than those for the pressure and volumes it is

important to note that in many cases the enthalpy values reported in ASHRAE

o are not measured cal or imetr ically but are derived from p-V-T and C

P

information.

The estimation of heat capacity involves a second derivative of the equation

of state and is thus a particularly stringent test of physical consistency.

The equation of state has been integrated over the volume to give the

Eelmholtz free energy function, from which all the other functions, including

heat capacity, arise. Associated with that integration is a temperaturedependent constant which can be evaluated by knowning the perfect gas heat

capacity of the fluid(s) of interest. As a test of the equation of state let

us compare values of C for the saturated liquid state of R152a from several

sources. The values calcualted by the equation of state use the perfect gas

82
properties presented by Chen et al. (1975). Also considered are the values of

C

D appearing in the ASHRAE (1981) tables, an experimentally-based data

correlation published by the National Engineering Laboratory (Cartwright,

1981) in the United Kingdom, and a single experimental value of C measured by

Radermacher (1983). These four sets of information are compared in figure

7.3. The four sets agree well at the high- temperature end of the range. At

lower temperatures, however, the ASHRAE values are consistently below the

experimental information of NEL and Radermacher; this is an indication of the

shortcomings of data correlation schemes that are in common use. The most

striking feature of figure 7.3 is that the values of C predicted by the

equation of state compare well with experimental information even though they

have been evaluated without reference to any previously measured liquid heat

capacities.

Because of the data smoothing and correlation necessary to generate complete

sets of properties such as the ASHRAE tables, the use of original literature

data to generate the equation of state coefficients would have been

preferable. The ASHRAE tables, while not always based on the most recent

refrigerant data, represent an extensive compilation in a uniform format and

were thus a convenient source for the correlation of coefficients. The

formulation for 1,1,2-trif luorotrichloroethane (R113) given in ASHRAE is that

of Mastroianni, et al. (1978). The consequences of basing the equation of

state correlation on tabular data will be examined for this fluid.

Mastroianni, et al. report the saturation pressure of 0.9999 pure R113 at 26

temperatures with an estimated accuracy of 0.1 percent for pressure and 0.02

to 0.05 °F (0.01 to 0.03K) for the temperature of the thermostated bath*,

83
0.14

0.12

o

£

Q.

o

0.10

0.08

0.06

240

T

Refrigerant 152a

Saturated liquid

1

—Equation of state

NEL tabulation • ASHRAE tabulation

Measurement by

Radermacher

1 1

260 280 300

TEMPERATURE (K)

320

Figure 7.3 A comparison of the molar heat capacity at constant pressure

for saturated liquid R152a given by several sources

84
saturated liquid density at 9 temperatures with an accuracy of 0.1 percent;

and the p-V-T behavior of superheated vapor along 7 isochors at 0.1 to 1.2

times the critical density (yielding 42 data points) with a volume uncertainty

of 0.01 percent and temperature and pressure uncertainties as stated above.

They also measured the critical temperature but obtained the critical pressure

and volume by extrapolation. Using these measurements they fit the saturation

pressure to a 6 parameter expression for -dn (p..+) as a function of T to

Stt t obtain an RMS error of 0.28 percent. The 9 liquid density measurements were

represented by a 5 parameter expression with an RMS error of 0.12 percent.

The p-V-T behavior of the superheated vapor was represented with an RMS error

of 0.20 percent using the Martin-Hou equation of state (involving 13

parameters plus the critical temperature). These correlations were combined

with a 5 parameter expression for the perfect gas heat capacity to generate a

complete set of thermodynamic tables and a pressure-enthalpy chart (which are

reproduced in ASHRAE)

.

We now compare the literature data with the values predicted by the hard

sphere equation of state using a and b parameters based on ASHRAE data. The

agreement for the saturation pressure and saturated liquid volume (quantities

which were used in the correlation of a and b) shown in figure 7.4 are

excellent for reduced temperatures below 0.9 (439 K) . The RMS error over the

temperature range of the correlation (240-430 K) is 0.038 percent for liquid

volumes and 0.29 percent for saturation pressures. These deviations are only

slightly greater than those given in Table 7.1; they are also nearly identical

to those of Mastroianni, et al . over the same temperature range. Thus the use

of tabular data for R113 does not appear to degrade the fit of the equation of

state for conditions which are within the limits of both the correlation which

85
1.0

0.5 - 6

LU

QZ

<

<

I-

<

Q

Z

LLJ LU

H -0.5

LU

CD

O

< -1.0

> -?.o

UJQ

I-

-4.0

Z

LU

O

-6.0

cc

LU

a. -8.0

10.0

-12.0

:&

0.5

REDUCED TEMPERATURE (T/Tc)

0.6 0.7 0.8 0.9 1.0

o o

o°oo° o

Scale change

Saturated liquid volume

o Saturation pressure

J_

220 260 300 340 380 420 460

TEMPERATURE (K)

500

Figure 7.4 Percent difference between values of saturation pressure and

saturated liquid volume correlated by equation of state and

experimental data of Mastroianni et al . (1978) for R113.

86
generated the tabular data and the correlation of the CSD equation of state

parameters.

The equation of state parameters were generated using saturation data from

240 to 430 K. As shown in figure 7.4 the predicted saturation pressure and

liquid volume show serious deviations (as great as 10 percent) outside these

bounds. This behavior is an indication of the inability of the CSD equation

of state to quantitatively reproduce the critical region rather than a

consequence of merely extrapolating beyond the limits of the correlation.

A comparison of experimental and predicted pressures as functions of

temperature and volume in the superheated vapor range is shown in figure 7.5.

Apart from the critical region, the equation of state exhibits commendable

accuracy upon extrapolation to higher temperatures and/or away from

saturation. Note that the behavior at the critical temperature is accurately

represented as long as the actual critical region is avoided. Even though the

parameters were fit using only saturation data below 43 CK the equation of

state predicts the pressure of superheated vapor to within 1.2 percent at

temperatures as high as 532 K for reduced volumes between 1.43 and 10.0. The

hard sphere equation of state represents the behavior of fluids in a

physically meaningful way and thus, with care, it can be extrapolated within

reasonable bounds in the absence of data.

The correlation discussed up to this point was based on the saturation data of

ASHRAE. We now discuss a correlation of the equation of state based directly

on the data of Mastroianni, et al . This data was used with a non-linear least

squares technique as discussed in Section 6 to simultaneously generate the

87
CO

b

ui

Q

z

<

<

Q

z

LU

111 <:

iLU

z

o

I-

<

>

LU

Q

Z

111 O

cc

LU

Q.

-12

0.85

REDUCED TEMPERATURE (T/Tc ) 0.90 0.95 1.00 1.05 1.10

2

-

-4 -

-6

-8

-10

Reduced

volume

o10.3

• 3.32

a 2.00

1.43

1.11

O1.00

?0.84

T

* A *

400 420 440 460 480 500 520 540

TEMPERATURE (K)

Figure 7.5 Percent difference be tween values of pressure in the superheated

vapor region as a function of temperature and reduced volume

(V/Vc

) correlated by the equation of state and experimental

data of Mastroianni, et al. (1978) for R113.

88
coefficients aQ

, a1? &

2

, b

Q

, bj, b

2 (the ASHRAE-based values were used as the

initial estimates to start the iteration). Saturation data at reduced

temperatures above 0.9 were excluded from the correlation as were superheated

vapor data with reduced volumes below 1.5. The saturated liquid volumes were

given a greater weight than the pressure data to reflect the greater accuracy

of this measurement and also to give approximately equal total weights to each

of the three data types. As indicated in Table 7.2, the individual a. and hi

changed by as much as 31 percent between the two correlations. The values of

the actual a and b parameters, however, changed only slightly as indicated in

Table 7.2 for two representative temperatures. With the correlation based on

the literature data the fit of the superheated vapor is improved at the

expense of a slightly degraded fit for the saturation data. The total sum of

squares of the residuals decreased by 15 percent for the correlation based on

the literature data.

This comparison of the CSD equation of state with the original literature data

for R113 demonstrates the power of this expression to represent the p-V-T

behavior of a pure fluid. With the exception of the critical region, it

represents saturated liquid volumes to better than 0.1 percent and saturation

pressures and the pressures of superheated vapor on the order of one-half of

one percent. More importantly the CSD equation of state achieves these

accuracies for both liquid and vapor with a single expression involving only

o six adjustable parameters plus the three coefficients needed to represent C .

This is in contrast to the three separate expressions for R113 presented by

Mastroianni et al . involving a total of 25 parameters.

89
Table 7.2 Comparison of Data Correlations for R113 Based on Saturation Data

of ASHRAE and Saturation Plus Superheated Vapor Data of Mastroianni

et al.

Equation of State Parameters

a

Q (kJ m

3 /kmol 2

) a

x (kJ m

3 /kmol2

K)

a (kJ m

3 /kmol 2 K

2

) Dq (nr/kmol)

b

1 (n^/kmol E)

b« (m3 /kmol K

2

a @ T = 300 K (kJ m

3 /kmol 2

) a @ T = 400 K (kJ m

3 /kmol2

) b @ T = 300 K (m3 /kmol)

b @ T = 400 K (m3 /kmol)

Correlation Based on Data of

ASHRAE

»-3

7332.59

-2.20396 x 10^r -7.2656 x 10 7

0.230713

-1.87956 x 10_4

-1.06114 x 10 7

3545.8

2703.4

0.16478

0.13855

Mastroianni

7489.98

-2.34255 x 10_,

-5.03625 x 10 7

0.233100

,-3

-2.03232 x 10

-0.82063 x 10

3544.8

2707.3

0.16474

0.13 868

-4

-7

RMS deviation between EOS

Prediction and Data of Mastroianni (%)

1 iq , sa t (T = 282

0.038

425 K, 5 points)

0.057

sat

(T = 238 - 436 K, 23 points) 0.38

10.0, 19 points)

0.85

405 532 K, Vr = 2.0 -

0.51

0.59

Sum of squares of residuals: 16.61 14.19

90
Equilibrium Properties of the Mixture R13B1/R152a

As an example of a mixture which does not have an azeotrope we examine the

R13B1/R152a system. There is little published information about this mixture

and thus measurements of the saturation properties were undertaken at NBS.

The experimental details and measurements are presented elsewhere (Morrison

and Neal, 1985) but are briefly described here before proceeding to a discussion of the correlation of the data with the equation of state.

Samples for this study were prepared by distilling measured quantities of each

of the components from a gas buret into a stainless steel thermocompressor at

liquid nitrogen temperature. The amount of each of the components was determined in two ways: first by using the temperature-pressure-volume measurements from the gas buret and the values of the second virial coefficient

predicted by the equation of state described earlier in this paper (Eqn.

2.2); second, by weighing the thermocompressor after each successive addition of the components (Morrison and Kincaid, 1983) . The second virial coefficient for R152a agreed within experimental uncertainty with the values

listed by Dymond and Smith (1969); those values were derived from the measurements of Mears et al. (1955). Second virial coefficients for E13B1 were not

available in the literature. The two amount determinations typically agreed

with one another to 0.1%. The measurements included runs on each of the pure

refrigerant materials and on mixtures that were roughly 25, 33, 50, 67, and 75

mole percent R152a. The mixtures were moved from the thermocompressor into

the sample cell, constructed from a drawn sapphire tube having a volume of

approximately 7.0 cm (Davis, 1983). The volume accessible to the sample

could be changed by raising and lowering the mercury level in the sapphire

tube.

91
The cell was kept in a water bath whose temperature was measured with a quartz

crystal thermometer calibrated to an accuracy of 0.001 K with an NBScalibrated 25 ohm platinum resistance thermometer*, the temperature was

controlled to within 0.0003 E of a constant value. The total volume accessible to the sample and the volumes of the individual phases were determine by

measuring the distances between the top of the cell and the liquid-vapor

meniscus or the liquid mercury meniscus. The volume of the cell was cal ibrated to within 0.001 cm with triply distilled mercury. Pressures were

measured with a differential gauge calibrated to 0.2 kPa with a dead weight

gauge.

Measurements of the liquid and vapor volumes and the pressure were made at

five equally spaced temperatures between 15°C and 55°C. Several sets of

measurements were made by progressively enlarging the volume accessible to the

sample. The results most immediately determined from these data are the

liquid molar volume and the pressure on the bubble line. Although no samples

of either phase were taken during the experiments, there is sufficient information in these data to locate the dew point curve (Knobler and Scott, 1982).

The bubble point pressures are shown as open circles in figure 7.6. The

filled circles represent the composition of the dew phase in equilibrium with

the 50% mixture calculated from the data. The molar volumes of the saturated

liquid are shown in figure 7.7.

The equation of state was generated by first fitting the properties of the

pure materials in the fashion refered to in Section 6. A value of 0.089 for

the mixing parameter was found to optimize the pressure correlation along

92
3.0 -

2.5

CO

Q. 2.0

2

111

DC3 1.5

CO

CO

LJJ DC

a. 1.0

0.5

0.2 0.4 0.6 0.8 1.0

COMPOSITION (mol frac R152a)

Figure 7.6 Saturation pressure for the mixture R13B1/R152a. The open

circles, 0, are experimental bubble points and the closed

circles, #, dewpoints. The solid and dashed curves are the

predicted bubble and dew curves respectively beginning at the

lowest pressure and temperature, 15.902°C and incresing in

temperature to 25.323°C, 35.323°C, 45.134°C, and 55.033°C.

93
0.2 0.4 0.6 0.8 1.0

COMPOSITION (mol frac R152a)

Figure 7.7 Saturated liquid volumes for the mixture R13B1/R152a. The open

circles are the measured values and the lines, the values

predicted by the equation of state. The temperatures are identical to those of Figure 7.6.

94
the babble curve. We elected to use only the pressure data and not include

the volume data to find the mixing parameter because of the higher precision

of the pressure measurement.

The pressures and liquid volumes given by the equation of state are shown in

figures 7.6 and 7.7 as solid and dotted curves for the bubble and dew lines

respectively. The pressures are correlated to +/- 5 kPa; the predicted

volumes for the mixture are consistently greater than the measured volumes.

One should note the dramatic narrowing of the two phase region on the R13B1

side of the phase diagram, especially at the highest temperature. This

behavior does not indicate the onset of an azeotrope, rather it is an

indication of the nearness of the R13B1 critical point at 67°C.

Azeotropes and the Equilibrium Properties of R12/R22

In the above discussion, we examined the simplest kind of behavior one can

expect from a mixture of two liquids, complete miscibility and the monotone

variation of properties from one component to the other. We now consider an

example of the next most complex behavior case, azeotropy. Our motivation is

twofold, first, to discuss the situations when one can expect an azeotrope,

and second, to examine a refrigerant mixture for which there is calorimetric

data.

Let us consider the situation when one would expect an azeotrope. An

azeotrope is inevitable when the saturation curves of the two components cross

the p-T projection. Such a point of apparent intersection is called a

'Bancroft point' (Rowlinson and Swinton, 1982b). The mixture R12/R152a, also

known as RSOO, shows just such a behavior. In mixtures where there is a

95
Bancroft point, the azeotrope can exist over the entire composition range from

one pure component to the other. In R12/R152a the locus of azeotropes is

intercepted by the solid phase and the critical locus (Pennington, 1952). A

Bancroft point is not required for an azeotrope to exist. Since any kind of

non-ideal behavior causes a deviation from the Raoult law prediction, the

closer the vapor pressures of the two components, the more likely there will

be an extremum in the T-x or p-x lines and, hence, an azeotrope. Closeness of

vapor pressures is often associated with closeness of critical temperatures,

since, with a few notable exceptions, the critical pressure of many materials

(including halogenated hydrocarbons) are nearly the same, about 4 MPa, and

saturation lines are roughly parallel in a p-T projection. One can thus

conclude that two materials with nearly the same critical temperature are

likely to show azeotropic behavior. A sharp demarcation cannot be drawn,

however. A certain critical temperature difference needn't guarantee the

presence or absence of an azeotrope; the situation depends upon the detailed

molecular character of the components in the mixture.

The mixture R12/R22, also known as RS01, has an azeotrope that emerges from

the pure R22 axis of the phase diagram at approximately 320 E and moves into

the mixture region as the temperature is lowered until, at 253 K, the

azeotropic composition is approximately 10 mol percent R12 (Spauschus, 1962).

The variation of the azeotropic composition with temperature is not unusual.

The amount by which it varies with the temperature will however, be different

from one azeotropic mixture system to another. The azeotrope in R12/R22 is a

positive azeotrope (Rowlinson and Swinton, 1983b); that is, the boiling curve

will have a maximum pressure when measured at constant temperature or,

conversely, a minimum temperature when measured at constant pressure. This is

96
the kind of azeotrope one would typically expect from refrigerant mixtures.

The opposite kind of behavior is usually encountered when the two components

have a strong attractive interaction, such as the hydrogen bonding encountered

in chloroform/acetone (Karr et al. , 1951).

The mixture R12/E22 shows azeotropic behavior over a range of temperatures and

compositions as shown in fig. 7.8. The azeotropic point is typically found by

locating an extremism in the bubble point curve; at such a extremum, the second

law of thermodynamics requires the liquid and vapor phases to have the same

composition and the dew and bubble lines to be tangent (Bett et al., 1975).

Because of the flatness of these lines around an azeotrope, the uncertainty in

the composition is typically large. The data of Eiseman (1957) , for example

determine the azeotrope only to within 5 mole percent as shown by the error

bar in figure 7.8. The azeotropic composition can be determined more

accurately by a distillation process (Pennington, 1952).

The equation of state parameters for pure R12 and R22 were determined by using

data from the ASHRAE Tables (1981) as described above. The mixing parameter

was determined by fitting the bubble point predicted by the equation of state

to the constant pressure boiling temperatures measured by Eiseman (1957) over

the full range of compositions. The values of the mixing parameter were

slightly composition dependent. Equally weighting all values over the entire

composition range yielded an average of 0.041 with a standard deviation of

0.009.

The first test of the equation of state is the comparison of its prediction of

the temperature dependence of the azeotropic composition. The solid curve in

97
~ 0.24

CM

co 0.20

o

E

O

CO

o

Q.

o

o

o

o

cc

0.16 -

0.12 -

0.08 -

H 0.04 - O

UJ

N

<

I i — Equation of state

Data of: — Spauschus

a Loffler

o Pennington

• Eiseman

Madinabeitia

-

a

D \>

Nd S

n

1 A

200 240 280

TEMPERATURE (K)

320

Figure 7.8 The variations of the azeotropic composition of the R12/R22

mixture as reported by several authors and predicted by the

equation of state.

98
figure 7.8 shows that this prediction falls within the experimental

measurements of the azeotrope. As we have noted previously, there is

considerable uncertainty in those points. The dashed line in the figure is

the locus suggested by Spaschus (1962) from measurements he made on a grid of

composition- temperature conditions. Our prediction, which arises from data at

a single pressure nearly coincides with Spauschus's experimental locus.

The equation of state contains enough information to evaluate enthalpies as

long as the temperature is fixed. (For differences in enthalpy due to

temperature changes, one also needs the perfect gas heat capacities.) Neil son

and White (1959) have measured the enthalpy change associated with the

complete evaporation of mixtures of R12/R22 over the entire composition range

at 222.0 K. Their data, plotted in figure 7.9 are integral quantities

because the composition of the liquid changed as the evaporation proceeded.

The prediction of a closely related quantity, the enthalpy change associated

with the complete vaporization of the liquid mixture at fixed composition to

its vapor at the same fixed composition, is shown by the solid curve in figure

7.9. The agreement between these two closely related quantities is on the

order of +/- 0.5%. The equation of state also predicts nearly quantitatively

the curvature associated with the composition. By using these evaporation

data, Neilson and White were also able to evaluate a quantity closely related

to the enthalpy of mixing along the saturation line, H

e

. The comparison of

their calculated values and the predicted values are shown in figure 7.10. The

equation of state overpredicts this quantity. One should note however the

small magnitude of this quantity. Furthermore, the experimental and equation

of state quantities are not exactly the same: the experimental value should

underpredict slightly. This comparison with calorimetric data is encouraging.

99
21.6

I- 21.2

<

N

O o

5t

S 20.8

> a.

O

l-

<

Lil

20.4 -

20.0

I I I I

- o Data of Neilson and White -

^ Equation of State

_ \ -

- o\ -

- >v o

A

^\° o o >**

1 1 1 1 0.2 0.4 0.6 0.8

COMPOSITION (mol frac R12)

1.0

Figure 7.9 The heat of complete vaporization at 222. 00K for the mixture

R12/R22 as measured by Neilson and White (1959) and predicted by

the equation of state.

100
0.8

o

E

J*

o

o

CO

<

LU

0.6 -

0.4 -

O 0.2

1 1 1 1 o Data of Neilson and White

Equation of State

/ o

/ °

/ o

o

/° °\

^111 1 0.2 0.4 0.6 0.8

COMPOSITION (mol frac R12)

1.0

Figure 7.10 The heat of mixing for the R12/E22 mixtures as reported hy

Neilson and White (1959) and predicted by the equation of state,

101
That such intricate details about this mixture are accessible through so

little information is a demonstration of the power and versatility of the

Carnahan-Starl ing-DeSantis equation of state.

102
8. CONCLUDING REB1ARKS

This Note presents the development of an equation of state based on a hardsphere reference fluid*, the expression used is a modification of the CarnahanStarl ing-DeSanti s (1969) hard-sphere fluid proposed by DeSantis, et al (1976).

The resulting expression is simple, based on a good physical model and has a

physical interpretation of its parameters. The complete set of thermodynamic

functions derived from this equation of state represent with thermodynamic

consistency the behavior of both liquid and vapor phases of pure fluids and

binary mixtures. The CSD equation of state has been demonstrated to represent

the p-V-T properties of the 11 pure refrigerants investigated very well, often

within the accuracy to which the experimental quantities are known. It is

likely to do equally well for other similar materials (e.g., simple hydrocarbons and hydrocarbon derivatives) and simple inorganic molecules.

The representation of mixture properties with a single interaction parameter

(in addition to pure component information) was demonstrated, detailed aspects

of mixture behavior were accurately predicted even though the correlation of

the interaction parameter was based on limited data. There is, however, a

wide range of behaviors possible in mixtures and this is an area that requires

continued study. The present mixing rules need to be applied to many more

fluid pairs. In addition, the modification of the Carnahan-Starl ing expression to account for a mixture of hard spheres of differing diameters (discussed in Section 3) should be investigated.

The CSD equation of state concisely represents properties with thermodynamic

consistency and permits a rational estimation or extrapolation of quantities

103
for which no measurements are available. It is not, however, intended to

replace experimental data or preclude the need for experimental data covering

the entire range of interest.

The representation of both liquid and vapor phases with a single, physically

well-founded equation of state is the desirable state of affairs. Therefore,

this, or similar, equations of state should be considered as the basis for the

preparation of property tabulations for inclusion in handbooks, etc. Such a

preparation must be based on a correlation of carefully evaluated experimental

data (either original or literature values). Most of the numerical

coefficients presented in this work are based on a correlation of tabulated

data because of project time limitations and in some cases the unavailability

of the original data. For this reason, calculated values based on these

coefficients are not meant to replace compilations such as ASHRAE (1981).

At present, the primary use of the equation of state is in studies where the

properties of a number of fluids or mixtures are required. A single set of

computerized property routines would be usable for different fluids simply by

changing the numerical coefficients. Because of the small number of coefficients required, it would be of particular use with fluids for which limited

data is available.

104
REFERENCES

American Society of Heating, Refrigerating and Air Conditioning Engineers,

Inc., ASHRAE Handbook of Fundamentals, New York, (1981).

J. A. Barker and D. Henderson, Mol. Phys. . 21 , 787 (1971).

K. E. Bett, J. S. Row 1 in son, and G. Saville, Thermodynamics for Chemical

Engineers , MIT Press, Cambridge, pp. 377-380 (1975).

N. F. Carnahan and K. E. Starling, J. Chem. Phys. , 5JL, 635 (1969).

G. Cartwright, NEL Report No. D(PP)SR5, East Kilbride, United Kingdom,

National Engineering Laboratory (1981).

S. S. Chen, S. A. Rodgers, J. Chao, R. C. Wilhoit and B. J. Zwolinski, J.

Phys. Chem. Ref. Data , 4, 441 (1975).

H. A. Davis, Rev. Sci. Instrmn. , 54 , 1412 (1983).

R. Desantis, F. Gironi, and L. Marrelli, Ind. Eng. Chem., Fundam. , 15 , 183

(1976) .

J. H. Dymond and E. E Smith, The Virial Coefficients of Gases, Oxford

University Press, Oxford (1969).

B. J. Eiseman, J. Am. Chem. Soc , 79 , 6087 (1957).

J. F. Ely and H. J. M. Hanley, NBS Technical Note 1039, Washington (1981).

J. R. Fox, Fluid Phase Equilibria , 14 , 45 (1983).

A. Fredenslund, Guehling, J., and Rasmussen, P., Vapor Liquid Equilibria Using

UNIFAC, Elseview, New York (1977).

Z. I. Geller, W. F. Czoikowski and A. W. Egorow, Kholodil' naia Tekhnika i Tekhologiia, Kiev , 15, 70 (1972).

K. C. Gubbins and C. B. Twu, Chem. Eng. Sci. , 33_, 863 and 879 (1978).

L. Haar, J. S. Gallagher, and G. S. Kell, NBS/ NRC Steam Tables, Hemisphere

Publishing Corp., Washington 1984.

G.H. Hackstein, Thesis, Technical University of Dresden (1976).

D. Henderson, Practial Calculations of the Equation of State of Fluids and

Fluid Mixtures Using Perturbation Theory and Related Theories, in

Equations of State in Engineering and Research, Advances in Chemistry, v. 128, American Chemical Society, Washington, pp. 1-30 (1970).

D. L. Henderson and P. J. Leonard, in Treatise on Physical Chemistry, Vol.

VIIIB, Chapter 7, Academic Press, New York (1971).

105
J. 0. Eirshfelder, C. F. Curtiss, and R. B. Bird, Molecular Theory of Gases

and Liquids , pp. 162-170, Wiley, New York, (1954).

JANAF Thermochemical Tables, 2nd ed., D. R. Stull and H. Prophet, project

directors, NBS Circular 37, Washington (1971).

A. E. Earr, H. E. Schiebel, W. M. Bowes and D. F. Othmer, Ind. Eng. Chem. , 43 , 961 (1951).

H. Kammerl ingh-Onnes, Comm. Phys. Chem. Lab. Leiden, Nos. 71 and 74

(1901).

C. M. Knobler, Volumetric Properties of Gaseous Mixtures, in Specialist

Periodical Reports , Chemical Thermodynamics, vol. 2, M.L. McGlashan, ed.,

The Chemical Society, London, pp. 199-237 (1978).

C. M. Knobler and R. L. Scott, J. Chem. Phys. , 73, 5390 (1980).

J. L. Lebowitz, Phys. Rev. . 133, 895 (1964).

H. J. Loffler, Kaltetechnik , 12, 256 (1960).

M. J. Mastroanni, R. F. Stahl, and P. N. Sheldon, J. Chem. Eng. Data . 23 , 113

(1978).

D. Madinabeitia, thesis. The Effect of Pressures and Temperature on Azeotropic

Composition in the R22 and R12 Binary System. University of Delaware

(1960).

A. Mansoori and T. W. Leland, Jr., Faraday Transactions II , 68 , 320 (1972).

J. G. Mayer, J. Chem. Phys. , 5, 67 (1937).

W. B. Mears, R. F. Stahl, S. R. Orteo, R. C. Shair, L. F. Kells, W. Thompson

and H. McCann, Ind. Eng. Chem. . £7_, 1449 (1955).

J. Mollerup and A. Fredenslund, J. Chem. Eng. Data , 21, 299 (1976).

G. Morrison and J. M. Kincaid, AIChE Journal , 96A, 471 (1983).

G. Morrison, The Importance of Including the Liquid Phase in Equations of

State for Non- Azeotropic Refrigerant Mixtures, ASHRAE Transactions, 91,

Part 1 (1985).

G. Morrison and T. B. Neal, (to be published (1986)).

E. F. Neil son and D. White, J. Phys. Chem. , 63_, 1363 (1959).

D. Y. Peng and D. B. Robinson, Ind. Eng. Chem., Fundam. . 15 , 59 (1976).

W. A. Pennington, Ind. Eng. Chem. . 4_4_, 23 97 (1952).

106
W. A. Pennington, World Refrig. , 8, 151 (1957).

D. R. Pesuit, Ind. Eng. Chem. . Fnndam. , 17 , 235 (1978).

K. S. Pitzer and R. F. Curl, J. Am. Che m. Soc , 77, 3427 (1955); 79, 2369,

(195 9).

F. H. Ree and W. G. Hoover, J. Chem. Phys. , 40, 93 9 (1964).

R. C. Reid, J. M. Prausnitz and T. K. Sherwood, Properties of Gases and

Liquids , 3rd ed. , McGraw-Hill Book Company, (1977).

H. Reiss, Adv. Chem. Phys. . 9, 1 (1965).

J. S. Rowlinson and F. L. Swinton, Liquids and Liquid Mixtures, London,

Butterworth Scientific (1982).

J. V. Sengers, R. S. Basu and J. M. H. Level t Sengers, Representative

Equations of State for the Thermodynamic and Transport Properties of

Fluids near the Gas-Liquid Critical Point, NASA Contract Report No. 3424 (1981).

G. Soave, Chem. Eng. Sci. . 27, 1197 (1972).

H. 0. Spauschus, ASHRAE Journal . 4, 49 (1962).

H. E. Stanley, Introduction to Phase Transitions and Critical Phenomena,

Chapter 3, Oxford University Press, Oxford (1971).

F. P. Stein and P. C. Proust, J. Chem. Eng. Data , 16 , 389 (1971).

L. Verlet and D. Levesque, Physic a , 36 , 254 (1967).

L. Verlet, Physica , 30, 95 (1964).

K. P. Wallis, P. Clancy, J. A. Zollweg, and W. B. Street, J. Chem.

Thermodynamics , 16 , 811 (1984).

H. W. Wooley, International J. Thermophysics , 4, 51 (1983).

107
o &3 +> co -«S

<Ho 2 V •*^ ea o o CS

©

©

© ts VO f> es

© o o a +j es *• f~ © © en r~ en es en r- vo •«

C3 0. •H

UL,

•>» en eN en en en cs ^r en en en i O 1 o i O 1 © 1 o o o 00 u i ft © © 1 © </> • ea * © 00 © >»• © 00 en © rH cs © 11 p. +J cs es 1-1 cs i-H CS rH cs cs •* CS es cs o s co

»-> a>

Q

u

H in

a

a

O

m

m

ft

ec) CS 43 CO H ^-v CS a +->

T3 <s^ *.i oa « cs; a <*H ai «>H » »

<a o > * » » ~ « » u. i -s_^ > • « > ft «s> > > a 4)a. J^ > "ft >r a £ >

<u ft ft

CIoas

c~ *4- iH 00 00 CS

© OS CS en

CS rH "Of ft ©\© VO VO en •« VO vo

u 43 t- vo 0\ 00 •fl- ft r- vo

© ft •A

so o\ o r~ ft irt •* ft

© rH o It

V

O *4- o • r- • es • CS • • • 3

Oh

H • • o • en • ft • 00 00 rH en1 ft1

i—i I I rH1 on1 rH 001 t~ rHI rH1

<+Ho

rH ft »*• 00 "* rH en es "* o oo

in i-( en o\ t- <s r-

© ft en ^ ft 00

+J 43 «o »-(

T •* o\ CS s

cs <<r cs

T •H T *4» * ** CS "* VO en vo 0\ © =H O t~ 00 00 00 OS i-H CS © <f\ 00 o\ 1

rH

CS «S r-© o S 1 1 I 1 1 1 1 1 l 1 I

SI u

o\ ft i—i

o\ i-H rH t-

© 00 VO rH

u w-» ft t •«*• o 00 en

© rH o cs o NO r- •h ON VO VO rH rH en o «•*> <*H O vo en 00 ON. © en en en vo vo vo

43 r* «"» eN en

© rH o en

© ^, •H

</> cs en rH rH

*>a o o o ©

©

©

© o ©©O <u •HO

•H •«* 00

© 00 r~ en CS en en CS CS

*H CS 0\ i-i en *** vo ft ft vo rH •* en

>H CO rH en •») 00 i—t

vo ft en ^, vo •V

<U v© rH t~ es o os en on

©

O cs r~ o v> vo ft ft 00 00 ON «o ON rH en

u I—1 o1 o1

i-iIi-H1 -H1 rH1 © o1 en 1 1 H< ON o 00 i-H ft VO VO ft c~ VO 00

rH vo en eS on 00

©

© ft <r> ft ru 03 vo CS 00 e- ON t~ CS cs m © r-

-H en •<r r- rH 0>i rH 00 00 «* nn

T 00

4= o es r- •* c~ 00 en vo en 00

m ft

r* rH H 1 CS1 en1 cs1 1 esI I es1 1 oI oI

•* CS en

©

© on en 00 ft

© r-

«o <H i-i T-t vo >o OS e»\ en© en o rH *4* 00 00 en * ft 0\ rH o ^, ea r- CS o\ es o\ rH CS 00 r- ov ft o\ ft es t~ en <o © •*• r~ ON es * en CS CS rH es CM r- os CS cs

M

ft •au

•H

IAcou« +

ftH«O

ea

«N +

•a

4>V

IHftM

-H *

-v ,0

XI

hHQZe

rH

(*H M

es

-HOS

en

rH

caen cs2 2 Oh

43 CO 'O- «S CS

rH <er <o a 5

rH

Oh

.a

091Z

eau

<uftB

si J3a

ea

43

e«

00 H-> 44

en

w s « w 43

43 *

108
Table A.2. Coefficients for Curve-fit of Pur e- Component Perfect-gas Heat

Ca pa c i ty

Refrigerant c c

l 104

c 2

Temp. Range of

Data Fit (K) Source

Rll 22.0418 0.260895 -2.45319 200-400 JANAF

R12 17.5387 0.248546 -2.16271 200-400 JANAF

R13 13.93 00 0.232181 -1.82929 200-400 JANAF

R13B1 19.9537 0.216394 -1.70241 200-400 JANAF

R14 11.0629 0.209740 -1.40992 200-400 JANAF

R22 17.0547 0.161633 -0.91256 200-400 JANAF

R23 20.4760 0.106183 -0.12189 200-400 JANAF

R113 76.2637 0.119641 0.71879 240-420 ASHRAE

R114 20.7005 0.464035 -4.175 89 240-470 ASHRAE

R142b 23.7611 0.231706 -1.06534 250-500 Mears, et al

R152a 22.2804 0.154009 -0.03067 200-400 Chen, et al.

NOTE: Cp(kJ/k mol K) = cQ + CjT + ^T2

where temperature T is expressed in Kelvin

109
es *-, ,—

*~^ •ft vo

g

V. \© o\ l>

t- r~ i-t o\ ^-^ 0\ •ft ^v 0\ N-^ i-H <S

i-H 0\ •H i-l

w t-

•^r i-H t- s_^

a <3N >^ 0\ o a ft

a •H ft •M •H ~s o a *—

a oo o « (A « m a •p u u •H § a © •H M o 3 h •H a

•M i-H

M 09 s i-HU iH o o •H *-> o © OS 4> CO s w CO ss a. H ©

ao •ft i-H •w © o © o vo ©\

© •M © o o H » 1 1 1 1 1 1 1 p< eS © o o © •ft

«ft oI-H

©

a o o o o • • o o o O oa u 00 ON en o en en en u >_> cs en t- ON f-t en t- « en cs «s N en en cs n » I 1 1 1 1 1 1 « CO on «s © •ft en en en

o<

a oo en o "ft •t •ft ON §<s <s «s <s IN ts <s »H Ho a

a,

•ft on

NO#

en

noen

H^ t

•m -* ~ H

H

<S «V> CvJ

w

w 0) oH

asa«u©

SOn

«HeM8U©§u£ao

en<W

r-l

.£> eS

•M *» «H O *

ft •>

ft Z

es

<«H

e» cS ft§ f

•w

O

CS o

<A

N

ft 09

•o

a

n -^ * ON m © © g a o a o * "H © O

60 cs Q

ts oi-H ©o ooo oo

•s«

ft o o o

ON t-l 0\ •ft •ft 00 « 00 ^f oo en en en 00

00 o o o o o o o «u o o o o o o o 4) ^

CS e •ft CS 3 a «s en <s •ft ts M <-H & 2

I-l OS

a

a

a

•H

n ***^ «*^ **^ •*<^ «»» "«» as en cs crj en CS es en

i-H OS

u

a

i-H 03

aa 2

•a i-iCS M +* •o a a ©

CS

a

M i-i W n

g> o ft # H©© • Jd c -w

(H

4>

a

-a •w

TJu

©

u

4-1 (-> 14 H

-H CS »H u

0)

U) —

§ © K

tM

!-HCS 60a

•H

© U) CO —

a)

M

n cs ©S

es

i-i S-S a a

•H •H • o

ft

CS M

a

H <4H

CS©

•w© <«

•a 9 o©

-H

©

u

CS J5

a

a

•M 3 © <*-t a J O

S,

H

a e © o •a

ft -»H M M

a

© CS o O -H Ji IH •W o •o •H

CS

H

W

© •w

a

a

OS §

<I-I u> o

•o

H

u

© o a B!

i-i

© 1 a

i-H 3 » CS •M *> CS © N

•a o «

»- ft

M

•s (M 1 » © •w es-w

1? «H a CS •M •a a a o CS s

•H

•M V)

I-l CSw

cd iH

a > iH < © o

ft

©

110
APPENDIX B : FORTRAN SUBROUTINES

This appendix contains listings of the FORTRAN subroutines developed to

implement the equation of state for binary refrigerant mixtures as well as a

sample run which uses the subroutines to compute a table of properties for an

example mixture.

Notes on the use of these routines:

1) Language. All of the routines are written in ASCII standard FORTRAN

(FORTRAN 77). Generic names are used for all calls to intrinsic

functions (e.g., EXP and LOG)

.

2) Argument lists. Inputs and Outputs for each of the subroutines are

described in comment statements at the beginning of each routine.

3) Units. The subroutines as presented here require all inputs and compute

all outputs in SI units with the following multipliers:

composition: mole fraction

temperature: K

volume: m°/kg mol (equivalent to £/g mol)

pressure: kPa

enthalpy: kJ/kg mol

entropy, heat capacity and gas constant: kJ/kg mol K

The routines will work with any other set of consistent units if the

values of the gas constant and equation of state coefficients are

converted in the BLOCK DATA routine. The user is reminded, however, that

111
English units (ftJ

, psia, Btu, lb, etc.) are not consistent and would

require insertion of conversion factors into each of the property

subroutines.

4) Convergence tolerances and machine precision. The convergence tolerance

for the iteration loops in the subroutines BUBLT and VIT are set by the

value of TOLR in BLOCK DATA. The BUBLT routine contains nested iteration

loops. As the calculation proceeds from the inner to outermost loop the

convergence tolerance must be progressively relaxed/ this is accomplished

by scaling the convergence tolerance of each loop to a multiple of TOLR.

A value of 10 is presently used for TOLR; this value worked well on a

Sperry-Univac 1100/82 mainframe computer with a single precision word of

—3 fi 36 bits (approximately 8 decimal digits of accuracy) and a range of 10

to 10 . This yielded a final precision of about one part in 10 for the

calculation of saturation pressure. With a machine of lesser precision,

a larger value of TOLR may be required. Greater precision can be

obtained by reducing the value of TOLR but this may require converting

the entire set of subroutines to double precision. (The programs are

presented in single precision except for subroutine VIT.)

5) Warning messages. The logical unit to which any warning messages are

written is currently set to 6 in the BLOCK DATA routine; this may be

changed to suit different systems.

112
1 SUBROUTINE FITAB (T.PE, VLE, WE.WP.WL.WV.A.B.PC.VLC, WC)

2 C

3 C THIS SUBROUTINE CALCULATES THE *A* AND 'B* PARAMETERS OF THE

4 C EQUATION OF STATE WHICH BEST FITS A SET OF EXPERIMENTAL SATURATION

5 C DATA FOR A PURE COMPONENT AT A GIVEN TEMPERATURE. THE REQUIRED

6 C DATA ARE SATURATION PRESSURE AND LIQUID AND VAPOR MOLAR

7 C VOLUMES AT TEMPERATURE T. WEIGHTING FACTORS ARE ASSOCIATED WITH

8 C EACH EXPERIMENTAL QUANTITY AND MAY BE ADJUSTED TO ACCOUNT FOR THE

9 C RELATIVE RELIABILITY OF THE VARIOUS EXPERIMENTAL QUANTITIES. THE

10 C INDIVIDUAL A AND B VALUES GIVEN BY THIS ROUTINE WOULD THEN BE FIT

11 C BY A SEPARATE ROUTINE TO AN APPROPRIATE FUNCTION OF TEMPERATURE

12 C FOR USE WITH THE EQUATION OF STATE.

13 C

14 C INPUTS:

15 C T - TEMPERATURE (K)

16 C PE - EXPERIMENTAL SATURATION PRESSURE (KPA)

17 C VLE - EXPERIMENTAL LIQUID PHASE MOLAR VOLUME (M**3/KMOL)

18 C WE - EXPERIMENTAL VAPOR VOLUME (M**3/KMOL)

19 C WP - WEIGHTING FACTOR FOR PRESSURE DATA

20 C WL - WEIGHTING FACTOR FOR LIQUID VOLUME DATA

21 C WV - WEIGHTING FACTOR FOR VAPOR VOLUME DATA

22 C

23 C OUTPUTS:

24 C A - EQUATION OF STATE PARAMETER ASSOCIATED WITH INTERMOLECULAR

25 C ATTRACTION (KJ M**3/KMOL**2)

26 C B - EQUATION OF STATE PARAMETER ASSOCIATED WITH THE MOLECULAR

27 C VOLUME (M**3/KMOL)

28 C PC - SATURATION PRESSURE CALCULATED BY EQUATION OF STATE USING

29 C ABOVE VALUES OF A AND B (KPA)

30 C VLC - CALCULATED LIQUID VOLUME (M**3/KMOL)

31 C WC - CALCULATED VAPOR VOLUME (M**3/KMOL)

32 C

33 C OTHER SUBROUTINES REFERENCED:

34 C PLIMIT - DETERMINES INITIAL BOUNDS ON PRESSURE ITERATION

35 C VIT - ITERATION FOR MOLAR VOLUME GIVEN T, P 36 C

37 C NOTE: IT IS NOT NECESSARY TO REFERENCE BCONST WHEN USING THIS

38 C ROUTINE. ALL NECESSARY COMMON BLOCKS ARE INITIALIZED HERE.

39 C THE GAS CONSTANT IS SET (IN SI UNITS) IN THIS ROUTINE.

40 C

41 IMPLICIT REAL (A-H.O-Z)

42 COMMON /TOL/ TOLR, ITMAX, LUP 43 COMMON /RDATA4/ R 44 DIMENSION AA(2) ,BB(2) ,VL(2,2) , W(2,2) ,P(2,2) ,PG(3)

45 1 ,PL(3),FP(2)

46 LOGICAL LPPOS, LPNEG, LV1CON, LV2CON

47 C

48 C STATEMENT FUNCTIONS FOR GIBBS FREE ENERGY AND THE DERIVATIVE

49 C OF GIBBS FREE ENERGY WITH RESPECT TO THE B PARAMETER

50 C

51 G(T,V,A,B)=R*T*(-LOG(V)+0.25*B/(V-0.25*B)**3

52 1 »((8.0*V-2.25*B)*V+0.1875*B»B))+A/B*LOG(V/(V+B))-A/(V+B)

53 DGDB(T,V,A > B)=R*T/(V-0.25*B)»*3*((-2.25*V+0.375*B)+

54 1 ((2.0«V-0.5625*.B)*V+0.046875*B*B)*(V+0.5*B)/(V-0.25*B))+

55 1 A*((1.0/(V+B)-ALOG(V/(V+B))/B)/B+1 .0/(V+B)**2)

56 R=8.314

57 TOLR=1.0E-7

58 ITMAX=20

59 LUP=6

60 PE2=PE*PE

61 VLE2=VLE«VLE

62 WE2=WE*WE

63 RT=R*T

64 C

65 C USE EXPERIMENTAL LIQUID MOLAR VOLUME AS INITIAL GUESS FOR B AND

66 C ENTER ITERATION TO FIND A, B WHICH EXACTLY SATISFY VOLUME DATA

67 C

68 B=VLE

69 DO 200 IT=1 .100

70 YL=0.25*B/VLE

113
71 YV=0.25*B/VVE

72 A=RT*((1 .0+(1 .0+(1 .0-YL)*YL)*YL)/(VLE*(1 .0-YL)**3)

73 1 -(1 .0+(1 .0+(1 .0-YV)*YV)*YV)/(WE*(1 .0-YV)**3))/

74 1 (1 .0/(VLE*(VLE+B))-1 .0/(WE*(WE+B)))

75 FG=G(T,VLE,A.B)-G(T,WE,A,B)

76 IF (ABS(FG).LT.TOLR) GOTO 240

77 B2=B-FG/(DGDB(T,VLE,A,B)-DGDB(T,WE,A,B))

78 IF (B2.GT.4.0*VLE) B2=(B+4.0*VLE)/2.0

79 IF (B2.LT.0.0) B2=0.5*B

80 B=B2

81 200 CONTINUE

82 240 AA(2)=A

83 BB(2)=B

84 DO 260 1=2,1,-1

85 DO 260 J=3-I,2

86 VL(I,J)=VLE

87 260 VV(I,J)=VVE

88 C 89 C ENTER ITERATION TO FIND A, B WHICH BEST FIT P AND V DATA

90 C

91 DO 400 IT=1 .ITMAX 92 C 93 C FIND VOLUMES AND PRESSURE PREDICTED BY EQUATION OF STATE FOR CURRENT

94 C GUESS OF A, B AND FOR SLIGHTLY DIFFERENT A, B (FOR PURPOSES OF

95 C NUMERICAL DIFFERENTIATION)

96 C 97 AA(1)=1 .001*AA(2)

98 BB(1)=1 .001»BB(2)

99 DO 360 1=2,1 ,-1

100 DO 360 J=3-I,2

101 C

102 C CALL SUBROUTINE TO DETERMINE THE UPPER AND LOWER BOUNDS

103 C ON PRESSURE FOR WHICH THERE ARE BOTH LIQUID AND VAPOR

104 C SOLUTIONS OF THE EQUATION OF STATE

105 C 106 CALL PLIMIT (T.AA( I ) ,BB(J) .VLOW.VUP, PLOW, PUP)

107 C

108 C SET INITIAL GUESSES FOR PRESSURE NEAR THE UPPER AND

109 C LOWER BOUNDS. IF THE LOWER BOUND FOR PRESSURE IS NEGATIVE

110 C RESET IT TO A SMALL POSITIVE VALUE.

111 C 112 IF (PLOW. LE. 0.0) THEN

113 VLOW=0.8*BB(J)

114 TC=AA(I)/(BB(J)*4.398909*R)

115 PC=0.1049995*R*TC/BB(J)

116 PL0W=1 .0E-12*PC

117 PG(1)=PL0W

118 ELSE

119 PG(1)=PLOW+0. 0001* (PUP-PLOW)

120 END IF

121 PG ( 2 )=PUP-0. 0001* (PUP-PLOW)

122 PL(1)=AL0G(PG(1))

123 PL(2)=ALOG(PG(2))

124 VL(I,J)=0.9*VLOW

125 VV(I,J)=1 .1*VUP

126 KP=1

127 LPPOS=. FALSE.

128 LPNEG=. FALSE.

129 C

130 C STARTING WITH INITIAL VALUES OF PRESSURE CLOSE TO THE UPPER

131 C AND LOWER BOUNDS (FOUND BY SUBROUTINE PLIMIT) ITERATE ON

132 C LOG (P) UNTIL THE GIBBS FREE ENERGY OF BOTH PHASES ARE EQUAL.

133 C A COMBINATION OF SECANT AND REGULI-FALSI METHODS IS USED

134 C FOR THE ITERATION.

135 C 136 C ENTER ITERATION FOR SATURATION PRESSURE

137 C

138 DO 300 ITP=1 , ITMAX

139 LV1C0N=. FALSE.

140 LV2C0N=. FALSE.

114
141 C EVALUATE VOLUMES FOR CURRENT VALUES OF A, B. P 142 CALL VIT (T ,PG(KP) , AA( I ) ,BB( J) , VL( I , J) . . TRUE. , LV1CON)

143 CALL VIT (T,PG(KP) ,AA( I ) ,BB(J) , W( I . J) , . FALSE. , LV2CON)

144 C PRESSURE ITERATION HAS CONVERGED WHEN GIBBS FREE ENERGY OF

145 C LIQUID AND VAPOR PHASES ARE EQUAL

146 FP(KP)=G(T ( VL(I > J),AA(I),BB(J))-G(T.W(I,J),AA(I),BB(J))

147 IF (FP(KP).LT.0.0) THEN

148 LPNEG=.TRUE.

149 FPNEG=FP(KP)

150 PNEG=PL(KP)

151 ELSE

152 LPPOS=.TRUE.

153 FPPOS=FP(KP)

154 PPOS=PL(KP)

155 END IF 156 IF (ITP.LE.1) THEN

157 KP=2

158 ELSE

159 DGDPL=(FP(2)-FP(1))/(PL(2)-PL(1))

160 IF (ABS(FP(KP)/(PL(KP)*DGDPL)).LT.TOLR) GOTO 340

161 C NEXT GUESS FOR LOG (P) GIVEN BY SECANT METHOD

162 PL(3)=PL(2)-FP(2)/DGDPL

163 C IF NEXT GUESS FOR LOG (P) IS FURTHER FROM SOLUTION THAN

164 C PREVIOUS BEST GUESS, USE REGULI-FALSI METHOD FOR NEXT GUESS

165 IF ((PL(3).GT.MAX(PNEG,PPOS) .OR. 166 1 PL(3) .LT.MIN(PNEG.PPOS)) .AND. LPNEG .AND. LPPOS)

167 1 PL(3)=PPOS-FPPOS*(PPOS-PNEG)/(FPPOS-FPNEG)

168 PL(1)=PL(2)

169 PL(2)=PL(3)

170 FP(1)=FP(2)

171 PG(2)=EXP(PL(2))

172 END IF 173 300 CONTINUE

174 340 CONTINUE

175 P(I,J)=PG(KP)

176 360 CONTINUE

177 C

178 C EVALUATE DERIVATIVE OF LIQUID AND VAPOR VOLUMES AND SATURATION

179 C PRESSURE WITH RESPECT TO THE A AND B PARAMETERS AND SOLVE SYSTEM

180 C OF EQUATIONS TO ARRIVE AT NEW GUESSES FOR A, B

181 C

182 DA=AA(2)-AA(1)

183 DB=BB(2)-BB(1)

184 DVLDA=(VL(2,2)-VL(1,2))/DA

185 DVLDB=(VL(2,2)-VL(2.1))/DB

186 DWDA=(VV(2,2)-W(1 ,2))/DA

187 DWDB=(W(2,2)-W(2,1))/DB

188 DPDA=(P(2,2)-P(1 ,2))/DA

189 DPDB=(P(2,2)-P(2,1))/DB

1 90 Q1 1=WL* (DVLDA/VLE) *»2+WV* (DWDA/WE) **2+WP* (DPDA/PE) *»2

191 Q1 2=WL»DVLDA*DVLDB/VLE2+WV*DVVDA*DWDB/WE2+WP*DPDA*DPDB/PE2

1 92 Q22=WL» (DVLDB/VLE) *»2+WV* (DWDB/WE) * »2+WP* (DPDB/PE) **2

193 C1=WL*DVLDA*(VL(2,2)-VLE)/VLE2+WV*DWDA*(W(2 > 2)-WE)/VVE2+

194 1 WP*DPDA*(P(2,2)-PE)/PE2

1 95 C2=WL*DVLDB* (VL(2 , 2)-VLE)/VLE2+WV*DWDB* ( W(2 , 2)-VVE)/VVE2+

196 1 WP»DPDB*(P(2.2)-PE)/PE2

197 DET=Q1 1 *Q22-Q1 2*Q12

198 DELA=(C1*Q22-C2*Q12)/DET

199 DELB=(C2*Q11-C1»Q12)/DET

200 IF (ABS(DELA).GT. 0. 2*ABS(AA(2) )

) THEN

201 AA(2)=AA(2)-0.2*SIGN(AA(2),DELA)

202 ELSE

203 AA(2)=AA(2)-DELA

204 END IF 205 IF (ABS(DELB).GT. 0.2*ABS(BB(2))) THEN

206 BB(2)=BB(2)-0.2*SIGN(BB(2),DELB)

207 ELSE

208 BB(2)=BB(2)-DELB

209 END IF 210 C ITERATION HAS CONVERVED WHEN SUCCESSIVE GUESSES FOR A AND B

115
211 C ARE WITHIN A CONVERGENCE TOLERANCE.

212 IF (ABS(DELA/AA(2)).LT.10.0*TOLR .AND. 213 1 ABS(DELB/BB(2)).LT.10.0*TOLR) THEN

214 A=AA(2)

215 B=BB(2)

216 PC=P(2,2)

217 VLC=VL(2,2)

218 WC=W(2,2)

219 RETURN

220 END IF

221 400 CONTINUE

222 WRITE (LUP.1200)

223 A=AA(2)

224 B=BB(2)

225 PC=P(2.2)

226 VLC=VL(2,2)

227 WC=W(2,2)

228 RETURN

229 1200 FORMAT (1X, ' ITERATION FOR A AND B DID NOT CONVERGE')

230 END

116
1 SUBROUTINE FITF (T,XL ( PE.VLE,VVE,XVE,WP,WL,WV,WX, F.PC.VLC,

2 1 VVC.XVC)

3 C

4 C THIS ROUTINE DETERMINES THE INTERACTION PARAMETER WHICH BEST FITS

5 C A SET OF EXPERIMENTAL SATURATION DATA AT A GIVEN TEMPERATURE AND

6 C LIQUID COMPOSITION. THE INDIVIDUAL VALUES OF THE INTERACTION

7 C PARAMETER DETERMINED BY THIS SUBROUTINE MAY THEN BE AVERAGED OVER

8 C MULTIPLE DATA SETS OR FIT AS A FUNCTION OF TEMPERATURE EITHER

9 C MANUALLY OR IN THE CALLING PROGRAM. WEIGHTING FACTORS ARE

10 C ASSOCIATED WITH EACH EXPERIMENTAL QUANTITY. IF A PARTICULAR

11 C QUANTITY IS UNAVAILABLE, THE CORRESPONDING WEIGHTING FACTOR SHOULD

12 C BE SET TO ZERO. A MINIMUM OF ONE MEASURED QUANTITY (IN ADDITION TO

13 C TEMPERATURE AND LIQUID COMPOSITION) IS REQUIRED FOR THE

14 C DETERMINATION OF THE INTERACTION PARAMETER.

15 C

16 C INPUTS:

17 C T - TEMPERATURE (K)

18 C XL - LIQUID PHASE COMPOSITION (MOLE FRACTION)

19 C PE - EXPERIMENTAL SATURATION PRESSURE (KPA)

20 C VLE - EXPERIMENTAL LIQUID MOLAR VOLUME (M**3/KM0L)

21 C WE - EXPERIMENTAL VAPOR MOLAR VOLUME (M**3/KMOL)

22 C XVE - EXPERIMENTAL VAPOR PHASE VOMPOSITION (MOL FRAC)

23 C WP - WEIGHTING FACTOR FOR PRESSURE DATA

24 C WL - WEIGHTING FACTOR FOR LIQUID VOLUME DATA

25 C WV - WEIGHTING FACTOR FOR VAPOR VOLUME DATA

26 C WX - WEIGHTING FACTOR FOR VAPOR COMPOSITION DATA

27 C

28 C OUTPUTS:

29 C F - VALUE OF THE INTERACTION PARAMETER WHICH BEST FIT THE

30 C EXPERIMENTAL DATA (DIMENSIONLESS)

31 C PC - SATURATION PRESSURE CALCULATED BY EQUATION OF STATE

32 C USING ABOVE VALUE OF F (KPA)

33 C VLC - CALCULATED LIQUID VOLUME (M**3/KMOL)

34 C VVC - CALCULATED VAPOR VOLUME (M**3/KMOL)

35 C XVC - CALCULATED VAPOR COMPOSITION (MOL FRAC)

36 C

37 C OTHER SUBROUTINES REFERENCED:

38 C ESPAR - COMPUTATION OF EQUATION OF STATE PARAMETERS

39 C BUBLT - CALCULATE SATURATION CONDITIONS AT GIVEN T. XL

40 C BUBLT IN TURN REFERENCES:

41 C PLIMIT - DETERMINES INITIAL BOUNDS ON PRESSURE ITERATION

42 C VIT - ITERATION FOR MOLAR VOLUME GIVEN T, P 43 C ZXLSF - ONE-DIMENSIONAL MINIMIZATION ROUTINE CONTAINED IN 44 C PROPRIETARY IMSL MATH/STATISTICS LIBRARY. AN

45 C EQUVILENT ROUTINE FROM THE USER'S COMPUTER CENTER

46 C MAY BE USED BY MODIFYING LINES 69 - 85 AND THE EXTERNAL

47 C FUNCTION SFUNC AS APPROPRIATE.

48 C

49 C NOTE: THE SUBROUTINE BCONST MUST BE CALLED ONCE BY THE MAIN

50 C PROGRAM FOR EACH MIXTURE UNDER CONSIDERATION

51 C 52 C COMMON BLOCK DPASS IS USED TO PASS INFORMATION TO THE EXTERNAL

53 C FUNCTION SFUNC WHICH CALCULATED THE EXPRESSION TO BE MINIMIZED.

54 C

55 COMMON /DPASS/ TP,XLP,PEP,VLEP,VVEP.XVEP.WPP,WLP,WVP,WXP,

56 1 PCP.VLCP.VVCP.XVCP

57 EXTERNAL SFUNC

58 DATA ST EP/0 . 004/ , SBOUND/0 . 4/ , SACC/ 1 . E-6/ , MAX FN/40/

59 C FILL UP COMMON BLOCK FROM ARGUMENT LIST OF SUBROUTINE

60 TP=T

61 XLP=XL

62 PEP=PE

63 VLEP=VLE

64 WEP=VVE

65 XVEP=XVE

66 WPP=WP

67 WLP=WL

68 WVP=WV

69 WXP=WX

70 F=0 . 04

117
71 IERR=0

72 C

73 C PARAMETERS FOR THE MINIMIZATION ROUTINE:

74 C SFUNC - EXTERNAL FUNCTION SUBROUTINE WHICH EVALUATES THE

75 C EXPRESSION TO BE MINIMZED AS A FUNCTION OF F 76 C F - VALUE OF THE INDEPENDENT VARIABLE (IN THIS CASE THE

77 C INTERACTION PARAMETER) WHICH MINIMIZES SFUNC

78 C STEP - INITIAL INCREMENT BY WHICH TO CHANGE F IN THE ITERATION

79 C TO FIND MINIMUM

80 C SBOUND - MAXIMUM AMOUNT THAT F CAN BE VARIED ABOVE OR BELOW

81 C INITIAL VALUE

82 C SACC - ABSOLUTE ACCURACY TO WHICH F IS TO BE DETERMINED

83 C MAXFN - MAXIMUM NUMBER OF ALLOWED CALLS TO SFUNC

84 C I ERR - ERROR FLAG; A VALUE > 128 INDICATES AN ERROR

85 C 86 CALL ZXLSF (SFUNC, F, STEP, SBOUND, SACC, MAXFN, I ERR)

87 C 88 C TRANSFER CALCULATED QUANTITES FROM COMMON BLOCK

89 C 90 PC=PCP

91 VLC=VLCP

92 VVC=WCP

93 XVC=XVCP

94 END 95 C 96 C 97 FUNCTION SFUNC (F)

98 C 99 C THIS FUNCTION EVALUATES THE SUM OF SQUARES DEVIATION BETWEEN

100 C THE EXPERIMENTAL AND CALCULATED QUANTITIES

101 C

102 LOGICAL LCRIT.LTRUE

103 COMMON /DPASS/ TP,XLP,PE,VLE,VVE,XVE,WP,WL,WV,WX,

104 1 PCP,VLCP,VVCP,XVCP

105 COMMON /RDATA1/ AA0.AA1 ,AA2,AB0,AB1 .AB2.BA0.BA1 ,BA2,

106 1 BB0.BB1 ,BB2,F0,F1

107 LTRUE=.TRUE.

108 F0=F

1 09 T=TP

110 XL=XLP

111 C

112 C CALL SUBROUTINE ESPAR TO REEVALUATE A, B PARAMTERS USING CURRENT

113 C GUESS FOR F 114 C 115 CALL ESPAR (-1 .T.XL.A.B)

116 C COMPUTE SATURATION PROPERTIES AT GIVEN T, XL

117 CALL BUBLT (T ,XL,XVC,PC,VLC,VVC, LTRUE, LCRIT)

118 SFUNC=WP*((PE-PC)/PC)**2+WL*((VLE-VLC)/VLC)*»2+

119 1 WV»((VVE-VVC)/VVC)"2+WX*(XVE-XVC)»*2

1 20 PCP=PC

121 VLCP=VLC

1 22 WCP=VVC

1 23 XVCP=XVC

124 RETURN

125 END

118
1 SUBROUTINE BCONST (IR1

,

IR2.F0.F1 ) 2 C

3 C THIS ROUTINE ACCESSES THE CURVE FIT COEFFICIENTS TO THE EQUATION

4 C OF STATE PARAMETERS (STORED IN BLOCK DATA BDESC) FOR THE

5 C REFRIGERANT PAIR OF INTEREST. THE REFERENCE STATES FOR ENTHALPY

6 C AND ENTROPY ARE ALSO COMPUTED. THIS SUBROUTINE MUST BE CALLED

7 C BEFORE ANY OTHER PROPERTY ROUTINES ARE REFERENCED AND ALSO IF

8 C THE MIXTURE OR THE VALUES OF THE INTERACTION COEFFICIENTS F0, F1 9 C ARE CHANGED.

10 C

11 C INPUTS:

12 C IR1. IR2 - CODE NUMBERS FOR PURE COMPONENTS

13 C F0, F1 - COEFFICIENTS TO THE CURVE FIT FOR THE INTERACTION

14 C PARAMETER: F = F0 + F1 *T

15 C

16 C OUTPUTS (VIA COMMON BLOCKS):

17 C A - ARRAY OF A COEFFICIENTS FOR THE TWO PURE COMPONENTS

18 C B - ARRAY OF B COEFFICIENTS FOR THE PURE COMPONENTS

19 C CP - ARRAY OF PURE COMPONENT CP0 COEFFICIENTS

20 C HR - TWO ELEMENT ARRAY OF PURE COMPONENT REFERENCE

21 C ENTHALPIES; THESE ARE EQUAL TO THE SATURATED LIQUID

22 C ENTHALPY AT THE REFERENCE TEMPERATURE MINUS THE PERFECT

23 C GAS ENTHALPY AT THE REFERENCE TEMPERATURE

24 C SR - REFERENCE ENTROPIES; EQUAL TO THE DIFFERENCE BETWEEN

25 C THE SATURATED LIQUID AND PERFECT GAS ENTROPIES AT THE

26 C REFERENCE TEMPERATURE

27 C TC - PURE COMPONENT CRITICAL TEMPERATURES

28 C TREF - REFERENCE TEMPERATURES AT WHICH HR AND SR ARE COMPUTED

29 C WM - PURE COMPONENT MOLECULAR WEIGHTS

30 C

31 C OTHER SUBROUTINES REFERENCED:

32 C BUBLT - COMPUTE SATURATED LIQUID AND VAPOR CONDITIONS

33 C HCVCP - COMPUTE ENTHALPY AT REFERENCE STATE

34 C ENTROP - COMPUTE REFERENCE ENTROPY

35 C 36 IMPLICIT REAL (A-H.O-Z)

37 DIMENSION COEFF(9 ,20) ,CRIT(5,20) ,WM(2) ,TC(2) ,A(3,2) .6(3,2)

,

38 1 CP(3,2),HR(2),SR(2),VR(2),TREF(2),HZERO(20),SZERO(20)

39 INTEGER IR(2)

40 CHARACTER*6 HREF(20)

41 COMMON /ESDATA/ COEFF.CRIT

42 COMMON /HREF1/ HREF

43 COMMON /RDATA1/ A.B.FF0.FF1

44 COMMON /RDATA2/ WM.TC

45 COMMON /CPDATA/ CP

46 COMMON /REF/ TREF.HR.SR.VR

47 COMMON /HSZERO/ HZERO.SZERO

48 IR(1)=IR1

49 IR(2)=IR2

50 IR1=ABS(IR1)

51 IR2=ABS(IR2)

52 FF0=F0

53 FF1=F1

54 DO 100 KR=1 ,2

55 IF (IR(KR).GT.0) THEN

56 C IF IR IS NEGATIVE READ IN COEFFICIENTS FROM DATA FILE

57 ELSE IF (IR(1).EQ.ABS(IR(2)) .AND. KR.EQ.2) THEN

58 HREF(IR1)=HREF(IR2)

59 DO 66 KC=1 ,5

60 66 CRIT(KC,IR2)=CRIT(KC,IR1)

61 HZERO(IR2)=HZERO(IR1)

62 SZERO(IR2)=SZERO(IR1)

63 DO 68 KC=1 ,9

64 68 COEFF(KC.IR2)=COEFF(KC,IR1)

65 ELSE

66 IRK=ABS(IR(KR))

67 READ (*,1000) HREF(IRK)

68 READ (»,*) (CRIT(KC,IRK),KC=1 ,5),

69 1 HZERO(IRK),SZERO(IRK),

70 1 (COEFF(KC,IRK),KC=1 ,9)

119
71 END IF 72 IR(KR)=ABS(IR(KR))

73 WM(KR)=CRIT(1 , IR(KR))

74 TREF(KR)=CRIT(2,IR(KR))

75 TC(KR)=CRIT(3.IR(KR))

76 HR(KR)=0.0

77 SR(KR)=0.0

78 DO 100 KC=1,3

79 A(KC,KR)=COEFF(KC,IR(KR))

80 B(KC,KR)=C0EFF(KC+3,IR(KR))

81 100 CP(KC,KR)=COEFF(KC+6,IR(KR))

82 CALL ESPAR (-2 ,TREF(1 ) ,0.0,AB,BB)

83 C 84 C CALL BUBBLE POINT ROUTINE TO CALCULATE SATURATED LIQUID AND VAPOR

85 C VOLUMES AND THEN CALL ENTHALPY AND ENTROPY ROUTINE TO DETERMINE

86 C REFERENCE VALUES. THE HZERO AND SZERO ALLOW AN ARBITRARY VALUE

87 C TO BE ASSIGNED TO THE SATURATED LIQUID H OR S AT THE REFERENCE

88 C TEMPERATURE.

89 C 90 CALL BUB LT (TREF(1 ), 1 .0 ( XV,P,VR(1 ) ,VV, . TRUE. .. FALSE.

)

91 CALL BUB LT (TREF(2) ,0.0,XV,P,VR(2) ,W. .TRUE. ,. FALSE.

)

92 DO 160 1=1 ,2

93 XL=FLOAT(2-I)

94 CALL HCVCP (1 ,TREF(I ) ,VR( I ) ,XL.HR( I ) .CV.XCP)

95 HR(I)=HR(I)-HZERO(IR(I))

96 SR( I )=ENTROP(TREF( I ) , VR( I ) ,XL)-SZERO( IR( I )

)

97 160 CONTINUE

98 RETURN

99 1000 FORMAT (A6)

100 END

120
1 SUBROUTINE BUBLT (T.XL.XV.P.VL.VV, LBUB, LCRIT)

2 C

3 C GIVEN TEMPERATURE AND COMPOSITION OF ONE PHASE THIS ROUTINE

4 C CALCULATES THE SATURATION PRESSURE, THE COMPOSITION OF THE OTHER

5 C PHASE AND THE LIQUID AND VAPOR MOLAR VOLUMES.

6 C

7 C INPUTS: , 8 C T - TEMPERATURE (K)

9 C ONLY ONE OF: XL - LIQUID COMPOSITION (MOLE FRACTION)

10 C OR: XV - VAPOR COMPOSITION (MOLE FRACTION)

11 C LBUB - LOGICAL VARIABLE

12 C IF LBUB=.TRUE. LIQUID COMPOSITION IS GIVEN (COMPUTE

13 C BUBBLE POINT)

14 C IF LBUB=. FALSE. VAPOR COMPOSITION IS GIVEN (COMPUTE

15 C DEW POINT)

1 6 C

17 C OUTPUTS:

18 C XL OR XV - COMPOSITION OF CALCULATED PHASE

19 C P - SATURATION PRESSURE (KPA)

20 C VL - LIQUID MOLAR VOLUME (M**3/KM0L)

21 C VV - VAPOR MOLAR VOLUME (M«.*3/KMOL)

22 C LCRIT - ERROR FLAG; IF LCRIT=.TRUE. THE INPUT TEMPERATURE

23 C EXCEEDS THE CRITICAL TEMPERATURE OF THE PURE COMPONENT

24 C OR THE PSEUDO-PURE COMPONENT CORRESPONDING TO THE

25 C MIXTURE COMPOSITION AND NO CALCULATIONS ARE DONE.

26 C 27 C OTHER SUBROUTINES REFERENCED:

28 C VIT - ITERATION FOR MOLAR VOLUME

29 C PLIMIT - DETERMINES INITIAL BOUNDS ON PRESSURE AND VOLUME

30 C ESPAR - COMPUTATION OF EQUATION OF STATE PARAMETERS

31 C

32 C GENERAL NOMENCLATURE FOR FIRST LETTER OF VARIABLE NAMES

33 C A,B - EQUATION OF STATE PARAMETERS

34 C F - MIXING PARAMETER

35 C T - TEMPERATURE

36 C P - PRESSURE

37 C V - MOLAR VOLUME

38 C X - COMPOSITION

39 C G - GIBBS FREE ENERGY

40 C U - CHEMICAL POTENTIAL

41 C Y - COMBINATION OF VARIABLES USED IN EQUATION OF STATE

42 C TOL - CONVERGENCE TOLERANCE FOR ITERATION LOOPS

43 C I, J - INDEX VARIABLES FOR ITERATION AND DO LOOPS

44 C L - LOGICAL VARIABLES SUCH AS NON-CONVERGENCE FLAGS

45 C 46 C GENERAL NOMENCLATURE FOR SECOND OR THIRD LETTER OF VARIABLES

47 C A,B - COMPONENTS OF MIXTURE; COMPOSITION IS MOLE FRACTION A 48 C L - LIQUID PHASE

49 C V - VAPOR PHASE

50 C 1 - PARENT PHASE (PHASE WITH SPECIFIED COMPOSITION)

51 C 2 - INCIPIENT PHASE

52 C (FOR EXAMPLE UA1 REFERS TO CHEMICAL POTENTIAL OF COMPONENT A

53 C IN PHASE 1)

54 C 55 C 56 IMPLICIT REAL (A-H.O-Z)

57 LOGICAL LBUB, LCRIT, LV ICON, LV2CON, LXCON, LXPOS, LXNEG. 58 1 LPPOS , LPNEG , LPPCON

59 COMMON /ESPAR1/ AA.AB.BA.BB, F.C1 ,D1 ,C2,D2

60 COMMON /RDATA4/ R

61 COMMON /TOL/ TOLR, ITMAX, LUP 62 DIMENSION PP(3) , FP(2) ,XX2(3) , FX2(2) ,PL(3)

63 C

64 C STATEMENT FUNCTIONS FOR GIBBS FREE ENERGY AND CHEMICAL POTENTIAL

65 C NOTE THAT SINCE ONLY DIFFERENCES OF G AND U ARE USED IN THE PROGRAM

66 C ANY TERMS WHICH WOULD CANCEL ARE OMITTED. THE EXPRESSION FOR U 67 C IS DIVIDED BY RT TO OBTAIN A DIMENSIONLESS QUANTITY.

68 C

69 G(T,V,A,B)=R»T*(-LOG(V)+0.25*B«((8.0»V-2.25*B)*V+0.1875*B»B)

70 1 /(V-0.25*B)»*2/(V-0.25*B))+A/B*LOG(V/(V+B))-A/(V+B)

121
71 U(T,X,V.A,B,AA,AB,BI , F,Y)=(Y*(4.0-3.0*Y)+BI*(4.0-2 .0»Y)*Y/

72 1 (B*(1 .0-Y)))/(1 .0-Y)**2+(BI»A*LOG(1 .0+B/V)/B-BI»A/(V+B)

73 1 +2.*(X*AA+(1 .0-F)»(1 .0-X)*SQRT(AA*AB))»LOG(V/(V+B)))/(R*T*B)

74 1 -LOG(V)

75 C

76 LCRIT=.FALSE.

77 C

78 C COMPUTE PURE COMPONENT E.S. COEFFICIENTS, THE MIXING PARAMETER,

79 C AND THE E.S. COEFFICIENTS FOR PHASE 1 80 C

81 IF (LBUB) THEN

82 X1=XL

83 XV=XL

84 ELSE

85 X1=XV

86 XL=XV

87 END IF 88 XB1=1.0-X1

89 CALL ESPAR (0.T.X1 ,A1 ,B1

)

90 C

91 C DETERMINE IF INPUT TEMPERATURE EXCEEDS CRITICAL POINT;

92 C IF SO, SET ERROR FLAG AND RETURN

93 C 94 TC=A1/(B1*4.398909*R)

95 IF (T.GT.0.99*TC) THEN

96 LCRIT=.TRUE.

97 WRITE (LUP.1010)

98 RETURN

99 END IF 100 C

101 C ENTER ITERATION FOR PSEUDO-PURE COMPONENT. THIS ITERATION

102 C YIELDS THE FINAL RESULT FOR A PURE COMPONENT AND PROVIDES

103 C A STARTING GUESS FOR THE PRESSURE OF A MIXTURE

104 C

105 C CALL SUBROUTINE TO DETERMINE THE UPPER AND LOWER BOUNDS

106 C ON PRESSURE FOR WHICH THERE ARE BOTH LIQUID AND VAPOR

107 C SOLUTIONS OF THE EQUATION OF STATE

108 C

109 CALL PLIMIT (T.A1 ,B1 .VLOW.VUP, PLOW, PUP)

110 C

111 C SET INITIAL GUESSES FOR PRESSURE NEAR THE UPPER AND

112 C LOWER BOUNDS. IF THE LOWER BOUND FOR PRESSURE IS NEGATIVE

113 C RESET IT TO A SMALL POSITIVE VALUE.

114 C 115 IF (PLOW. LE. 0.0) THEN

116 VLOW=0.8*B1

117 PC=0.1049995*R*TC/B1

118 PLOW=1 .0E-12*PC

119 PP(1)=PLOW

120 ELSE

121 PP(1)=PLOW+0. 0001* (PUP-PLOW)

122 END IF 123 PP (2 )=PUP-0. 0001* (PUP-PLOW)

124 PL(1)=LOG(PP(1))

125 PL(2)=LOG(PP(2))

126 VL=0.9*VLOW

127 VV=1.1*VUP

128 J=1

129 LPPOS=. FALSE.

130 LPNEG=. FALSE.

131 LPPCON=. FALSE.

132 C 133 C STARTING WITH INITIAL VALUES OF PRESSURE CLOSE TO THE UPPER

134 C AND LOWER BOUNDS (FOUND BY SUBROUTINE PLIMIT) ITERATE ON

135 C LOG (P) UNTIL THE GIBBS FREE ENERGY OF BOTH PHASES ARE EQUAL.

136 C A COMBINATION OF SECANT AND REGULI-FALSI METHODS IS USED

137 C FOR THE ITERATION.

138 C

139 DO 400 IT=1 .ITMAX

140 LV1CON=. FALSE.

122
141 LV2C0N=. FALSE.

142 CALL VIT (T,PP( J) ,A1 ,B1 , VL, . TRUE. . LV1C0N)

143 CALL VIT (T,PP(J) ,A1 ,B1 ,W, . FALSE. . LV2C0N)

144 GL=G(T,VL,A1 ,B1)

145 GV=G(T,W.A1 ,B1)

146 FP(J)=GL-GV

147 IF (FP(J).LT.0.0) THEN

148 LPNEG=.TRUE.

149 FPNEG=FP(J)

150 PNEG=PL(J)

151 ELSE

152 LPPOS=.TRUE.

153 FPPOS=FP(J)

154 PPOS=PL(J)

155 END IF 156 IF (IT.LE.1) THEN

157 J=2

158 ELSE

159 DGDPL=(FP(2)-FP(1))/(PL(2)-PL(1))

160 IF (DGDPL.EQ.0.0) GOTO 440

161 IF (ABS(FP(J)/(PL(J)*DGDPL)).LT.TOLR) GOTO 440

162 C NEXT GUESS FOR LOG (P) GIVEN BY SECANT METHOD

163 PL(3)=PL(2)-FP(2)/DGDPL

164 C IF NEXT GUESS FOR LOG (P) IS FURTHER FROM SOLUTION THAN

165 C PREVIOUS BEST GUESS, USE REGULI-FALSI METHOD FOR NEXT GUESS

166 IF ((PL(3).GT.MAX(PNEG,PPOS) .OR. 167 1 PL(3).LT.MIN(PNEG,PP0S)) .AND. LPNEG .AND. LPPOS)

168 1 PL(3)=PPOS-FPPOS*(PPOS-PNEG)/(FPPOS-FPNEG)

169 PL(1)=PL(2)

170 PL(2)=PL(3)

171 FP(1)=FP(2)

172 PP(2)=EXP(PL(2))

173 END IF 174 400 CONTINUE

175 C IF ITERATION HAS NOT CONVERGED, SET ERROR FLAG. 176 LPPCON=.TRUE.

177 C 178 C END OF PSEUDO-PURE COMPONENT ITERATION

179 C 180 C FOR A PURE COMPONENT THE ABOVE ITERATION GIVES THE FINAL RESULT

181 C

182 440 IF (X1*XB1.LE.TOLR) THEN

183 IF (LV1CON) WRITE (LUP.1050)

184 IF (LV2CON) WRITE (LUP.1055)

185 IF (LPPCON) WRITE (LUP.1020)

186 P=PP(J)

187 RETURN

188 END IF 189 C 190 C ENTER ITERATION FOR MIXTURE

191 C

192 C THE MIXTURE ITERATION CONSISTS OF TWO CONCENTRIC ITERATION

193 C LOOPS WHICH VARY THE SATURATION PRESSURE OF THE MIXTURE AND THE

194 C COMPOSITION OF THE COMPUTED PHASE TO GIVE EQUAL CHEMICAL

195 C POTENTIALS FOR EACH OF THE COMPONENTS BETWEEN THE TWO PHASES.

196 C THE INITIAL GUESS FOR THE PRESSURE IS GIVEN BY THE PSEUDO-PURE

197 C ITERATION ABOVE; THE INITIAL GUESS FOR COMPOSITION IS THAT X2=X1

198 C

199 C ASSIGN INITIAL VALUES OF LIQUID AND VAPOR VOLUMES FROM ABOVE

200 C ITERATION TO PHASE 1 AND 2 VOLUMES.

201 C 202 IF (LBUB) THEN

203 V1=VL

204 V2=VV

205 ELSE

206 V1=VV

207 V2=VL

208 END IF 209 PP(1)=PP(J)

210

123
211 C BEGIN ITERATION FOR SATURATION PRESSURE OF MIXTURE

212 C

213 J=1

214 X2CONV=X1

215 LPNEG=. FALSE.

216 LPPOS=. FALSE.

217 DO 800 ITP=1 .ITMAX 218 XX2(1)=X2CONV

219 LXCON=. FALSE.

220 LV1CON=. FALSE.

221 CALL VIT (T,PP(J) ,A1 ,B1 ,V1 , LBUB, LV1CON)

222 C

223 C IF VOLUME ITERATION HAS NOT CONVERGED, TRY A NEW PRESSURE AND

224 C RETURN TO THE BEGINNING OF THE ITERATION

225 C

226 IF (LV1CON .OR. LXCON) THEN

227 PP(2)=0.5*(PP(1)+PP(2))

228 GOTO 800

229 END IF 230 C COMPUTE CHEMICAL POTENTIALS FOR PHASE 1 231 Y=0.25*B1/V1

232 UA1=U(T,X1 ,V1 ,A1 ,B1 .AA.AB.BA, F.Y)

233 UB1=U(T,XB1 ,V1 ,A1 ,B1 .AB.AA.BB, F.Y)

234 C

235 C ENTER INNER ITERATION LOOP (FOR COMPOSITION OF PHASE 2)

236 C

237 JJ=1

238 LXNEG=. FALSE.

239 LXPOS=. FALSE.

240 DO 600 IT=1 .ITMAX

241 LV2CON=. FALSE.

242 XB2=1 .0-XX2(JJ)

243 C COMPUTE EQUATION OF STATE COEFFICIENTS FOR PHASE 2 244 CALL ESPAR (0,T,XX2( JJ) ,A2,B2)

245 CALL VIT (T,PP( J) .A2.B2.V2, .NOT. LBUB, LV2CON)

246 C

247 C IF VOLUME ITERATION HAS NOT CONVERGED, TRY A NEW PRESSURE

248 C AND RETURN TO THE START OF THE PRESSURE ITERATION.

249 C

250 IF (LV2CON) THEN

251 PP(2)=0.5*(PP(1)+PP(2))

252 GOTO 800

253 END IF 254 C COMPUTE CHEMICAL POTENTIALS OF PHASE 2 255 Y=0.25*B2/V2

256 UA2=U(T,XX2(JJ),V2,A2,B2,AA,AB,BA,F,Y)

257 UB2=U ( T , XB2 , V2 , A2 , B2 , AB , AA , BB , F , Y)

258 C

259 C CALCULATE THE COMPOSITION OF PHASE 2 FROM THE COMPOSITION

260 C OF PHASE 1 AND THE CHEMICAL POTENTIALS. THE INNER ITERATION

261 C LOOP HAS CONVERGED WHEN THE CALCULATED COMPOSITION EQUALS

262 C (WITHIN A CONVERGENCE TOLERANCE) THE GUESSED VALUE OF X2.

263 C

264 ZA=X 1 * EXP ( UA 1 -UA2

265 ZB=XB1»EXP(UB1-UB2)

266 C=ZA+ZB

267 X2CALC=ZA/C

268 FX2(JJ)=X2CALC-XX2(JJ)

269 IF (ABS(FX2(JJ)).LT.TOLR) THEN

270 X2CONV=XX2(JJ)

271 GOTO 640

272 END IF 273 C UPDATE POSITIVE OR NEGATIVE BOUNDS FOR USE WITH REGULI-FALSI

274 IF (FX2(JJ).LT.0.0) THEN

275 LXNEG=.TRUE.

276 FXNEG=FX2(JJ)

277 XNEG=XX2(JJ)

278 ELSE

279 LXPOS= . TRUE

280 FXPOS=FX2(JJ)

124
281 XP0S=XX2(JJ)

282 END IF 283 C 284 C UPDATE THE GUESS FOR X2. THE COMPOSITION COMPUTED ABOVE IS 285 C USED FOR THE SECOND GUESS. A COMBINATION OF SECANT AND

286 C REGULI-FALSI METHODS IS USED FOR THIRD AND SUBSEQUENT GUESSES.

287 C 288 IF (IT.LE.1) THEN

289 JJ=2

290 XX2(2)=X2CALC

291 ELSE

292 IF (FX2(1).EQ.FX2(2)) THEN

293 X2CONV=XX2(JJ)

294 GOTO 640

295 END IF 296 XX2(3)=XX2(2)-FX2(2)*(XX2(2)-XX2(1))/(FX2(2)-FX2(1))

297 IF (LXPOS .AND. LXNEG) THEN

298 IF (XX2(3).LT.MIN(XPOS,XNEG) .OR. XX2(3) .GT.MAX(XPOS.XNEG)) THEN

299 XX2(3)=XPOS-FXPOS*(XPOS-XNEG)/(FXPOS-FXNEG)

300 END IF 301 END IF 302 XX2(1)=XX2(2)

303 XX2(2)=XX2(3)

304 FX2(1)=FX2(2)

305 END IF 306 XX2(JJ)=MIN(1 .0,MAX(0.0,XX2(JJ)))

307 600 CONTINUE

308 C IF INNER ITERATION LOOP HAS NOT CONVERGED. SET ERROR FLAG

309 LXCON=.TRUE.

310 C

311 C END OF ITERATION LOOP FOR PHASE 2 COMPOSITION

312 C 313 640 FP(J)=1 .0-C

314 C 315 C OUTER (PRESSURE) ITERATION HAS CONVERGED WHEN C = 1.000

316 C (I.E. WHEN THE CHEMICAL POTENTIALS OF EACH COMPONENT ARE

317 C THE SAME IN BOTH PHASES).

318 C 319 IF (ABS(FP(1)).LT.100.*TOLR) GOTO 840

320 C

321 C PROVIDED THAT THE X2 ITERATION HAS CONVERGED FOR THE CURRENT

322 C GUESS OF PRESSURE, UPDATE THE POSITIVE AND NEGATIVE

323 C BOUNDS FOR USE WITH THE REGULI-FALSI METHOD.

324 C 325 IF (.NOT.LXCON) THEN

326 IF (FP(J).LT.0.0) THEN

327 LPNEG=.TRUE.

328 FPNEG=FP(J)

329 PNEG=PP(J)

330 ELSE

331 LPPOS=.TRUE.

332 FPPOS=FP(J)

333 PPOS=PP(J)

334 END IF 335 END IF 336 C 337 C COMPUTE NEW GUESS FOR SATURATION PRESSURE.

338 C 339 IF (ITP.LE.2 .OR. FP(1 ) . EQ. FP(2)) THEN

340 PP(1)=PP(J)

341 FP(1)=FP(J)

342 IF (LBUB) THEN

343 PP(2)=PP(J)*C

344 ELSE

345 PP(2)=PP(J)/C

346 END IF 347 J=2

348 ELSE

349 PP(3)=PP(2)-FP(2)*(PP(2)-PP(1))/(FP(2)-FP(1))

350 IF ((PP(3).GT.MAX(PNEG,PPOS) .OR.

125
351 1 PP(3).LT.MIN(PNEG,PP0S)) .AND. LPNEG .AND. LPPOS)

352 1 PP(3)=PPOS-FPPOS*(PPOS-PNEG)/(FPPOS-FPNEG)

353 PP(1)=PP(2)

354 PP(2)=PP(3)

355 FP(1)=FP(2)

356 END IF 357 800 CONTINUE

358 WRITE (LUP.1040)

359 840 P=PP(J)

360 C

361 C ASSIGN RESULTS FOR PHASES 1 AND 2 TO LIQUID AND VAPOR PHASES

362 C DEPENDING ON WHETHER THE DEW OR BUBBLE POINT WAS CALCULATED.

363 C 364 IF (LBUB) THEN

365 XV=XX2(JJ)

366 VL=V1

367 VV=V2

368 ELSE

369 XL=XX2(JJ)

370 VL=V2

371 W=V1

372 END IF 373 C 374 C PRINT WARNING MESSAGES FOR ANY CASES OF NON-CONVERGENCE OCCURING

375 C ON FINAL CALL TO EACH ITERATION AND RETURN.

376 C 377 IF (LV1CON) WRITE (LUP.1050)

378 IF (LV2CON) WRITE (LUP.1055)

379 IF (LXCON) WRITE (LUP.1060)

380 RETURN

381 1010 FORMAT (1X, 'CRITICAL POINT OF PURE OR PSEUDO-PURE MATERIAL',

382 1

' EXCEEDED IN BUBLT')

383 1020 FORMAT (IX.'PURE MATERIAL PRESSURE ITERATION IN BUBLT'.

384 1 ' DID NOT CONVERGE')

385 1040 FORMAT (1X, 'MIXTURE PRESSURE ITERATION IN BUBLT DID NOT',

386 1 ' CONVERGE')

387 1050 FORMAT (1X, 'VOLUME ITERATION FOR PARENT PHASE DID',

388 1

' NOT CONVERGE')

389 1055 FORMAT (1X. 'VOLUME ITERATION FOR INCIPIENT PHASE DID',

390 1

' NOT CONVERGE')

391 1060 FORMAT (1X, 'COMPOSITION ITERATION IN BUBLT DID NOT CONVERGE')

392 END

126
1 FUNCTION ENTROP(T,V,X)

2 C

3 C COMPUTE SPECIFIC ENTROPY OF A SINGLE PHASE TWO-COMPONENT MIXTURE

4 C AS A FUNCTION OF TEMPERATURE, SPECIFIC VOLUME, AND COMPOSITION

5 C

6 C INPUTS:

7 C T - TEMPERATURE (K)

8 C V - SPECIFIC VOLUME (M**3/KMOL)

9 C X - COMPOSITION (MOLE FRACTION)

10 C

11 C OUTPUT:

12 C S - SPECIFIC ENTROPY (KJ/KMOL K)

13 C

14 C OTHER SUBROUTINES REFERENCED BY ENTROP:

15 C ESPAR - COMPUTATION OF EQUATION OF STATE PARAMETERS

16 C

17 C 18 IMPLICIT REAL (A-H.O-Z)

19 COMMON /ESPAR1/ AA, AB.BA.BB, F.C1 ,D1 ,C2 ,D2

20 COMMON /REF/ TREFA.TREFB.HRA.HRB.SRA.SRB.VRA.VRB

21 COMMON /CPDATA/ CPA0.CPA1 ,CPA2,CPB0,CPB1 ,CPB2

22 COMMON /HSPURE/ HPA,HPB,SPA,SPB,CP0A,CP0B

23 COMMON /RDATA4/ R 24 CALL ESPAR (1 .T.X.A.B)

25 B4=0.25*B

26 S=X*(SPA-SRA+R*LOG(V/VRA))+(1 .0-X)*(SPB-SRB+R*LOG(V/VRB))

27 1 +(C1*B-A»D1)/B**2*LOG((V+B)/V)+A*D1/B/(V+B)

28 1 -R*B4/(V-B4)**2*(4.0*V-3*B4)

29 1 -R*T»D1»0.5*V/(V-B4)**3*(2.0*V-B4)

30 IF (X.GT.0.0 .AND. X.LT.1.0) THEN

31 S=S-R*(X»LOG(X)+(1 .0-X)*LOG(1 .0-X))

32 END IF 33 ENTROP=S

34 RETURN

35 END

127
1 SUBROUTINE HCVCP(IQ,T ,V,X,H,CV,CP)

2 C

3 C GIVEN TEMPERATURE, MOLAR VOLUME AND COMPOSITION COMPUTE ENTHALPY

4 C AND/OR HEAT CAPACITY AT CONSTANT VOLUME AND/OR PRESSURE AS SPECIFIED

5 C BY OUTPUT QUALIFIER IQ. (SINGLE PHASE ONLY)

6 C

7 C INPUTS:

8 C IQ - OUTPUT QUALIFIER

9 C =1 COMPUTE ENTHALPY ONLY

10 C =2 ENTHALPY AND CONSTANT VOLUME HEAT CAPACITY

11 C =3 ENTHALPY AND HEAT CAPACITY AT CONSTANT VOLUME AND PRESSURE

12 C =4 COMPUTE HEAT CAPACITY AT CONSTANT VOLUME ONLY

13 C =5 HEAT CAPACITY AT CONSTANT VOLUME AND AT CONSTANT PRESSURE

14 C T - TEMPERATURE (K)

15 C V - MOLAR VOLUME (M**3/KM0L)

16 C X - COMPOSITION (MOLE FRACTION)

17 C

18 C OUTPUTS:

19 C H - MOLAR ENTHALPY (KJ/KMOL)

20 C CV - HEAT CAPACITY AT CONSTANT VOLUME (KJ/KMOL K)

21 C CP - HEAT CAPACITY AT CONSTANT PRESSURE (KJ/KMOL K)

22 C 23 C OTHER SUBROUTINES REFERENCED BY HCVCP:

24 C ESPAR - COMPUTATION OF EQUATION OF STATE PARAMETERS

25 C 26 C 27 IMPLICIT REAL (A-H.O-Z)

28 COMMON /ESPAR1/ AA.AB.BA.BB, F.C1 ,D1 ,C2,D2

29 COMMON /REF/ TREFA.TREFB.HRA.HRB.SRA.SRB.VRA.VRB

30 COMMON /HSPURE/ HPA.HPB,SPA.SPB ( CP0A,CP0B

31 COMMON /RDATA4/ R 32 CALL ESPAR (IQ.T.X.A.B)

33 B4=0.25^B

34 VB=V+B

35 VBL=LOG(V/VB)

36 VB4=V-B4

37 VB43=VB4**3

38 RT=R*T

39 IF (IQ.LE.3) THEN

40 C COMPUTE MOLAR ENTHALPY AS A FUNCTION OF T, V, X

41 H=X*(HPA-HRA)+(1 . 0-X) * (HPB-HRB)

42 1 +((A+(A«D1/B-C1)*T)*VBL+A*(D1*T-B)/VB)/B

43 1 +2.0*RT»V»(2.0»V-B4)*(B4-0.25»D1»T)/VB43

44 END IF 45 IF (IQ.GE.2) THEN

46 C COMPUTE CONSTANT VOLUME MOLAR HEAT CAPACITY

47 D12=D1*D1

48 CV=X*(CP0A-R)+(1 .0-X)*(CP0B-R)

49 1 +(R*V*((0.375*D12*T/VB4+0.5*D2*T+D1)*(B4-2.0»V)

50 1 +0.125*D12*T)/VB43

51 1 +((1.0/VB+VBL/B)*(A»D2*B+2.0*(C1*D1*B-A*D12))/B

52 1 -C2*VBL-A*D12/VB»*2)/B)*T

53 IF (IQ.EQ.3 .OR. IQ.EQ.5) THEN

54 C COMPUTE MOLAR HEAT CAPACITY AT CONSTANT PRESSURE USING CV

55 Y=B4/V

56 DPDT=2.0*R/VB4*(-1 .0+(-0.25*T*D1+(V*V»(1 .0+0.75*T*D1/VB4))

57 1 /VB4)/VB4)+(R+(-C1+A*D1/VB)/VB)/V

58 DPDV=(-RT*(1 .0+(4.0+(4.0+(-4.0+Y)*Y)*Y)*Y)/(1 .0-Y)**4

59 1 +A*(2.0*V+B)/VB**2)/V**2

60 CP=CV-DPDT*DPDT*T/DPDV

61 END IF 62 END IF 63 RETURN

64 END

128
1 BLOCK DATA BDESC

2 C

3 C THIS ROUTINE INITIALIZES THE COMMON BLOCKS CONTAINING INFORMATION

4 C ABOUT THE PURE COMPONENTS. IT IS NOT REFERENCED DIRECTLY BY ANY

5 C OTHER SUBROUTINE BUT MUST BE INCLUDED IN THE EXECUTABLE ELEMENT.

6 C DATA ARRAYS ARE DIMENSIONED TO ACCOMODATE ADDITIONAL

7 C PURE COMPONENTS.

8 C

9 C EXPLANATION OF CONSTANTS:

10 C COEFF(I.J) - FOR REFRIGERANT J, COEFFICIENTS OF A, B, CP0

11 C CURVE FITS:

12 C A = A0 * EXP(A1»T + A2*T*T) (KJ M**3/KMOL**2)

13 C B = B0 + B1*T + B2*T»T (M**3/KMOL)

14 C CP0 = C0 + C1»T + C2*T»T (KJ/KMOL K)

15 C (STORED IN ORDER A0, A1 , A2 ,B0,B1 ,B2 ,C0,C1 ,C2)

16 C CRIT(I.J) - FOLLOWING INFORMATION FOR REFRIGERANT J: 17 C 1=1- MOLECULAR WEIGHT

18 C 2 - REFERENCE TEMPERATURE FOR ENTHALPY AND ENTROPY (K) 19 C 3 - CRITICAL TEMPERATURE (K)

20 C 4 - CRITICAL PRESSURE (KPA)

21 C 5 - CRITICAL VOLUME (M**3/KMOL)

22 C HREF(J) - REFRIGERANT NAME (ASHRAE DESIGNATION)

23 C HZERO(J) - VALUE OF SATURATED LIQUID ENTHALPY OF REFRIGERANT

24 C J AT ITS REFERENCE TEMPERATURE (KJ/KMOL)

25 C SZERO(J) - VALUE OF SATURATED LIQUID ENTROPY AT REFERENCE

26 C TEMPERATURE (KJ/KMOL K)

27 C R - GAS CONSTANT (KJ/KMOL K)

28 C TOLR - RELATIVE CONVERGENCE TOLERANCE FOR ITERATION LOOPS

29 C SHOULD BE AT LEAST 10 TIMES LARGER THAN MACHINE PRECISION

30 C ITMAX - MAXIMUM ITERATION COUNT FOR ITERATIVE LOOPS

31 C LUP - LOGICAL UNIT TO WHICH ANY WARNING MESSAGES ARE WRITTEN

32 C 33 IMPLICIT REAL (A-H.O-Z)

34 DIMENSION COEFF(9,20) ,CRIT(5,20) ,HZERO(20) ,SZERO(20)

35 CHARACTERS HREF(20)

36 COMMON /ESDATA/ COEFF.CRIT

37 COMMON /HREF1/ HREF

38 COMMON /HSZERO/ HZERO.SZERO

39 COMMON /RDATA4/ R 40 COMMON /TOL/ TOLR, ITMAX, LUP

41 DATA R /8.314/

42 DATA TOLR /1 .0E-7/

43 DATA ITMAX. LUP /20,6/

44 C 45 C DATA FOR R11, R12, R13, R13B1 , R14, R22 , R23, R113, R114,

46 C R142B. R152A FOLLOW.

47 C 48 C

49 C R11, TRICHLOROFLUOROMETHANE

50 C

51 DATA HREF(1) /'R11 */

52 DATA (CRIT(I,1),I=1 ,5) /137 . 37 ,233. 15,471 . 2 ,4467. ,0.247/

53 DATA HZERO(1) ,SZERO(1) /0. 0,0.0/

54 DATA (COEFF(I , 1 ) , 1=1 ,9) /4971 .54,-2. 24669 E-3, -5. 1 1943E-7

,

55 1 0.176659,-1 .74531 E-4, -3. 4971 7E-8,

56 1 22. 041 8, 0.260895, -2. 4531 9E-4/

57 C

58 C R12, DICHLORODIFLUOROMETHANE

59 C 60 DATA HREF(2) /'R12 1 /

61 DATA (CRIT(I ,2), 1=1 ,5) /120. 91 ,233. 15,384.95,4180. ,0.217/

62 DATA HZERO(2),SZERO(2) /0. 0,0.0/

63 DATA (COEFF(I,2),I=1,9) /3524. 12,-2. 77230E-3. -6.731 80E-7,

64 1 0.153755,-1 .84195E-4, -5. 03644E-8,

65 1 17.5387,0.248546,-2. 16271 E-4/

66 C

67 C R13, CHLOROTRIFLUOROMETHANE

68 C

69 DATA HREF(3) /'R137

70 DATA (CRIT(I,3),I=1 ,5) /104. 46, 233. 15,302.0.3921 . .0. 181/

129
71 DATA HZER0(3),SZER0(3) /0. 0,0.0/

72 DATA (COEFF(I,3),I=1.9) /2298. 13,-3. 41828E-3.-1 .52430E-6.

73 1 0.128141 ,-1 .84474E-4.-1 .07951E-7,

74 1 13.9300,0.232181 ,-1 .82929E-4/

75 C 76 C R13B1, BROMOTRIFLUOROMETHANE

77 C

78 DATA HREF(4)/'R13B1

V

79 DATA (CRIT(I,4),I=1 ,5) /148. 91 ,233. 15,340.2,4017. ,0.200/

80 DATA HZER0(4),SZER0(4) /0. 0,0.0/

81 DATA (COEFF(I,4),I=1 ,9) /2728. 10,-2.79791 E-3,-1 .50848E-6,

82 1 0.139949,-1 .82428E-4, -7. 75898E-8,

83 1 19.9537,0.216394,-1 .70241E-4/

84 C 85 C R14, TETRAFLUOROMETHANE

86 C 87 DATA HREF(5) /'R14'/

88 DATA (CRIT(I,5),I=1 ,5) /88. 00, 200. 00, 227. 5, 3795. ,0. 141/

89 DATA HZERO(5),SZERO(5) /0. 0,0.0/

90 DATA (C0EFF(I,5),I=1 .9) /1393.60.-4.81985E-3.-1 .89167E-6,

91 1 0.100601 ,-1 .94974E-4.-1 .35408E-7,

92 1 11.0629.0.209740,-1 .40992E-4/

93 C 94 C R22, CHLORODIFLUOROMETHANE

95 C 96 DATA HREF(6) /

, R22'/

97 DATA (CRIT(I,6),I=1,5) /86. 47 ,233. 15,369. 3,5054. ,0. 169/

98 DATA HZERO(6),SZERO(6) /0. 0,0.0/

99 DATA (COEFF(I,6),I=1,9) /2514. 59.-2 .38706E-3.-1 .83653E-6,

100 1 0.113681,-1 .16201E-4, -9. 24562E-8,

101 1 17.0547.0. 161633.-9. 12559E-5/

102 C 103 C R23, TRIFLUOROMETHANE

104 C 105 DATA HREF(7) /'R23'/

106 DATA (CRIT(I,7).I=1 ,5) /70. 01 .233. 15.299. 1 .4900. ,0. 133/

107 DATA HZERO(7),SZERO(7) /0. 0,0.0/

108 DATA (COEFF(I,7).I=1,9) /2025.93.-4.68206E-3.9.95524E-7,

109 1 0.103137.-2.29653E-4.1.55760E-7.

110 1 20. 4760,0. 106183. -1.21892E-5/

111 C

112 C R113, 1,1 ,2-TRICHLOROTRIFLUOROETHANE

113 C 114 DATA HREF(8) /'R113'/

115 DATA (CRIT(I.8),I=1 .5) /187. 38, 233. 15.487.5,3456. .0.329/

116 DATA HZERO(8),SZERO(8) /0. 0,0.0/

117 DATA (COEFF(I,8),I=1 ,9) /7332.59 ,-2. 20396E-3.-7 .26656E-7,

118 1 0.230713,-1 .87956E-4.-1 .0611 4E-7,

119 1 76.2637,0.119641 ,7.1 8786E-5/

120 C

121 C R114, 1,2-DICHLOROTETRAFLUOROETHANE

122 C 123 DATA HREF(9) /'R114'/

124 DATA (CRIT(I,9),I=1 ,5) /170. 92, 233. 15,419.03,3304. ,0.307/

125 DATA HZERO(9),SZERO(9) /0. 0,0.0/

126 DATA (COEFF(I,9),I=1 ,9) /9771 . 35.-5. 85557E-3, 3. 9941 3E-6,

127 1 0.306318. -7. 96444E-4. 7. 81059E-7,

128 1 20.7005,0.464035.-4.1 7589 E-4/

129 C 130 C R142B. 1-CHLORO-1 , 1-DIFLUOROETHANE

131 C

132 DATA HREF(10) /'R142B'/

133 DATA (CRIT(I.10),I=1 ,5) /100. 49 .233. 15,410.3,4120. .0. 231/

134 DATA HZERO(10),SZERO(10) /0. 0,0.0/

135 DATA (COEFF(I,10),I=1 ,9) /2990.00,-5.40563E-4,-4. 12642E-6,

136 1 0. 146006. -8. 92503E-5, -1 .80562E-7,

137 1 23.7611.0.231706.-1 .06534E-4/

138 C

139 C R152A, 1 , 1-DIFLUOROETHANE

140 C

130
141 DATA HREF(11) /'R152A'/

142 DATA (CRIT(I . 11 ) , 1=1 ,5) /66. 05,233. 15,386. 7, 4492. ,0. 181/

143 DATA H2ER0(11) ,SZER0(11) /0. 0,0.0/

U4 DATA (C0EFF(I,11),I=1,9) /2254.37.-5.87778E-4.-4.37432E-6,

145 1 0.116521 ,-9. 04883E-5.-1 .14563E-7,

146 1 22.2804,0. 154009, -3. 06670E-6/

1 47 END

131
1 SUBROUTINE ESPAR (IQ.T.X.A.B)

2 C

3 C THIS ROUTINE CALCULATES THE EQUATION OF STATE PARAMETERS AND THEIR

4 C TEMPERATURE DEVI VAT IVES AS A FUNCTION OF TEMPERATURE AND COMPOSITION

5 C AS NEEDED BY THE OTHER PROPERTY ROUTINES. BASED ON THE VALUE OF THE

6 C INPUT QUALIFIER THE NECESSARY PARAMETERS ARE CALCULATED EXCEPT THAT

7 C IF THE TEMPERATURE AND COMPOSITION ARE UNCHANGED FROM THE LAST CALL

8 C THE PREVIOUS VALUES ARE USED. THE TEMPERATURE DEPENDENCE OF THE

9 C A, B. AND CP0 PARAMETERS ARE CONTAINED ENTIRELY WITHIN ESPAR AND

10 C THE STATEMENT FUNCTIONS (IN BUBLT) FOR GIBBS FREE ENERGY

11 C AND CHEMICAL POTENTIAL; ALTERNATE EXPRESSIONS FOR 'A' AND 'B' 12 C REQUIRE CHANGING ONLY THESE ROUTINES.

13 C

14 C INPUTS:

15 C IQ - INPUT QUALIFIER

16 C =0 COMPUTE ONLY A AND B

17 C >= 1 ALSO COMPUTE TEMPERATURE DERIVATIVES OF A AND B

18 C >= 2 ALSO COMPUTE SECOND DERIVATIVE OF A AND B AND

19 C IDEAL GAS HEAT CAPACITY

20 C = 1 , 2 OR 3 ALSO COMPUTE CONSTANTS FOR PURE COMPONENT ENTHALPY

21 C AND ENTROPY

22 C T - TEMPERATURE (K)

23 C X - COMPOSITION (MOLE FRACTION COMPONENT A)

24 C

25 C OUTPUTS:

26 C A - 'A' PARAMETER FOR MIXTURE AT T, X 27 C B - 'B

1 PARAMETER FOR MIXTURE AT T, X 28 C

29 C OUTPUTS (VIA COMMON BLOCKS):

30 C AA - 'A* PARAMETER FOR PURE COMPONENT A

31 C AB - *A' PARAMETER FOR PURE COMPONENT B 32 C BA - 'B' PARAMETER FOR PURE COMPONENT A

33 C BB - 'B' PARAMETER FOR PURE COMPONENT B 34 C F - MIXTURE INTERACTION PARAMETER

35 C DADT - TEMPERATURE DERIVATIVE OF A

36 C DBDT - TEMPERATURE DERIVATIVE OF B 37 C D2ADT2 - SECOND DERIVATIVE OF A WITH RESPECT TO TEMPERATURE

38 C D2BDT2 - SECOND DERIVATIVE OF B WITH RESPECT TO TEMPERATURE

39 C HPA - INTEGRAL OF CP0 WITH RESPECT TO TEMP FOR PURE A

40 C HPB - INTEGRAL OF CP0 WITH RESPECT TO TEMP FOR PURE B

41 C SPA - INTEGRAL OF (CP0 - R)/T WITH RESPECT TO TEMP FOR PURE A

42 C SPB - INTEGRAL OF (CP0 - R)/T WITH RESPECT TO TEMP FOR PURE B 43 C CP0A - PERFECT GAS HEAT CAPACITY FOR COMPONENT A (KJ/KMOL K)

44 C CP0B - PERFECT GAS HEAT CAPACITY FOR COMPONENT B (KJ/KMOL K)

45 C

46 C

47 IMPLICIT REAL (A-H.O-Z)

48 COMMON /ESPAR 1/ AA.AB.BA.BB, F, DADT , DBDT, D2ADT2.D2BDT2

49 COMMON /RDATA1/ AA0.AA1 , AA2.AB0.AB1 .AB2.BA0.BA1 ,BA2,

50 1 BB0 BB1 BB2 F0 F1 5

COMMON /CPDATA/ CPA0 , CPA 1 , CPA2 , CPB0 , CPB 1 , CPB2

52 COMMON /HSPURE/ HPA,HPB.SPA.SPB,CP0A,CP0B

53 COMMON /REF/ TREFA,TREFB,HRA,HRB,SRA,SRB,VRA,VRB

54 COMMON /RDATA4/ R 55 SAVE TLAST0.TLAST1 .TLAST2 ,TLAST3 , XLAST0.XLAST1 .XLAST2 , ALAST

,

56 1 B LAST , SQAB. XB, XB2, XBX.DAA, DAB. DBA, DBB.DSQAB.F0 LAST, F1 LAST

57 DATA TLAST0.TLAST1 .TLAST2 .TLAST3 , XLAST0.XLAST1 .XLAST2 /7*-999./

58 IF (IQ.LT.0) THEN

59 IQ=ABS(IQ)

60 GOTO 100

61 END IF 62 A=ALAST

63 B=BLAST

64 IF (T.NE.TLAST0) GOTO 100

65 IF (F0.NE.F0LAST .OR. F1 .NE. F1 LAST) GOTO 110

66 IF (X.NE.XLAST0) GOTO 120

67 IF (IQ.LE.0) RETURN

68 IF (T.NE.TLAST1 .OR. X.NE.XLAST1) GOTO 200

69 IF (T.NE.TLAST3) GOTO 230

70 IF (IQ.LE.1) RETURN

132
71 IF (T.NE.TLAST2 .OR. X.NE.XLAST2) GOTO 300

72 RETURN

73 100 AA=AA0*EXP((AA1+AA2*T)*T)

74 AB=AB0*EXP((AB1+AB2*T)*T)

75 BA=BA0+(BA1+BA2*T)*T

76 BB=BB0+(BB1+BB2»T)*T

77 TLAST0=T

78 110 F=F0+F1*T

79 F0LAST=F0

80 F1 LAST=F1

81 F=F0+F1»T

82 120 X2=X«X

83 XB=1.0-X

84 XB2=XB*XB

85 XBX=X*XB

86 SQAB=SQRT(AA*AB)

87 A=X2»AA+2.0*XBX*(1 .0-F)*SQAB+XB2*AB

88 B=X*BA+XB*BB

89 ALAST=A

90 BLAST=B

91 XLAST0=X

92 IF (IQ.LE.0) RETURN

93 200 DAA=AA*(AA1+2.0*AA2»T)

94 DAB=AB*(AB1+2.0*AB2»T)

95 DBA=BA1+2.0»BA2*T

96 DBB=BB1+2.0*BB2*T

97 DSQAB=0 . 5 ( AA * DAB+AB * DAA ) /SQAB

98 DADT=X2»DAA+XB2»DAB+2.0»XBX»((1 .0-F)*DSQAB-F1 *SQAB)

99 DBDT=X*DBA+XB*DBB

100 TLAST1=T

101 XLAST1=X

102 230 IF (IQ.LE.3) THEN

103 HPA=(CPA0+(0.5*CPA1+CPA2/3.0*T)*T)*T

104 HPB=(CPB0+(0.5»CPB1+CPB2/3.0*T)»T)*T

1 05 SPA=(CPA0-R) * L0G(T/TREFA)+CPA1 • (T-TREFA)+0 . 5*CPA2*

106 1 (T*T-TREFA**2)

1 07 SPB=(CPB0-R) * L0G(T/TREFB)+CPB1 * (T-TREFB)+0 . 5*CPB2*

108 1 (T*T-TREFB**2)

109 TLAST3=T

110 END IF

111 IF (IQ.LE.1) RETURN

112 300 D2AA=AA*((AA1+2.0*AA2*T)**2+2.0*AA2)

113 D2AB=AB*((AB1+2.0»AB2»T)»»2+2.0*AB2)

114 D2SQAB=(0.5»(AA*D2AB+AB*D2AA)+DAA*DAB-DSQAB*DSQAB)/SQAB

1 1

5

D2ADT2=X2*D2AA+XB2*D2AB+2 . 0»XBX» (

( 1 . 0-F) «D2SQAB-2 . 0* F1 *DSOAB)

116 D2BDT2=2.0*(X*BA2+XB*BB2)

117 CP0A=CPA0+(CPA1+CPA2*T)*T

118 CP0B=CPB0+(CPB1+CPB2*T)»T

119 TLAST2=T

120 XLAST2=X

121 RETURN

122 END

133
1 SUBROUTINE PLIMIT (T, A.B.VL, VU.PLOW.PUP)

2 C

3 C GIVEN TEMPERATURE AND EQUATION OF STATE PARAMETERS, THIS

4 C ROUTINE CALCULATES THE UPPER AND LOWER BOUNDS ON PRESSURE

5 C FOR WHICH THERE ARE BOTH LIQUID AND VAPOR SOLUTIONS TO THE

6 C EQUATION OF STATE. IT CARRIES OUT TWO BISECTION METHOD

7 C ITERATIONS TO FIND THE POINTS WHERE THE DERIVATIVE OF PRESSURE

8 C W.R.T. VOLUME IS ZERO.

9 C

10 C INPUTS:

11 C T - TEMPERATURE (K) 12 C A,B - EQUATION OF STATE PARAMETERS AT TEMPERATURE T

13 C

14 C OUTPUTS:

15 C PLOW - LOWER BOUND ON PRESSURE (PLOW CAN BE NEGATIVE, THE

16 C CALLING PROGRAM MUST CHECK AND CORRECT FOR NEGATIVE

17 C PRESSURES)

18 C PUP - UPPER BOUND ON PRESSURE (KPA)

19 C VL - MOLAR VOLUME AT PLOW (M*»3/KMOL)

20 C VU - MOLAR VOLUME AT PUP (M**3/KMOL)

21 C 22 C OTHER SUBROUTINES REFERENCED:

23 C NONE

24 C

25 C 26 IMPLICIT REAL (A-H.O-Z)

27 COMMON /RDATA4/ R 28 COMMON /TOL/ TOLR, ITMAX, LUP

29 c 30 C STATEMENT FUNCTIONS FOR THE EVALUATION OF PRESSURE AS A

31 C FUNCTION OF V AND THE DERIVATIVE OF PRESSURE W.R.T

32 C VOLUME AS A FUNCTION OF V

33 C

34 P(RT,V,Y,A,B)=(RT»(1 .0+(1 .0+(1 .0-Y)*Y)*Y)/(1 .0-Y)**3

35 1 -A/(V+B))/V

36 DP(RT , V , A , B , B4 , B42W-RT* (B42*B42+(-4 . 0»B42*B4+(4 . 0»B42

37 1 +(4.0«B4+V)*V)*V)*V)/(V-B4)**4+A*(2.0*V+B)/(V+B)»*2)/V**2

38 C 39 B4=0.25*B

40 B42=B4*B4

41 RT=R*T

42 C

43 C STARTING AT A VOLUME OF 12.0*B4 (WHICH HAS A POSITIVE SLOPE

44 C FOR ALL 'REASONABLE' VALUES OF A, B, T) REDUCE THE VOLUME

45 C UNTIL A NEGATIVE SLOPE OF P W.R.T. V IS FOUND AND THEN BEGIN

46 C BISECTION METHOD TO FIND LOWER BOUND ON VOLUME AND PRESSURE.

47 C

48 VC=12.0272727*B4

49 V=VC

50 DO 100 IT=1.ITMAX

51 DPDV=DP(RT,V ^,6,64,642)

52 IF (DPDV.LE.0.0) GOTO 116

53 VPOS=V

54 V=0.5*(V+B4)

55 100 CONTINUE

56 116 VNEG=V

57 DO 120 IT=1,20

58 V L=0 . 5 * ( VN EG+VPOS

)

59 DPDV=DP(RT.VL.A,B.B4,B42)

60 IF (DPDV.LT.0.0) THEN

61 VNEG=VL

62 ELSE

63 VPOS=VL

64 END IF 65 120 CONTINUE

66 Y=B4/VL

67 PLOW=P(RT,VL,Y.A.B)

68 C

69 C STARTING AT V = 2*A/RT INCREASE V UNTIL A NEGATIVE

70 C SLOPE IS FOUND; USE WITH V = 12.0*B TO BEGIN BISECTION

134
71 C ITERATION FOR UPPER BOUND

72 C

73 VPOS=VC

74 V=2.0*A/RT

75 DO 160 IT=1 , ITMAX 76 DPDV=DP(RT,V,A,B.B4,B42)

77 IF (DPDV.LE.0.0) GOTO 164

78 VPOS=V

79 V=2.0*V

80 160 CONTINUE

81 164 VNEG=V

82 DO 180 IT=1 ,20

83 VU=0.5*(VNEG+VPOS)

84 DPDV=DP(RT,VU,A,B,B4 > B42)

85 IF (DPDV.LT.0.0) THEN

86 VNEG=VU

87 ELSE

88 VPOS=VU

89 END IF 90 180 CONTINUE

91 Y=B4/VU

92 PUP=P(RT,VU,Y,A,B)

93 RETURN

94 END

135
1 SUBROUTINE VIT (T.P.A.B.VS. LLIQI . LVCON)

2 C

3 C GIVEN TEMPERATURE, PRESSURE, AND EQUATION OF STATE

4 C PARAMETERS, THIS ROUTINE CALCULATES THE LIQUID OR VAPOR

5 C MOLAR VOLUME THAT SATISFIES THE EQUATION OF STATE.

6 C

7 C INPUTS:

8 C T - TEMPERATURE (K)

9 C P - PRESSURE (KPA)

10 C A.B - EQUATION OF STATE PARAMETERS AT TEMPERATURE T

11 C VS - INITIAL GUESS FOR VOLUME. IN ABSENCE OF BETTER

12 C GUESSES SUGGESTED VALUES ARE:

13 C LIQUID: VS=0.8*B

14 C VAPOR: VS=R*T/P

15 C LLIQI - LOGICAL VARIABLE

16 C IF LLIQI = .TRUE. COMPUTE LIQUID VOLUME

17 C IF LLIQI = .FALSE. COMPUTE VAPOR VOLUME

18 C NOTE: IF EITHER THE TEMPERATURE OR THE PRESSURE IS ABOVE

19 C THE CRITICAL VALUE, ONLY ONE SOLUTION EXISTS AND THE

20 C VALUE OF LLIQI HAS NO EFFECT.

21 C

22 C OUTPUTS:

23 C VS - MOLAR VOLUME (M**3/KG MOL)

24 C LVCON - ERROR FLAG; IF LVCON = .TRUE. THE ITERATION HAS

25 C NOT CONVERGED

26 C

27 C OTHER SUBROUTINES REFERENCED:

28 C NONE

29 C

30 C (FOR EXPLANATION OF NOMENCLATURE SEE BUBLT)

31 C

32 C NOTE: THIS ROUTINE IS WRITTEN IN DOUBLE PRECISION EXCEPT

33 C THAT THE ARGUMENTS ARE SINGLE PRECISION

34 C

35 IMPLICIT DOUBLE PRECISION (A-H.O-Z)

36 LOGICAL LLIQ. LVCON, LLIQI

37 REAL T,P,A,B,R,VS,TOLR,TC,PC

38 COMMON /RDATA4/ R 39 COMMON /TOL/ TOLR, ITMAX. LUP 40 LVCON=. FALSE.

41 LLIQ=LLIQI

42 V=VS

43 VL=LOG(V)

44 PL=LOG(P)

45 RT=R»T

46 B4=0.25*B

47 B4L=LOG(B4)

48 IF (VL.LT.B4L) VL=B4L+0.5

49 TC=A/(B*4.398909*R)

50 PC=0.02386944*A/B»*2

51 VCL=LOG(12.0272727*B4)

52 IF (P.GT.PC) THEN

53 LLIQ=.TRUE.

54 ELSE IF (T.GT.TC) THEN

55 LLIQ=. FALSE.

56 END IF 57 C

58 C ENTER NEWTONS METHOD ITERATION FOR VOLUME. FOR LIQUIDS

59 C (OR FLUIDS ABOVE THE CRITICAL PRESSURE) THE ITERATION

60 C IS CARRIED OUT IN TRANSFORMED COORDINATES OF LOG (V). FOR

61 C VAPOR (OR FLUIDS AT SUPERCRITICAL TEMPERATURES BUT PRESSURES

62 C BELOW THE CRITICAL VALUE) THE ITERATION IS IN TERMS OF

63 C LOG (V) AND LOG (P) . THE ITERATION HAS CONVERGED WHEN

64 C THE PRESSURE CALCULATED FROM THE EQUATION OF STATE AGREES

65 C WITH THE INPUT PRESSURE.

66 C

67 DO 100 IT=1 .ITMAX

68 IF (VL.GT.VCL .EQV. LLIQ) VL=VCL

69 VLS=VL

70 Y=B4/V

136
71 C 72 C CALCULATE PRESSURE AS A FUNCTION OF VOLUME AND THE

73 C DERIVATIVE OF THE PRESSURE W.R.T. LOG (VOLUME).

74 C 75 P2=(RT*(1 .0+(1 .0+(1 .0-Y)»Y)*Y)/(1 .0-Y)**3-A/VB)/V

76 DPDLV=RT/V*(-1 .0+(-4.0+(-4.0+(4.0-Y)*Y)*Y)*Y)/(1 .0-Y)»»4

77 1 +A»(2.0*V+B)/(V*VB*VB)

78 IF (LLIQ) THEN

79 IF (DPDLV.GE.0.0) THEN

80 VL=0.5*(B4L+VLS)

81 ELSE

82 FVDP=(P2-P)/DPDLV

83 IF (ABS(FVDP/P).LT.0.001*TOLR) THEN

84 VS=EXP(VL-FVDP)

85 RETURN

86 ELSE

87 VL=VL-FVDP

88 IF (VL.LE.B4L) VL=0.5*(B4L+VLS)

89 END IF 90 END IF

91 ELSE

92 IF (DPDLV.GE.0.0 .OR. P2.LE.0.0) THEN

93 VL=VL+0 .

5

94 ELSE

95 FVDPL=(LOG(P2)-PL)*P2/DPDLV

96 IF (ABS(FVDPL).LT.0.001*TOLR) THEN

97 VS=EXP(VL)

98 RETURN

99 END IF 100 VL=VL-FVDPL

101 IF (ABS(VL-VLS).GT.1.5) VL=VLS+SIGN(1 .0D0.VL-VLS)

102 IF (VL.LT.VCL) VL=0.5»(VLS+VCL)

103 END IF 104 END IF 105 V=EXP(VL)

106 100 CONTINUE

107 LVCON=.TRUE.

1 08 VS=V

109 RETURN

110 END

137
Exampi e Run

The following example main program uses the property routines to calculate

tables of saturation properties for a specified refrigerant mixture. Compiled

versions of the following routines must be linked together before execution:

MAIN PROGRAM GRID

BCONST

BUBLT

ENTROP

HCVCP

BDESC

ESPAR

PLIMIT

VIT

The program reads the following information from data 'cards': (values used

for the example are shown in parenthesis):

IR1,IR2 - code numbers for components of mixture (4,11)

F - mixture interaction parameter (0.0902)

TB,TE,TDELT - beginning and ending temperature and temperature interval

for which to calculate properties (260., 340., 40.)

XB,XE,XDELT - beginning and ending composition and interval for which to

calculate properties (0.0, 1.0, 0.1).

A listing of the main program and sample output follows.

138
1 PROGRAM GRID

2 C THIS PROGRAM SERVES TO TEST THE PROPERTY ROUTINES BY CALCULATING

3 C TABLES OF DEW AND BUBBLE POINT PROPERTIES FOR A MIXTURE

4 C

5 C INPUT INFORMATION:

6 C IR1.IR2 - CODE NUMBERS FOR COMPONENTS OF MIXTURE

7 C F0 - MIXTURE INTERACTION PARAMETER

8 C TB.TE.DELT - BEGINNING 4 ENDING TEMPERATURE AND TEMPERATURE

9 C INTERVAL FOR WHICH TO CALCULATE PROPERTIES

10 C XB.XE.DELX - BEGINNING k ENDING COMPOSITION AND COMPOSITION

11 C INTERVAL FOR WHICH TO CALCULATE PROPERTIES

12 C

13 C OUTPUT INFORMATION (CALCULATED FOR EVEN INCREMENTS OF LIQUID

14 C COMPOSITION)

15 C XL - LIQUID COMPOSITION (MOLE FRACTION)

16 C XV - VAPOR COMPOSITION (MOLE FRACTION) IN EQUILIBRIUM WITH XL

17 C P - SATURATION PRESSURE (KPA)

18 C VL - SATURATED LIQUID VOLUME (M**3/KMOL) AT COMPOSITION XL

19 C VV - SATURATED VAPOR VOLUME (M**3/KMOL) AT COMPOSITION XV

20 C HL - LIQUID ENTHALPY (KJ/KMOL) AT VL. XL

21 C HV - VAPOR ENTHALPY (KJ/KMOL) AT W, XV

22 C SL - LIQUID ENTROPY (KJ/KMOL K) AT VL, XL

23 C SV - VAPOR ENTROPY (KJ/KMOL K) AT W. XV

24 C CVL - CONST VOLUME HEAT CAPACITY (KJ/KMOL K) FOR LIQUID AT VL. XL

25 C CVV - CONST VOLUME HEAT CAPACITY (KJ/KMOL K) FOR VAPOR AT W. XV

26 C CPL - CONST PRESSURE HEAT CAPACITY (KJ/KMOL K) FOR LIQUID

27 C CPV - CONST PRESSURE HEAT CAPACITY (KJ/KMOL K) FOR VAPOR

28 C 29 C OTHER SUBROUTINES REFERENCED:

30 C ENTIRE SET OF PROPERTY ROUTINES EXCEPT FOR FITAB. FITF

31 C

32 DIMENSION COEFF(9 .20) ,CRIT(5.20)

33 CHARACTER*6 HREF(20)

34 LOGICAL LCRIT

35 COMMON /ESDATA/ COEFF.CRIT

36 COMMON /HREF1/ HREF

37 COMMON /RDATA1/ A.B.F0.F1

38 COMMON /TOL/ TOLR. ITMAX. LUP 39 DIMENSION A(3.2) ,B(3.2)

40 OPEN (UNIT=LUP,FILE='OUTPUT')

41 READ (.) IR1 .IR2.F0.TB.TE.DELT.XB.XE.DELX

42 TC1=CRIT(3,IR1)

43 TC2=CRIT(3,IR2)

44 IF (TC1.GT.TC2) THEN

45 I1=IR2

46 I2=IR1

47 ELSE

48 I1=IR1

49 I2=IR2

50 END IF

51 CALL BCONST (1 1 . 12 , F0.0.0)

52 DO 400 T=TB,TE,DELT

53 WRITE (LUP.1000) T ,HREF( 1

1 ) ,HREF( 12) . F0

54 DO 380 XL=XB,XE+TOLR.DELX

55 CALL BUBLT (T , XL , XV, P.VL.W, . TRUE. . LCRIT)

56 IF (LCRIT) THEN

57 LCRIT=.FALSE.

58 WRITE (LUP, 1040)

59 GOTO 400

60 END IF

61 CALL HCVCP (3,T,VL,XL,HL,CVL,CPL)

62 CALL HCVCP (3.T , VV.XV.HV.CVV.CPV)

63 SL=ENTROP(T,VL,XL)

64 SV=ENTROP(T,VV,XV)

65 WRITE (LUP. 1020) XL, XV. P.VL.W

66 1 .HL.HV.SL.SV.CVL.CVV.CPL.CPV

67 380 CONTINUE

68 400 CONTINUE

69 STOP

70 1000 FORMAT (' 1 \ 'DEW/BUBBLE LINES AT T ='.F6.1.' K'/

139
COMPONENT B: ' ,A<5/ F =\F7.4//

P VL W

SL SV CVL

71 1 1X, "COMPONENT A: ' ,A6.

72 1 1X, 'MIXING COEFFICIENT,

73 1 1X. ' XL XV

74 1 ,

' HL HV SL SV CVL CW

75 1 CPL CPV/

76 1 1X,' (MOL FRAC A) (KPA) (L/MOL)

77 1 '(KJ/KG MOL) '.16('-'),' (KJ/KG MOL K) '.IBC-'))

78 1020 FORMAT (1X.2F9 .4, F9.2, F9.5, F9 .3.2F9 . 1 ,6F9.3)

79 1040 FORMAT (/1X. 'PSEUDO-PURE COMPONENT CRITICAL POINT EXCEEDED')

80 END

140
>O

COCDU,>CM<-^-<Oi-'*mir>

o>r-cMoooomf>-oi>--<i-co f*--*o>ooin--(OT-r^(Mr^

r<-)tOI^OOO>®®'-'>-CMCM

<oco<oto<or>r-*i^i^r^r^

IDNOi'-NIO'-IDSIO'-

r^*-rorO(OC0O>O5<OfOO> "-^••-r^roOTmo^-tom

(DOOMOKJiWOiOSGiG)

(OOlOID'-OlOlONincS >^*>-in*ifi>-iDiomo)0 >^COrOCO'-rOincOCO®CMCO

o

Qinminininininmininm

00rO00'-r~O>»OO5©«DrO

iDtomif)'-(D'-Tj-t^moo

•*oO'-'<i-r^a)CMT)-«DCDo> r^r^oooooooooi05CDO>05

<7><NCO<M»-0O®'*J-«-)l-*m

©®©©»-*-i-CMCMCMCM

in<oooc\iK)-*-<N-«*-'r-<0'«*-

-—vK) roa>®CM®l^^-tOrOif>

i^rorO0Ji/">0>iOC0®-tfC0CM

_J®^T-{N|(NrO(0'*'*'*m

0<0<£(OCO(0(0(OCOCO(OtO

>o

a.o

•^csit~~o)fM©ooi^o>oom cooooorofo<NOiO'-tDr^

cor^oo<No>"-G>i^T-r-»if) o>®^»O'i-r^®f")C0'-cM i-t-T-r-r-NMNIOlO

cMoo-*tDr^c^to©®in

<*-T-ooc»ooooioincMmoo

0)r^to<or^cM(MO)oo»-r^

cMr>.*ro<oifi'-inmio(o ioK>^ma)00'-'*00'>-»-

•r-^i-i-T-T-cNCMCMrorr)

o

OrOtOiOOOOOCM'-CMTl-T)- V'-O5rfr*®r^©oo<ocsi^

_joooio®'r-cMro^-(Dr^rs- o^in<oo)ro®t^.r-mtoi^oo

>-jooi-fomr^oo®fsimoocN O iC-ocM^cocostomr^cncM

<or^i^r^r^ooooooeoooo>

>

co

'-r^ooo>r^csi®'-o)CNi^ csi^tooom®<N^-'-'*r^

CMf>00r--l*)CM«-O><OCOr>

oom'-o>coi^tO'>i-K)'-oo

cooocor>r>.r^i^t^r>r^co

o>oo^®o>f^'— oor~»»-® inmoiooinooooNNoiN IDinOKDNNNMniOK)

I -I

oocN-^-ior^oocoi^to^-ir«-roo><ocoiococor^oo©

oomo)®oo<omcMoo<N<D

oomoo'*-cM®oo«5'*'-r>- cM®o)ooiooooooooor,s»

csiooi^"*if)'r-oor^'^'-®

••*C0'-rO'<*U">ir)<O00®rO '-<o*-'*<or>r>iO'*cMco CMCMK}K>r<->rOK>rOf>?OCM

^csio>ir>'-ir)®®i^ioo) tO'~u*><Dr^®®tocor~®

mo)K)(Dcsi®ooiDTi-ror)

<

CMm

cc

00

h- CM © Z © • LU O) ® 2 ®

(O O • CM 0-

II O

—l'*CMr')in»OCM<0OOCMCM'<l-

—I (DO)«-nif)NO)®MiOiO > <otor^r>r»r>r>oooocooo ©©©©©©©©©©©

<mmr~'-'*o)in'^r^<Ncr)

o. o_^<D OHO'-ifiNONtOintf) T-cMto-^^^minminiT)

o

^mooh-t^®oor^roco'*-'- —i\®o>moo®*-ioco®oom >-3rO«O©K)f>®ror>rOO>C0

o ^w© •^rO-*if)I^COO)'^CM-«l- oocooooooooooooocncAO)

>

CO

CO

X -Ji

o

"*®<-NnMO(DO)MD

<DTt-incoo)®cO'-r^i^o)

r»>T-'^cM"*o>K)co®a)to

oocoooi^i^r^h-r^r^cDco

to^toscMcoincomi^O) m^O)'-rO®^-«-<DIOCM t^O)(DPv«-®rO'-CM(Dm

•r-ir>00©CMrorOrOCM©r^

CMCMCMK)iOr>rororOK)CM

(£)®®Tl-(0^®l/)»-0)0>

(OThT-rOtOS^CM^OOto ^r^a>®mt>^io®T-®

•*l*}CM»-'-®®®CT>O>C0

CMCMCMCMCMCMCMCM*-*-'-

>

>

.mooNooiooiMMntos

rO"-i-~ir>a>®®iooO'*t-^i- co^s^r^cj)®®©)^^"

m<or-*r^r-.i^oooor>r»-i---

C7)CM®C7)®®OOlDCOCOlO

int^sroi^ifimcocMooco

<

CMm

m

I- CM ® 2 ® • LU O) ® 2 ®

® O • ro 0.

II

_i'*ifi'-oo®r-to-*ior^'* —'oo m * k) + io r- id io en m

> t^r^i^cooocoa>cj><7)0505 ®®®®®®®®®®®

<miO'-a)cocO'<iJ-'-<D(D<D

a. o_^-^rooor^'-'-®coio®mr^

CD00®CMtO'*^-tni£>U3CO

o

^r^r^*mo)CM'-oo®'-r^

—l\'J-lf)'-'-mt04DCMO)CMCO >-j®N*ifli/)i/)inioooin(0

o ^

^^ ®®®®®®®®®«— CM 0)OiO)O>O)CX>O)O}O)CX>O>

(D®(Dr-CM®IOO)COOOCM

K1IO* IDNOIIS® M/l'- mm^-cM'-cDooin®'*-®

»-»-o>r^mcM®cocoto®

oooor^i^r^r^r^tDtocotf)

ror^incos'j-'-t^cocMa)

T-rof^ior-»^-«o^-cMmcM r»'-T-mtor^inoo*®CT>

m®tomr^ooo)c»o)oo'*

a>»-o>cO'<-®if)-*'-oir-

>^-voo<£>®'«i-cM'<tr)'*m<ocM X-JOOCMCMOJtOI^S'J-'-m'-

oo>«3<or^®nr^®T»-cO'<i- lO'i-IOCMCMi-SOOCOCO

CMCMCMCMCMCMCMCM'.-*-'- .cO'-«DCM®o)tooo'*toin

00O)lT)'<l-O>®-<f'«l-CM00'* ® to r- cm ao r« lOifi>-O10O

*- t- t^ IO t^ y- * <o r» in cm © *-T-CMCMfOrOf)rOr«->ro

<

CMm

»- CM ® 2 ®

LU CT> ® 2 ®

* O •

II O II h- Lu

•^-CM'-^-inrOrO'OrOino) r^®i)-romo)'J-®i^m^- > •-0)MDifi"*"*'<l, niOIO >

—l(0 + lO(OMn'-IO(D(OOl

»

<o <7>tv-'-'-r^r^m'*®cM

—I iDCD^StDCMOlCDCMtDCO > 00CO<T)®®'-'-CMrOrOtO ©o®«-*-»-«-i-»-«-»-

^n o>'-'-'-®cn'i-cMT-u3 <®SO)ifiC\|^(O00'-'>i-if)

CL 0.i^m®mroo)'i-»oi^'-'^-i^ w® cMinro(oi^m®fO'-^-

COCMlf)00®CM'*CDr^COCO

I 00

Vi9

«

< »- K- C0 2

co r> lj UJ 1- M

2 CC O

< Lu LU LU -IK O

CD 2 O

CQ LU 3 2 O

CD O 2

Q O 2

^® mr^o>moo®if>coo>® ><®'*®r^"-to»-ootO'-© X ®(O®<X>--TfC0'-<£>CM® O®'t(D«)NSM0000)® <

CC T-

—12®®®®®®®®®®® Xv-^OG>®oo®®®®®® ©-cMfO^-miDr^ooo)®

< — t- C0 2

CO n LU

LU .- I-H 2 CC O

< ULU LU

-I I- O

CD 2

CD LU

Z> Z O

CD O 2

\Q_ —

* 2 X

Lu o •— Q O 2

^-v® loooomcMioiooi'-®

><©(OCMCMI"^rOlO"-CM"<»-® X S^COlD'-t^CMCO-*'-® o®>o^ir)(o<or^r^coc7)© <O®®®®®®®©©©® 12©©©©©©©©©©© X ^-^© o®®®®®®®®® ®'-cMtO'*if)U5r^ooo>®

< «- I- CD 2

CO rO LU

LU »- «-i 2 or o

< Li- LU LU

-i I- o

00 2 O

CD LU

Z> 2 O

CD O 2

\0- — gcB^

^s >-®(MinoO(Nonoif)® ><©®cMcor^o)^m<D^-®

X ®^-'*tO»-00<DfO'-®®

o®cMtO'*mif)(or-«cocr>® <

o: — o©©©©®®®©©©© 12©©©®®®®®®®® X^-'®®®®®®®®®®® ©^cMro'j-mKDr^coCTi®

o o

:

141
NBS-114A (rev. 2-8C) U.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO. KBS/TN--1226

2. Performing Organ. Report No. 3. Publication Date August 1986

4. TITLE AND SUBTITLE

Application of a Hard Sphere Equation of State to Refrigerants and Refrigerant

Mixtures

5. AUTHOR(S)

Graham Morrison and Mark McLinden

6. PERFORMING ORGANIZATION (If joint or other than NBS, see instructions)

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

Gcdther-bura, n:> ~ n "?~

7. Contract/Grant No.

8. Type of Report & Period Covered

Final

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

Electric Power Research Institute

P.O. Box 10412

3412 Hillview Avenue

Palo Alto, CA 94303

10. SUPPLEMENTARY NOTES

"2 Document describes a computer program; SF-185, FlPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bibliography or literature survey, mention it here) This note describes the application of the Carnahan-Starl ing-DeSantis equation

of state to halogenated hydrocarbon refrigerants and their mixtures. A complete and consistent set of thermodynamic functions is derived from the p-V-T

equation of state and the perfect (ideal) gas heat capacities. A thorough

discussion of reference states is included for both pure materials and their

mixtures. Although this model exhibits a critical point, it does not quantitatively represent properties in the critical region. Despite this limitation, this model can represent both liquid and gaseous mixtures away from

their own critical points, even at conditions near to and above the critical

points of their components.

Algorithms and FORTRAN routines for the use of this model are presented along

with the numerical coefficients for 11 pure refrigerants and 7 mixtures.

Routines for evaluating the coefficients from saturation data are included.

Several examples of the application of this equation of state are presented to

demonstrate its versatility. The average deviation from the tabulated saturation properties of the 11 pure refrigerants is 0.54% for pressures, 0.09% for

liquid volumes, and 0.50% for vapor volumes. It is shown to predict other

properties of pure materials well and to describe the detailed features of _ mixtures, both phase diagrams and thermodynamic properties.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolons)

equation of state; hard-sphere fluid; mixing rules; mixture; phase behavior;

reference states; refrigerant; thermodymamic properties

13. AVAILABILITY

\x~\ Unlimited

| | For Official Distribution. Do Not Release to NTIS

P^l Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402.

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

14. NO. OF

PRINTED PAGES

15. Price

USCOMM-DC 6043-P80

DATE DUE

3-23*9*

Demco. Inc. 38-233
NBSTechnical Publications

Periodical

Journal of Research—The Journal of Research of the National Bureau of Standards reports NBS research and development in those disciplines of the physical and engineering sciences in which the Bureau is active. These include physics, chemistry, engineering, mathematics, and computer sciences. Papers cover a broad

range of subjects, with major emphasis on measurement methodology and the basic technology underlying

standardization. Also included from time to time are survey articles on topics closely related to the Bureau's

technical and scientific programs. Issued six times a year.

Nonperiodicals

Monographs—Major contributions to the technical literature on various subjects related to the Bureau's scientific and technical activities. Handbooks—Recommended codes of engineering and industrial practice (including safety codes) developed in cooperation with interested industries, professional organizations, and regulatory bodies.

Special Publications—Include proceedings of conferences sponsored by NBS, NBS annual reports, and other

special publications appropriate to this grouping such as wall charts, pocket cards, and bibliographies.

Applied Mathematics Series—Mathematical tables, manuals, and studies of special interest to physicists,

engineers, chemists, biologists, mathematicians, computer programmers, and others engaged in scientific and

technical work.

National Standard Reference Data Series—Provides quantitative data on the physical and chemical properties

of materials, compiled from the world's literature and critically evaluated. Developed under a worldwide pro- gram coordinated by NBS under the authority of the National Standard Data Act (Public Law 90-396).

NOTE: The Journal of Physical and Chemical Reference Data (JPCRD) is published quarterly for NBS by

the American Chemical Society (ACS) and the American Institute of Physics (AIP). Subscriptions, reprints,

and supplements are available from ACS, 1155 Sixteenth St., NW, Washington, DC 20056.

Building Science Series—Disseminates technical information developed at the Bureau on building materials,

components, systems, and whole structures. The series presents research results, test methods, and perfor- mance criteria related to the structural and environmental functions and the durability and safety

characteristics of building elements and systems.

Technical Notes—Studies or reports which are complete in themselves but restrictive in their treatment of a

subject. Analogous to monographs but not so comprehensive in scope or definitive in treatment of the subject

area. Often serve as a vehicle for final reports of work performed at NBS under the sponsorship of other government agencies.

Voluntary Product Standards—Developed under procedures published by the Department of Commerce in Part 10, Title 15, of the Code of Federal Regulations. The standards establish nationally recognized re- quirements for products, and provide all concerned interests with a basis for common understanding of the

characteristics of the products. NBS administers this program as a supplement to the activities of the private

sector standardizing organizations.

Consumer Information Series—Practical information, based on NBS research and experience, covering areas of interest to the consumer. Easily understandable language and illustrations provide useful background

knowledge for shopping in today's technological marketplace.

Order the above NBS publications from: Superintendent of Documents, Government Printing Office,

Washington, DC 20402. Order the following NBS publications—FIPS and NBSIR s—from the National Technical Information Service, Springfield, VA 22161.

Federal Information Processing Standards Publications (FIPS PUB)—Publications in this series collectively constitute the Federal Information Processing Standards Register. The Register serves as the official source of

information in the Federal Government regarding standards issued by NBS pursuant to the Federal Property

and Administrative Services Act of 1949 as amended, Public Law 89-306 (79 Stat. 1127), and as implemented

by Executive Order 1 1717 (38 FR 12315, dated May 11, 1973) and Part 6 of Title 15 CFR (Code of Federal

Regulations).

NBS Interagency Reports (NBSIR)—A special series of interim or final reports on work performed by NBS

for outside sponsors (both government and non-government). In general, initial distribution is handled by the sponsor; public distribution is by the National Technical Information Service, Springfield, VA 22161, in paper copy or microfiche form.
U.S. Department of Commerce

National Bureau of Standards

Gaithersburg, MD 20899

Official Business

Penalty for Private Use $300

!' H

lFr«K

Mi

,>>-*v

>V<\S