A111D2 7flTDSfi

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

A1 11 02789058

Eaton, B. E/On the calculation of critic QC100 .U5753 N0.1313 1988 V198 C.1 NBS-P

Siimulaiing America's Progress 1913-1988

NBS

.ICATIONS

NBS TECHNICAL NOTE 1313

U.S. DEPARTMENT OF COMMERCE / National Bureau of Standards

OC

100

U5753

//1313

1988

C.2
TM he Nauonal Bureau of Standards' was established by an act of Congress on March 3, 1901. The Bureau's overall ^t 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 to assure international competitiveness and leadership of U.S.

industry, science arid technology. NBS work involves development and transfer of measurements, standards and related science and technology, in support of continually improving U.S. productivity, product quality and reUability, innovation and underlying science and engineering. The Bureau's teclmical 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,

industry, and commerce; provides advisory and research services to other Government agencies; conducts physical and chemical research; develops,

produces, and distributes Standard Reference Materials; provides

calibration services; and manages the National Standard Reference Data

System. The Laboratory consists of the following centers:

• Basic Standards^

• 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. Tlie Laboratory consists of the following centers:

• Applied Mathematics

• Electronics and Electrical Engineering^

• Manufacturing Engineering

• Building Technology

• Fire Research

• Chemical Engineering^

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 divisions:

Information Systems Engineering

Systems and Software

Technology

Computer Security

Systems and Network

Architecture Advanced Computer Systems

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-cutting 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'

• Polymers

• Metallurgy

Reactor Radiation

'Headquarters and Laboraiories at Gaithersburg, MD, unless otherwise noted; mailing address

Gailhcrsburg, MD 20899. 'Some divisions within the center arc located at Boulder, CO 80303.

'Located u Boulder, CO, with some elements at Gaithersburg, MD
Research Information Center

National Bureau of Standards

Gaithersburg, Maryland 20899

'.^ ~. ^ '-^ On the Calculation of Critical Points ^^^^ by the Method of Heidemann and Khalil

Brian E. Eaton

Thermophysics Division

Center for Chemical Engineering

National Engineering Laboratory

National Bureau of Standards

Boulder, Colorado 80303-3328

^s.*'-'

^ ^-i^

•"•tAU O* *

U.S. DEPARTMENT OF COMMERCE, C. William Verity, Secretary

NATIONAL BUREAU OF STANDARDS, Ernest Ambler. Director

Issued March 1988
National Bureau of Standards Technical Note 1313

Natl. Bur. Stand. (U.S.), Tech Note 1313, 52 pages (Mar. 1988)

CODEN:NBTNAE

U.S. GOVERNMENT PRINTING OFFICE

WASHINGTON: 1988

For sale by the Superintendent of Docunnents, U.S. Government Printing Office, Washington, D.C. 20402
CONTENTS

Page

1 . Introduction 1

2. Conditions for the Local Stability of a Homogeneous Phase 1

3. Alternative Formulations of the Conditions for a Critical Point 5 3.1 The Quadratic Form 5 3.2 The Cubic Form 12 3.3 Computational Considerations 1

i)

4. The Method of Heidemann and Khalil 15

4,1 Analytical Expression for the Volume Derivative of the Cubic Form ... 19

5. Implementation for a Two-Constant Cubic Equation of State 20

5.1 Functions of the Dimensionless Volume K and Their First-Order

Derivatives 26

6. How to Use the Critical Point Software 27

6.1 Verification of the Software 28

7. Summary of FORTRAN Subprograms 30

8. References 44

ill

ON THE CALCULATION OF CRITICAL POINTS BY THE METHOD OF HEIDEMANN AND KHALIL Brian E. Eaton

The formulation of critical point criteria by Heidemann and Khalil is analyzed

and contrasted with the original formulation of Gibbs. An extension to the solution

technique originally used by Heidemann and Khalil, and later improved by Michelsen

and Heidemann, is presented along with its detailed implementation for a general

two-constant cubic equation of state. Finally, FORTRAN software developed for these

computations is carefully documented.

Key words: algorithm; critical locus; Helmholtz free energy; mixtures; Peng-Robinson

equation of state; thermodynamics

1 . Introduction

The topic of this technical note is the computation of critical points in multicomponent

mixtures using equations of state. The thermodynamic conditions for a critical point were orglnally formulated by Gibbs \3}i but because of their complexity, they were not solved

using equations of state for multicomponent mixtures for another 100 years. Peng and

Robinson \10} computed critical points in multicomponent mixtures using their cubic equation

of state and the conditions for criticality exactly as given by Gibbs \3}, i.e., Gibbs'

equations (206-208). Shortly thereafter, Heidemann and Khalil \4} proposed an alternative

computational procedure based on a mathematical reformulation of the original conditions of

Gibbs. The computational procedure of Heidemann and Khalil is substantially more efficient

than that employed by Peng and Robinson.

The purpose of this note is twofold. First, the criteria for a critical point will be

examined in some detail. Emphasis is placed on elucidating the similarities and differences

between the two mathematical formulations mentioned above. Second, a description of Heidemann

and Khalil' s solution technique using a general two-constant cubic equation of state is presented along with an extension which may lead to a more efficient iteration on the volume.

In the final section, the FORTRAN software developed for the computations is carefully

documented.

2. Conditions for the Local Stability of a Homogeneous Phase

To formulate the mathematical conditions for a critical point, we begin by considering

the stability of a homogeneous phase with respect to continuous change. Implicit in this

statement is the differentiation between an arbitrarily small (continuous) change in a phase

and a change which involves the appearance of a new phase (discontinuous). A phase may be

stable with regard to the former type of change and unstable with regard to the latter (i.e.,

a metastable state). Stability with respect to continuous changes is referred to as local

stability and that with respect to discontinuous changes is termed global stability.

Thermodynamic stability theory rests upon a hypothesis which, in its original form as stated by Gibbs, does not lead unambiguously to intuitive concepts and mathematical formalism;

an excerpt from Tisza \1iJ} (p. 41) clarifies this point:

Gibbs stated the extremum principle in two versions: in an isolated system the

entropy tends to a maximum at constant energy, or alternatively, the energy tends to

a minimum at constant entropy. Although these statements undoubtedly express

important truths, they lack precision to the point of being paradoxical. If an isolated system is not in equilibrium, we can associate no entropy with it, and if
it is in equilibrium, its entropy can no longer increase. Many authors have

grappled with this dilemma until a satisfactory solution was found in terms of the

composite system. Consider a system of two or more spatially disjoint parts

separated by adiabatic partitions. After reaching equilibrium each part has a definite entropy, the sum of which is associated with the composite system. The

relaxation of an internal constraint will, in general, trigger a process, namely the

redistribution, say, of energy at constant volume,

(DU-)^, + (DU")^„ =

leading to a new equilibrium with higher entropy. Here D denotes a finite change

that actually takes place in the system and for which all the conservation laws and

boundary conditions are satisfied. Thus the increase of entropy is perfectly well

defined since it is associated with a transition from a more restricted to a less

restricted equilibrium.

Development of the stability conditions from the extremum principles and the composite

system concept summarized below can also be found in Tisza [14] and in Modell and Reid [8].

It is convenient to consider the results in terms of the energy rather than the entropy since

the Legendre transforms of the energy result in the free energy expressions which are readily

calculated using equations of state.

The extensive state of a system can be represented in terms of c + 2 variables, where c is the number of independent components,

U(S, V, n^ n^_^ , N); (2.1)

the internal energy (U) is a function of the entropy (S), volume (V), number of moles of

component i (n^) and the total number of moles (N). It is generally assumed that the

stability of a system is an intrinsic property which does not depend on its size. This

assumption is not universally valid: for example, it is not valid for small droplets, where

surface tension effects are important, nor is it valid for large stars, where gravity must be

taken into account. However, size dependent effects will not be considered here; hence, the

size of the system will be fixed by holding, say, N constant. Variations in the energy are then expressed by

1=1 1

where the variables held constant in the partial differentiations follow from the functional

dependence expressed in (2.1). Introducing some additional nomenclature, we rewrite (2.2) as follows,

dU = I f-dX. = I U. dX. , (2.3)

1=1 1 1=1

where X^ is the i^*^ component of the vector X = [S, V, nj^, ..., ng--]]'^,

U. . — , and r . c . 1 .

The Uj_ represent the intensities of the system, that Is U^ is the i'''^ component of the vector

P = [T, -p, PI, ..., Uq-i ]T where T is temperature, p is pressure, and \ii is the chemical

potential of mixture component i. The superscript T represents the transpose operation, i.e.,

X and P are column vectors.
The extremum principle is expressed as follows. Consider a composite system, the

subsystems of which are identified by the superscripts a and 6. The initial state, in which

subsystems a and B are assumed to be in a state of unconstrained equilibrium (all the

intensities in a equal the intensities of g) , is characterized by

^jsystem

^

^^a ^ ^B _ (2_^)

The constraint that the composite system be isolated takes the form

x" + X^ = const., i = 1, .... r (2.5)

C( S and the relative sizes of a and 6 are specified by fixing the valves M and N . The allowable virtual processes used to test the stability of the system can be visualized as

an exchange of the quantities Xj^ between the subsystems a and B subject to the constraints of

(2.5). Thus, a virtual process, or a displacement from the equilibrium state, has the

property that

6X°' = - 6X^ , i = 1 r (2.6)

where the symbol 5 denotes an arbitrarily small displacement. According to the energy

extremum principle (as stated by Tisza), any virtual process which takes the composite system

from an initially unconstrained equilibrium to a constrained equilibrium must result in a higher energy for the composite. The response of the system's energy to such virtual

processes is expressed by a Taylor series expansion about the initial equilibrium state, i.e.,

= 5U" . i 6^ U" ^ i 63 U" + . . . (2.7)

2 D

. 5U^ . i 6^ U^ . 1 6^ U^ . . . . 2

where

6U =- I f- 6X. -. I U. Z. , (2.8)

1=1 1 1=1

r 2 r 6^U = I ^^ 6X 6X = I U Z Z . (2.9)

i,j = l ^^^^j ' ^ i,j = l 'J

' ^

^^" ^

. \, axixV ^^i^^j^^k ^ . \^ ^jk^i^j^k • (2.10)

i,J,k = l 1 J k

-^ i,J,k = l

"i = "3X7 ' '^ij = 3X71X7 • "ijk = 3X.9X.9X, ' ^""^ ^ = ^^i " 1 1 J 1 J k

Substituting (2.6, 2.8-10) into (2.7) gives

I (u" - uh z'^ . i I (U". ^ uf.) Z° Z"

1=1 1112 ._j^^ ij ij i J (2.11)
a S Since a and B were initially in equilibrium, U. = U.; thus the first order form is identically

zero. The quadratic and cubic forms are simplified by the following consideration. Since the

energy is represented by a first order homogeneous function,

U" [X^, X^...., X^] = k U^ [X^^ X^ X^^] (2.12)

i^ ft where k is the ratio of the sizes of two systems which can be denoted by N /N . Also,

X^ = k X^ . (2.13)

The following relationships between the partial derivatives of U and U follow from

differentiating (2.12) and making use of (2.13):

U" = U^ , (2.1H)

1 1

kU°'. = U^. , (2.15)

iJ . iJ

k^ U";*., = U^., . (2.16)

ijk ijk

Substituting into (2.11) gives

AU , = H^ I U". Z^Z'^ + 1^ I U", Z"Z"Z" + . . . (2.17) system 2 ";_. 13 ^ J 6 /- ijk 1 j k

The extremum principle requires that for an initially stable system, all virtual processes

must result in the change

AU ^ > . (2.18) system

ot Hence, the quadratic form, which dominates the expansion as the virtual displacements Z. are

made arbitrarily small, must be positive, i.e.,

6^U > . (2.19)

In (2.19) we divided by the positive constant (1 + k)/2, and dropped the superscript a since

it is irrelevant which of the subsystems is considered.

The limit of stability is determined by the condition

6^U = . (2.20)
This condition is of particular interest in this discussion because a critical point is a stable point on the stability limit. The stability of a state for which 5^11 = is determined

by the higher order forms. Consider the cubic form. If a virtual process exists such that

qSu > 0, then because the cubic form is an odd function of Z^, there also exists a process

(-Z^) for which 63u < 0. Therefore, the stability of a critical point requires that

o^U = and 6^U > . (2.21)

In a pathological case where 6^U = 0, a similar argument to that used for the cubic form would

require that a^U = and 6^U > 0.

In fact, it is normal procedure when computing critical points to consider only the

conditions &^{i = and &^U = 0, and to rely on experimental evidence to infer that a critical

point actually exists (i.e., is stable).

Before moving on to alternative f orm.ulations of the stability conditioi'is, there is one more point to be made concerning our composite system model. In the preceding derivation we found the relative sizes of the subsystems a and S, expressed by the variable k, to be

irrelevant to the stability conditions. There is a special value of k, however, which adds

valuable intuition about the kinds of virtual processes which may be used to determine

stability.

Consider the case where one of the subsystems, B, say, becomes very large. In the limit

as N^ ^ ", k ^ 0, and by eqs (2.15-16) the second and higher order derivatives of U*^ go to zero. Recognizing that the second order derivatives of U are the same as the first order

g derivatives of the intensities, U.. = implies that the intensities of 6 remain constant

during any virtual process. This corresponds to the concept of ideal reservoirs. For example, the temperature of a thermal reservoir remains constant when heat (entropy) is added

to or removed from it. Thus, a composite system consisting of a subsystem under study and an ideal reservoir leads naturally to the consideration of virtual processes in which the

intensities are held constant.

3. Alternative Formulations of the Conditions for a Critical Point

We will consider a critical point to be determined by the conditions

6^U = 6^U = . (3.1)

Reformulations will be considered first for the quadratic form and subsequently for the cubic

form.

3.1 The Quadratic Form

A stable system is given by the condition 6^11 > 0; in other words the quadratic form is

positive definite. We shall eventually find it useful to express the results concerning the

positive definiteness of the quadratic form in terms of matrix notations. Hence, the notation

is introduced here,

6^U = z"^ U Z , (3.2)

where Z_ is the column vector (5S, 6V, 6ni , ..., 6nj,_i)^,
3^U

U is a symmetric matrix of elements U. . = , = 1

J

d A . A . 1 J

and T indicates the transpose operation.

A common technique used to derive conditions which determine the definiteness of a quadratic

form is to diagonal ize U using orthogonal transformations. However, orthogonal

transformations do not turn out to be very useful in thermodynamics because the variables do not readily admit the definition of a metric (see, for example, V/einhold [15]). The most

useful procedure for diagonalizing U, from the viewpoint of thermodynamics, makes use of the

technique completing-the-square. The diagonalization of U and subsequent reformulations of

the stability conditions will now be considered in some detail for a binary system.

The first step in writing the quadratic form for a binary system as a sum of squares is

2 ^

-2,. „ ,- 1 f„ „ „ ^ M 1

where

Continuing to complete the square,

,(1) -.2 r .,,(1) ^

6^U -

or

where

(3.3)

"11 ^f

^ "22^ 4 ^ 2 41^ Z2Z3 . U^) Z^ , (3.4)

^ ^ ^1 ^u;; ^"12^2 ^"13^3^ • ^^-^^

"22^ =- "22 -U^' (3.6)

„(!) = ,, _

'"12 "l3 ,. 7^ "23 - "23 u^^ • (3-^^

"33' -^ "33 -uj;- ''-''

.. 2 .,(1) r_ ^"23 „1 , r,(i) ("23 ^ 1 .2

22 22

6^U = U^^ n^ * U^2^

^l

^ U^^^ Z^ , (3.10)
22

,(2) ^ „(1), ^"23

( ,2) ^33 33 ,,(1) ^^ * 22

Equation (3.10) implies that for 6^U to be positive definite, the conditions

ui°^ = "ii > . (3.13)

U^2^ > , (3.n)

U^^^ > , (3.15)

must be satisfied.

The first step toward reformulating the stability conditions is to recognize that in a stable system moving toward an unstable state, the first quantity which will go to zero, and

(2) thus indicate the limit of stability, is U__ . It is common to make use of this single

condition when searching for the stability limit of a system of interest and to use intuition

gained from experimental evidence to guarantee that the search is started from a stable state. This is a practical approach. The purpose of indicating the implicit assumptions normally

made is that when working with equations of state and Newton-Raphson search techniques on the

computer, we can unwittingly project the solution into a region where the lower order

stability conditions are not met and hence produce incorrect results.

(2) To prove that the condition U is the first to go to zero, consider (3.6) and (3-12)

which indicate a general recursion formula relating successive stability conditions and can be

written as:

r (k-2) ^^ ,,(k-l) „(k-2) ^ k(k-1)^ ,, .,. \k = "kk "7i^^2) • (3.16)

(k-l){k-1)

(k-2) (k-1 ) Under the assumption that the system is initially stable, U, >,,_,, U , and, hence,

U {k-2) (k-2) (k-2)

.. are all positive. If we let 'Jfi.-i wu-i ) 2° to zero, then assuming that U remains

,(k-2) ,2 ^ ^ ^ ,,('<"'')

J , ,) does not go to zero, U

K ^ K I / KK

(k-2) (k-1

have become zero before U, ,, , reaches zero. Therefore, U will become zero
(k-2) before U, ,, , does as the stability limit is approached.

\K i J \K \ J

If we assume that U > 0, and use the recursion formula (3-16), then by induction the

condition which will go to zero first, and thus indicate the limit of stability, will

(r-1 ) be U , where r equals the number of rows and columns in U. rr

(r-1 ) The fact that the stability limit is determined by the condition U =0 allows us to

rr deduce something about the types of virtual processes which may be used to test stability. In our model of the composite system we have imagined virtual processes (Z) as the exchange of

extensive quantities X^ between the subsystems a and B. The diagonalization of the quadratic

form for the binary system resulted in a transformation from the virtual processes (Z-], Z2,

Z3) to the processes ( m , ti2> Z3). In general, diagonalization will transform the processes

(Zi, Z2, ..., Zp) to the processes (ni, 112. •••> '^r-l > Zj,), Since the limit of stability is determined by the coefficient of Zp, the virtual processes used to test for stability need

only include the subset given by (ni = (i = 1, ..., r-1), Z^) . That is, we may arbitrarily

set the variations represented by the m to zero. Let us clarify exactly what types of

processes the m correspond to. Consider the connection between variations in the extensive and intensive variables for a binary system, i.e.,

"ss "sv

"vs "vv

Sn,

Vn,

n S n V n n

6S

6V

6n,

-P,

'IS 'IV In.

6S 6T

6V = -6P

«"l 6m^

(3.17)

A process at constant temperature implies that

"SS

^S ^ "SV '' ' hn^ ^"1 5T = , (3.18)

or, converting notations

"11 ^1 ' ^2 h ' ^3 h -

° (3.19)

and, since U-| 1 > for a stable system, we divide by it to obtain

1 ^1 ^117; ^"12 ^2^ "13 S^ . (3.20;

Comparing (3.20) and (3-5) leads us to the conclusion that m = is equivalent to a virtual

process at constant temperature.

Since the stability limit may be determined by using constant temperature processes, it

is natural to reformulate the quadratic form in terms of the Helmholtz free energy, i.e..

6h I A. , Z.Z. (3.21)
where

Diagonalizing leads to

(6V, dn^)"^ . (3.22)

2 2 (1)2 6^A = A^^ n^ + A^2 ^2 . (3.23)

And the limit of stability is determined by

(1 ) 12 4 = ^22 -17;- = °

• (3.2^)

The stability limit takes a simpler form in terms of A because the constant temperature

constraint is implicit in A which is the first Legendre transform of U with respect to S. Let us carry the present line of thought one step further. Consider virtual processes

which occur at constant temperature and pressure. The constraints placed on the variations of

the extensive variables are

^1 h ' "12^2^ U^3Z3 = 6T = ,

^21 h ^ "22 ^2 ^ "23 ^3 = -«P = •

Eliminating Z-| from these equations gives

Z + yZ = . (3.25)

"12

U - —-=- 22 U^^

Comparing (3.25) and (3.11) shows that n2 = is equivalent to a process at constant

temperature and pressure. Hence, the stability conditions can be expressed in the Gibbs free

energy formulation, i.e.,

6^0 = G^^ Z^ , (3.26)

where

^1 = ^"1 •

The limit of stability is given by the condition

G^^ = . (3.27)

For multicomponent mixtures we can generalize the conditions for the limit of stability

as follows. In the energy representation,

U^^'^^ = . (3.28)

rr
Considering virtual processes at constant T leads to

(r-l)(r-1) "

• (3.29)

and considering processes at constant T and p leads to

^(r-2)(r-2) " ° (3.30)

Beyond this point we run out of commonly used symbols to represent the Legendre transforms of

U, hence, the notation [U^] is introduced to represent the i^^ Legendre transform of U. Considering processes at constant T, p, p-| leads to

. (r-4)

u- ; (3.31)

(r-3)(r-3)

and considering processes where r-1 intensities are held constant, i.e., constant T, p, \i-\ ,

. . . , Vq-2 I leads to

(0) r-1 [U^"'] (3.32)

11

There are r equivalent formulations of the stability conditions for the particular

ordering of independent variables we have chosen. Many more formulations are possible by

taking different permutations of the order S, V, n-] , ..., nQ_i . Of all possible formulations,

there is really only one which interests us here; that is the formulation in terms of the

Helmholtz free energy because most equations of state are written with T, V, n-| , ..., n^ as the independent variables.

In addition to reformulating the stability conditions in terms of the Legendre transforms

of the energy, there is another important reformulation which we consider now. As before, the

motivation will come from working with the specific example of the binary system.

Consider the stability conditions given in eqs (3.13-15). The first condition can be

trivially rewritten

U

(0)

11 u11 M

• :3.33)

where A-| is the first discriminant of the matrix U. The m^'^ discriminant of a matrix

is the determinant of the submatrix obtained by deleting all elements which do not simultaneously lie in the first m rows and columns (see, for example, Hildebrand, p. 52).

The second condition can be rewritten.

U^^ ^ = —^ (U U - U^ ) 22 U ^11 22 12-^ (3.34)

or, due to the symmetry of U,

/I)

22

The third condition can be rewritten.

11

"l1 ^12

"21 "22

3.35)

10
if ' ("33

"l.

"23 U,

"22 U,,

^"ll"22 - "l2^^"ll"33 " "13^ - ^"ll"23 - ^12^13^'

^"11^22 -^12) "^^

The quadratic form for a binary system can then be written,

and the condition for the stability limit translates into

For the binary system, U is a 3 x 3 matrix, hence (3.40) becomes

In general.

(3.36)

(3.37)

(2) ^3 "33 = A^- (3.38)

2 2 "^2 2 ^3 2

5 U = A^ n^ .- n^ -^ Z3 , (3.39)

A^ = . (3.40)

A = det j U

I = . (3.41)

1 2 r-2 r-1

and the stability limit can be represented by

U^^'^^ = A^ = det I y I = . (3.43)

The condition for the stability limit expressed as Ap = is in the final form as given

by Gibbs and is in the form used for computation by Heidemann and Khalil. Equivalent forms

are obtained when working with a Legendre transformation of the energy representation. In general, the stability limit can be expressed

. (r-i-1)

[U^] = A^ . = . (3.44)

(r-i)(r-i)

Before moving on to a discussion of the cubic form, there is one final point to be made

regarding the quadratic form which is suggested by the condition for the limit of stability

written as a determinant. Consider the energy formulation. At the stability limit,

eq. (3.43) implies that U is a singular matrix. We have mentioned previously that U represents a transformation between the extensive and intensive system variables. When a sys11
tern is in a regular state (i.e., U is not singular) then a set of r extensive variables

uniquely specifies a set of r intensive variables (recall that x^+i has been fixed to specify

the size of the system). This is in accord with Gibbs' phase rule, i.e., the number of

independent intensive variables in a single phase system is e+1 = r. At the limit of

stability, the rank of U is reduced from r to r-1 ; hence, there are only r-l independent intensive variables. Since U is singular, there exists a non-trivial solution to the homogeneous equation

y

Z = . (S.'JS)

Thus, the virtual process which causes the quadratic form to be zero at the stability limit

defines a unique direction in Gibbs space within an arbitrary sign; this process is referred

to as a critical displacement. Furthermore, the critical displacement expressed by eq (3.^5)

represents a virtual process in which the r-1 independent intensities are held constant

(recall eq (3.17)).

3.2 The Cubic Form

The condition 6->U = is necessary, but not sufficient, to prove the stability of a

system on the stability limit. In this section two formulations of the condition 6^U = will

be presented. Analogous formulations in terms of the Legendre transforms simply involve fewer

independent variables.

The previous section established that, for a system on the stability limit, there exists

a virtual process, unique within an arbitrary multiplicative constant, called the critical

displacement for which the quadratic form becomes identically zero. Therefore, the critical

displacement is the only process which needs to be considered when evaluating the cubic form.

The cubic form is a special case of a general tri-linear form. Its linearity allows it

(2.10) to be written as

^'" - X h-i- t, ^ WIT ^i^J^^3.46)

k=l k i,j=1 1 J

or

&h = I \^ CS^U] . (3.^*7) k=1 k

or

6^U = Z-V (6^U) , (3.48)

12
where V is a gradient operator with components -^r— . The operator Z»V produces a directional

derivative, that is, the derivative in the direction of the vector 1_ in Gibbs space. (Gibbs

space is taken here to be spanned by the variables S, V, n-| , ..., n^--^.) Hence, eq (S.'iS)

indicates that the cubic form necessary to show the stability of a system on the stability

limit is the derivative of the quadratic form in the direction of the critical displacement.

To evaluate the cubic form, the direction of the critical displacement vector must be known;

this vector is determined for systems on the stability limit by solving eq (3-^5).

Gibbs formulated the second condition for a critical point in a way which avoided having

to explicitly determine the critical displacement vector. His formulation is more general

than is the specification that S^U = 0, but it contains a redundancy which is extremely costly

from a computational point of view. To formulate Gibbs' second condition we proceed as

follows.

The critical displacement is one in which r-1 of the r intensities are held constant. To

consider virtual processes of this type, choose any r-1 of the r equations represented by

eq {3-^5) to act as constraints on the r independent extensive variations 1^, i = 1, ..., r.

For definiteness choose the first r-1 equations. These equations impose the constraint that

3U the r-1 intensities -i-rr- , i oX .1

1 r-1, are constant in any virtual process. The only

independent variation is 6Xp which may be chosen to be arbitrarily small as required to

determine local stability. For the final equation consider an analog to the cubic form; that

is, rather than take a directional derivative of the quadratic form, differentiate the

discriminant A^ which is an equivalent indicator of the stability limit. Just as Z'V(62u)

must be zero for a critical point, so must Z_'V(Ap). Then our final set of equations is

U

11 U12

U

1r

"r-,,, V,,2

L

3A

r 3X.

dA

r 3X,

U , r-1 ,r

3A

r 3X

6X

6X

r-1

6X

. (3.iJ9)

A nontrivial solution to eq (3.^9) requires the matrix to have a zero determinant, which is

Gibbs' second condition for a critical point.

13
3.3 Computational Considerations

Let us first consider the computations required by Gibbs' formulation for a critical

point. The conditions are that Ap = (eq (S.'^S)) and that A^ = 0, where A^ is the

determinant of the matrix in eq (3- '^9). Both conditions involve r^^ order determinants,

however, the second condition contains derivatives of A^ within the determinant A^. An exact formula for the derivative of a determinant is derived by expressing the determinant as

a multiple matrix U, then (Courant and John [1])

A = det I U I = U, X U^ X . . . X U -U (3.50)

r 1=1 -1-2 -r-1 -r

where x is a vector cross product and • is a vector dot product. Differentiating (3.50) gives

3A^ 3U^ ^ = 33^ X U^ X ... X U^_^.U^

1 1

9U

. U^ X 3^ X ... X U^_^.U^ (3.51)

1

3U

. U^ XU^X ... XU^_^.^ ; 1

hence, the derivative of an r^^ order determinant equals the sum of r determinants, in each of

which one column is replaced by its derivative with respect to the appropriate independent

variable. Hence, to evaluate Gibbs' conditions for a critical point of a single state

requires the evaluation of (r^ + 2) r^^ order determinants. The most efficient way to compute

a determinant is to factor the matrix and then compute the product of the pivots. Matrix

factorization requires approximately r3/3 operations; therefore, the computational effort

required to evaluate Gibbs' conditions is proportional to r^.

The. efficiency of the method proposed by Heidemann and Khalil for computing critical

points comes from their realization that if eq (3.^5) is solved for the critical displacement

vector Z, then the cubic form may be evaluated directly. Direct evaluation of the cubic form

requires r3 operations. It will be shown later that for some equations of state, such as the

general two-constant cubic equation to be considered, the cubic form can be analytically

summed in such a way as to reduce the computational effort required for its evaluation to being proportional to r^. When the cubic form may be simplified in that way, the major

portion of the critical point calculation is required to factor the matrix U; thus the overall

effort remains proportional to r3. In addition to the difference in computational effort required by the two formulations

presented above, they require distinct iteration schemes which will be discussed here. The

conditions of criticality will be expressed in terms of the Helmholtz free energy, and the set

of independent variables is taken as (T, V, n-] n^-i). Since there are two conditions of

criticality, we will specify the composition of the mixture of interest and then solve two equations for the values of T, V at a critical point. Critical pressure is subsequently

obtained from the equation of state.

14
Gibbs' formulation results in the equations

Vl = ° • ^^-1 = ° • ^3.52)

These equations may be solved using a Newton-Raphson iteration scheme. The derivatives of the

determinants required for the iteration are normally computed numerically. Equation (3-52)

may be evaluated at any state point, and the critical point results will be correct as long as

the stability limit is approached from a stable region.

Heidemann and Khalil formulated the critical point as

A^_^ = , 6^A = . (3.53)

The solution procedure must be modified from a straightfoward Newton-Raphson technique because

the cubic form can be evaluated only at a state on the stability limit. Therefore, if the

volume is specified, the temperature on the stability limit may be determined by solving

A)---] = using a one-dimensional Newton-Raphson iteration. Then the critical displacement

vector is determined by solving A Z_ = 0, and the cubic form is evaluated. The search has now

been reduced to a one-dimensional search for the critical point along the limit of stability.

The derivative of the cubic form with respect to volume along the stability limit is determined numerically, and a Newton-Raphson Iteration employed to update the volume. This

procedure of nested one-dimensional searches is continued until a critical point is found. In addition to the efficiency of computing &^k rather than A^_-| , the formulation of Heidemann

and Khalil appears to have a larger radius of convergence. Heidemann and Khalil reported that

a single strategy for the initial guess was sufficient to compute critical points for the

large group of mixtures they considered.

iJ. The [lethod of Heidemann and Khalil

Heidemann and Khalil chose to work with the stability conditions written in terms of the

Helmholtz free energy; this is the logical choice when working with equations of state with T,

V, n-i n^ as the independent variables. Thus, all virtual processes are carried out at constant temperature. Furthermore, Heidemann and Khalil fixed the size of the system by

arbitrarily holding V constant, rather than N as was done in the theory sections of this note; this results in the independent variable set ni , n2, ..., n^, where c is the number of

components in the mixture. Their claim that this choice of variables "has the effect of

producing symmetrical quadratic and cubic forms" is misleading. The symmetry of the quadratic

and cubic forms is predicated on the order of differentiation in the second and third order

derivatives of the free energy being arbitrary. In general, the variables n-j , n2. .... n^ are convenient from the point of view of computations because all the derivatives of the free

energy will have the same form. Furthermore, for the particular equation of state to be

examined in this work, this choice of variables results in significant simplifications in the

analytical expressions of the quadratic and cubic forms.

In terms of the Helmholtz free energy (A), the stability limit is expressed by the

condition (dropping the subscript on Ar>_i )

A s det I Q I = (1.1)

where the components of the c x c matrix Q are

15
1 J

For a specified value of V, eq (4.1) is used to determine the value of T on the stability

limit. Because of the expense of computing -r— analytically, and because the computation of A

o 1 itself is the most time consuming part of the critical point computation, a quasi-Newton

iteration is used to solve [H.]). This technique involves computing -r;r numerically on the

oT first iteration, i.e.,

9A(T ) A(T + 6T) - A(T ) ^ = ° °- ' (1) 3) 3T 6T ^ '^'

where 6T = 10~° Tq, Tq is the initial temperature guess, and A(T) indicates the determinant

evaluated at T (V, n-\ , . , .

, n^ are fixed during the temperature iteration). On subsequent

iterations use

3A(T.) A(T.) - A(T. ,

)

'-'- (4. it) 3T T. - T. , 1 1-1

th where Tj^ is the temperature after the i^'' iteration which is obtained from

A(T.)

3T

The quasi-Newton iteration exhibits superlinear convergence, and requires only one evaluation

of A per iteration (after the first). The determinant computation is performed by subroutines

DSPFA and DSPDI from the LINPACK software (Dongarra et al . [2]). The first routine factors a symmetric indefinite matrix (stored in packed form) and the second computes the determinant as

a product of the pivots.

Once the iteration to find the stability limit has converged, the cubic form is evaluated: the first step in this process is to determine the critical displacement vector Z; i.e., solve

Q Z = . (4.6)

In general this step is efficiently performed by the inverse iteration technique (Peters and

Wilkinson [11]). In the implementation described here, Z is obtained from the LINPACK

subroutine DSPCO.

Because Z is determined only within an arbitrary multiplicative constant, it must be

subjected to further constraints to avoid arbitrary changes in the sign or magnitude of the

cubic form evaluated along the stability limit; such changes falsely indicate the location of

critical points. The "length" of Z is normalized by the constraint

z"^ Z = 1 . (4.7)

16

I
We determine the direction of !_ by specifying that its first component Snj be positive at the

initial guess of the volume. Constraining (>n^ to be positive to fix the direction of !_ cannot

be continued as we move along the stability limit searching for the zero of the cubic form.

Because it is possible for 6ni to go through a zero as Z_ changes with V along the stability

limit, constraining 6ni to always be positive would be equivalent to changing the direction of

Z at a value of the volume where 6ni goes through a zero. The technique used to ensure that

the direction of Z_ varies continuously as a function of V along the stability limit is the

following. The null vector from the previous iteration on V is saved and projected onto the

null vector computed for the current iteration. If the projection is positive we assume that

the direction of Z_ is changing continuously. This final condition on Z_ may be written

z'[_^ Z . > , (4.8)

where Zj^ is the null vector from the i^'^ iteration on V. If Z_]_ does not satisfy C^.S), it is

multiplied by -1 so that it does.

Michelsen and Heidemann [7] have recommended that discontinuities in Z be avoided by

making use of the condition

z"^ Z. > , (1.9)

where Iq corresponds to the initial guess Vq. Depending on how far Vq is from the critical

point, it is quite possible that Z could change enough along the stability limit that

eq (M.9) would result in an arbitrary change in direction. It seems safer, and no more

costly, to employ eq (4.8).

The summation of the cubic form can in general be reduced from c3 terms to c(c+l ) (c+2)/6

by making use of its symmetry properties, i.e.,

c c c c 6'A = I A.., Z.Z.Z, = y y y h.., a.., Z.Z.Z, (4.10)

.

, , ijk 1 J k .... . , ijk ijk 1 J k i,J,k = l ^ ^ k=j j = i 1 = 1

-^ -^ ^

where

1; i=j=k

hijk = 3; i=j^=k, or i = k|=j, or i=j=j=k

6; i=t=j=fk and i=|=k For the special case of the two-constant cubic equation of state to be considered here, the

cubic form has been simplified so that only a double summation is required for its evaluation.

Michelsen [6] suggested that the cubic form could be evaluated numerically because it is

the derivative of the quadratic form in the direction of the critical displacement vector;

hence,

- z'^Q(T,V,n + £Z)Z - z'^Q(T,V,n)Z

6^A S C = TT—ZTT (4.11) ¥

where n is a vector of the mole numbers, i.e., (n-], n2 Hq)',

e is a small number, i.e., 0(10~6)^

zT Q(T,V,n)Z = on stability limit,

lleZjl = e||z)| = £, since ||z|| = Z^Z = 1.

17
This expression for the cubic form is particularly useful when the required derivatives of A are not expressible analytically.

Once the stability limit is located and the cubic form is evaluated at the current value,

Vj^, the volume is updated by

C(V.)

^i.1 =^i -C^V ' ^'^•^2)

where C'(Vj^) represents the derivative of C with respect to the volume (dC/dV) or some estimate of it. The original method of Heidemann and Khalil approximated dC/dV numerically,

i.e.,

d C(V.) C (V. + 6V) - C (V.)

dV 6V ^^'':>J

where 6V = lO'^Vj^. This approximation is expensive because, when the volume is incremented,

eqs (^.1) and (M.6) must be solved again before evaluating C{y^ + 6V). A good approximation

can be made for T(Vj^ + 6V) with only an inexpensive calculation (see below) so that only one

or two iterations are required to solve (4.1). But we must keep in mind that the single most expensive calculation in this method is the factorization of Q; thus we should try

to minimize the number of evaluations of the cubic form. Michelsen and Heidemann [7]

presented a quasi-Newton method which requires only a single evaluation of C for each

iteration on V. In general, they used

C(V.) - C(V. ) C'(V.) = y \ y ^-^ . (4.14)

i i-1

Two values of C are needed to start this iteration. One option Is to use (4.13) for the first

iteration, then switch to (4.14) for subsequent iterations. Michelsen and Heidemann chose to use a modified form of (4.13) for the first iteration. To evaluate C(Vj^+6V) they make two approximations: first, that Z^i = Zq, and second, that T (Vj^ + 6V) is obtained from a Taylor

series expansion, i.e.,

T(V. + 6V) = T(V.) + SV IJ I , (4.15)

i

where T(Vj^+6V) is the first order approximation to the temperature on the stability limit at

(Vj^+6V), T(Vj^) is the temperature on the stability limit at Vj^, and dT/dV is the derivative

along the stability limit. To obtain an expression for this derivative, consider the

stability limit for a given mixture (fix n) to be parameterized in terms of V. Since the

stability limit is defined by Z^Q Z^ = 0, then

ri T '^^ T '^S T

'^- dv^i5^^=dV 5Z^Z d^^^^Q-=0

(4.16)

T dQ

= ^ d?^=°

18
Since Q Z = Z^Q = 0. By chain rule,

dQ 3Q 3Q .„ dV 9V 3T dV

^^•' 1

1

Substitution of (^.17) into (M.16) gives

dV - 3Q • ^''•^^^ ^ 3f^

These two approximations allow an estimate of C(Vj^+6V) which requires no additional

factorizations of Q. After the initial iteration on V, C(Vi) in (M.14) is evaluated exactly,

and the initial guess of temperature at a new volume V^ is made analogously to (4.15), i.e.,

dT

T (V.) = t[v. J + (v. - V. J -Si

1 ^ 1- V ^1 1-1 dV

(i|.19)

^i-1

Although the quasi-Newton method requires only half the number of evaluations of the

cubic form per iteration on V as does the full Newton-Raphson technique described initially,

its convergence is slower (superlinear rather than quadratic). In general, the cost of

computing critical points by various methods will depend on the cost of computing the

derivative C'(Vj^) relative to the convergence rate which can be attained with that

approximation. This must be determined by numerical experiments.

4.1 Analytical Expression for the Volume Derivative of the Cubic Form

It may be advantageous to retain the quadratic convergence of a Newton-Raphson iteration

if the derivative of the cubic form with respect to volume can be computed more efficiently.

In this section is presented an analytical derivation of dC/dV which requires one additional

factorization of an augmented Q matrix to be computed for each volume iteration. The

procedure is substantially more efficient than the forward difference approximation of dC/dV

(eq 4.13) which requires two evaluations of C at each volume iteration.

We begin, as previously for dT/dV, by considering the stability limit to be parameterized

in terms of V. Then

c dA

.

dZ

.

^= I [—

^

Z-Z-Z, +3A.. -;^ Z.Z, ] (4.20)

dV . .

, , dV 1 J k ijk dV j k-' i,j,k=l

where

dA . .

, 3A .

.

,

ijk _ ijk

dV (4.21)

The partials OA^j^-ZSV) and (3Aij|</3T) may be obtained by differentiation of the expression

dZ for Aj^ji^ obtained from the equation of state. The derivative — is obtained from the

following considerations. Differentiating the expression Q Z, which is zero on the stability

limit, leads to

19
dZdV

dQ

d^^ (4.22)

dZ This expression could be solved for — except that Q is singular on the stability limit.

dZ Therefore we must find another condition on ^rr which will be used to replace the row in Q

dV = containing the zero pivot. The appropriate condition is from a simple result of vector

calculus: the derivative of a vector of constant length is orthogonal to that vector, i.e.,

- dV (4.23)

This condition can be combined with (4.22) and the result written as the solution of c+1 equations in c unknowns,

' " dQ

Q dZ dV^

vT' dV~

= ______

_ L

(4.24)

5. Implementation for a Two-Constant Cubic Equation of State

Expressions for the derivatives of the Helmholtz free energy are required to evaluate A,

C, and their derivatives with respect to T and V. These expressions are given here for a general two-constant cubic equation of state (TCC EOS) following, for the most part, the

notation of Michelson and Heidemann [7].

The TCC EOS is written

NRT V-NbN^a

(V + 6 Nb)(V + &^Hb) ' (5.1)

where

—rr ) n.n .a. . n2 i.L ' J ^j

(5.2)

1=1 (5.3)

N =

I n. , 1=1

a.. = ^.. (a., a..)

1/2

(5.4)

a. . = n R^T^./P . fl + m.[l

11 a ci ci ' 1 /2i i2 (T/T .)''^]

ci (5.5)

m . = m(a). ) , 1 1

20
b . = n^ RT . /P . 1 b ci ci (5.6)

For the SRK equation of state,

6, = 1, 6„ = 0, Q = 0.M2748, n = 0.0866iJ , (5.7)

m(u.) = 0.480 + 1.574u. - 0.176a)/ (5.8)

and for the PR equation of state,

6, = 1 + /2, 6, = 1 - /2, n = 0,45724, Q = 0.07780 , (5.9)

I ^ cL D

m(a).) = 0.37464 + 1.54226 w. - 0.26992 w.' (5.10)

Heidemann and Khalil chose to proceed from the expression for fugacity

RT In f.= -

1

- 3P Sn.

1

RTn "i ^^ ^] dV + RT in ^;j—

T,V,n

(5.11)

Their choice of (n-| , n2. ..., nj,) as the independent variables results in all derivatives

containing differentiations with respect to a mole number. However, if one of the mole

numbers in the independent variable set were to be replaced with the volume, then one element

of g (i.e., —-) would contain no differentiations with respect to mole number. In that case,

3V those derivatives would have to obtained by differentiating the Helmholtz free energy rather

than the fugacity. To be consistent with any derivatives obtained from the fugacity, the free

energy must be relative to a reference state of pure ideal gases at the system temperature and

unit pressure, i.e. (see Prausnitz [12], p. 40),

A - A

CO 1=1 1 (5.12)

where

^°= I

i = 1

A. = Helmholtz free energy of pure i in an ideal gas state at the system T, and p = 1

When this expression is differentiated with respect to mole numbers,

3 A.

: 3 n. y? - RT (5.13)

T,V

21
where y = the chemical potential of pure i in an ideal gas state at the system

i T, and p = 1

.

The Helmholtz free energy for the TCC EOS is

c n.RTF

A-A° = RT l^ n. in^j^-^F^ , (5.14)

where Fj^ are functions of the dimensionless volume K {= V/Nb) which are summarized at the end

of this section. Differentiating (5.14) with respect to mole numbers gives

n.RTF e. RT in f. = RT Un -\^ ^ B. F^ ] - |

(,. p^ . _i F^] , (5.15)

where

c n . a. .

j=i

b.

The elements of Q are obtained from further differentiation of (5.15) with respect to mole

numbers , i.e.,

9 In f

.

N 6. . f^^ij = ^^^^-37^ —^T,V, n, = ^^t^^ ^ ^^i ^ 2j) ^1 ^ ^i^j ^1 J J k 1 (5.18)

where 6j_j is the Ki^onecker delta, and n^. held constant in the differentiation signifies that

all mole numbers are held constant except nj

.

To evaluate the cubic form one further differentiation of (5.18) with respect to the mole

numbers is required, i.e.,

? ^ 3^ in f ' ^Jk = " «^ (3?r-3-n-^T,V, n^

= RT[- ii-iii . (6.B. . S.B^ ^ 6.6,)F^2 ^ ^ e.B.B, F^^]

' (5.19)

-f {[2 (a.6.6, ^ u.B.B^ . a^B.B.) - 3 ^,^^\](F^^ F5)

- 2 B.B.g, F,, -- (a. .6, + a., 6. + a.,6.)F,} , 1 J k 4 a ij k Ik J jk 1 6'

'When this expression is substituted into the cubic form, algebraic simplification yields

2 2 ^ ?

c 9^ In f. c N Z-;* _ _ _

N RT I ^-g-i Z.7 Z^ = RT [ -

I -^ . 3 Z (6 F^)^ . 2(6 F^)^]

i,j,k=1 J k -^ 1 = 1 n.

^ (5.20)

22
+ ^

[(6a B^ - 3B^)(F2 + F^) - 2^'^?^ - 3 a 6 F^] ,

where

Z s ;_ Z. , (5.21)

i

a = I a.Z. , (5.22)

i

= ^ B. Z. , (5.23)

i

a s - A y a. . Z. Z. . (5.24)

a

^ j

ij 1 J

(There is a typographic error in the last term on the right side of eq (18) in the paper by

Michelsen and Heidemann [7].) To compute dT/dV along the stability limit, we need the partial

derivatives of qy with respect to T and V; these are given by

3 q 3F 8F N-^. RT [(6. . ,.) —.2 B.Bj F, ^]

(5.25)

9F a. . 3F 8F,

+ |-b'-ij3K [b.B.^ adK — ^ " (B.6. -a. 6. -a.BJxTT-] , ij ij J13K-'

where the derivatives of the Fj^ with respect to the dimensionless volume K are given at the

end of this section, and

3 q. . N 6. . „ B.8.

N -.7:;^^ = R i ^ + (6. + B.) F, + B.B.F, 1 + 4| -^ (F, + F,)

3T ^ n. '"^i ^^j' 1 '^i'^j 1 ^ 3T b ' 3

6'

3 a. , F^ 3(aa.) B . 3(aa.) B. JJ _5 _ r i_ _A + J_ _il F

(5.26)

9Tb '3Tb 3T b * 6

where

3 a.. (R T .)^ n. 1/2 - m.

^^ ^ P ^

' ^ci (T T .)^^^

ci ci

9 a,. 1/2 C. • 3 a.. 9 a.. M '-^-TT^ [^^ a,, . a,, ^r=^] (5.28)

3T ^ 1/2 ^ 3T jj "ii 9T

ii jj

9 (aa ) c 9a..

-9f^ = N .^, "j -9T^ (5.29)

23
^ ' ? .,i,

"'"j

^ • ''•'°'

The quadratic forms required to evaluate dT/dV can be simplified to: NZ^^Z = RT (2iZF^^ W B^p^ F^^

)

^f t^' ^K3 -^^K5 ^ ^^' ~2 a B) F^^^:

and

T

^? A - - -2 2 NZ-^Z = rUI — -^SBZF^ + B F^]

i i

b t^T ^ (^3 ' h^ - ^T ^ - 2 "t ^ ^6^ •

;5.3i)

(5.32)

where

8F.

^T " If ' ^^-^^^

c 5(aa.)

a„ = I .T.

^ Z , (5.35)

T .^^ 3T 1

c „ ^T^i,l=1 -ir^^i^J • ^5.36)

Dividing (5.31) by (5.32) gives dT/dK on the stability limit which is different from dT/dV by

a factor of Nb. To obtain the derivative of the cubic form with respect to K along the stability limit,

differentiating C (5.20) directly gives

2 3 ? dr c N Zf _ N^ ^ = RT^ [ -

I —^ + 3 Z (6 F^)^ + 2 (6 F^)^]

i=1 n

.

1

2 2

c " Z Z _ _ __

^ 3RT [ - I

2 ^ ^ ^^V ^ ^^^ ^K^^l ^ ^^^

""l^Kl i=l n

+ 2 6^ 6^ F^ ^ 28^ F^ F^^ ]

+ ^ 1(6 aai^ - 3a8^)(F^2 * ^k6^^ [^((aa)^^ 2^ + 2 aaSBj^] (5.37)

24
where

-3 ia^^^ + 3a6^6j^)] (F^ + F^)

-2 [(aj^e^ + 3ae^ Bj^)F^ + ag^ F^^^;

-3 [l(aa)^6 + ^^^K^^e * ^^^ ^K6 ^

^ '

T, =-f. (5.38)

dZ. ^^ dZ

^Ki^ diT ' ^ component of — , (5.39)

Z, -f = j Z^. . (5.n0)

1 = 1

1 = 1

^ da _ ^ dT re; ils>^ ^K dK 3T dK ' ^b.^^;

1=1

d(aa.) 3(aa.) ,_, }^_ _ i_ oT dK " 9T dK '

da.. a. . ._ y. ^ iJ dT

dK 9T dK

(5.iJ'4;

(5.46)

When tracing critical lines in regions where no data are available to indicate whether or not critical points actually exist, it is useful to be able to compute the mathematical

condition for stability, i.e., that the quartic form is positive. The quartic form for the

TCC EOS was evaluated by taking the directional derivative of the cubic form along the

critical displacement vector; that is, apply the operator ^ " ? Z. ,4- (5.47)

dn . , 1 dn

.

1 = 1 1

25
to eq (5.20). The result is

3 4 c N^ Z

i=1 r\i1

+ I

[ 12 (a6^ - 2 ai^ + S^)(?2 + F^) (5.48)

- 8 (2 aB^ - 6^) Fj^ + 5 6^ F^] .

5.1 Functions of the Dimensionless Volume K and Their First Order Derivatives

Let 'A =

6i - 52, then

F^ = 1/(K-1) ,

2 A "-K+e^ K+5^'' '

R = 1 r f_J_i^ -f ^-1^1

3 A

'- 'K+S '' ''K+5 '' ^ '

F - 1 r r-!L_i3 -f ^2 >3i

i) A

' ''K+5

' 'K+e '' J '

9 A

'- 'K+5 '' ^K+S^ J '

F^ = 4 In

'^ 5 A 'K+&^'

^6 = ^2 - ^5 •

dF. and their first order derivatives -rrr- = F,,. are

dK Ki

F = - F

2

KI "l '

K2 A (K+6^)^ (K+S^)

2 2

F = 1 [ - 1 . ?-^] , K3 A (K+6^)^ {K*6^)^

26
X4

H 4-

K5 (K-H6^)(K+5^)

F , = F K5 K2 :<5

•

5. How to Use the Critical Point Software

Two subroutines have Deen written to compute critical points. CUBCR3 uses a quasi-Newton volume iteration and CUBCRPT uses a Newton-Haphson volume iteration in which

the volume derivative of the cubic for:^ is evaluated analytically (eq 5.37). The calling

arguments are the same for either of these subroutines, i.e.,

CALL CUBCRPT (NTYPE, NC, MXC, X!i, T, ?, V, IPR, CONV),

On entry,

NTYPE

NC

MXC(I)

XN(I)

T

IPR

On return,

T

PVCONV

specifies type of TCC EOS. Current options are:

1 for SHK EDS (Soave [13])

2 for PR £03 (Peng and Robinson [10])

# of mixture components

identification number for each mixture component. The identification number

corresponds the position of the component's data in the arrays in subroutine

CUBDAT. Check a source listing to determine the pure components currently

available.

mole numbers (n, n2, ..., n^).

initial guess of critical temperature (K). If set to 0.0 an initial guess is

made internally.

initial guess of critical volume (em3), if set to 0.0 an initial guess is made

internally.

print flag for debugging. Set to for no output. Set to 1 for summary of

iterations written to unit 5.

critical temperature (K)

critical pressure (bar)

critical volume (cm3)

convergence flag (logical).

= .TRUE, if convergence achieved

= .FALSE, if iteration fails

Externals,

CUBPROP: CUBDAT, FVOL, STABLIM, CUBMXT , ZSUM, CFORM, S: CUBMX, CUBP.

LAS: SUMVC, DOT

FORTRAN: DABS, DSIGN

27

:k, zksum, dck.
Usage: The variables XN, T, P, and V are double precision. The subroutines which must

be linked to the driver are located in the following files: CUBCRIT, CUBPROP,

BIP, DSPMAT, and LAS. The software is structured so that all the pure component data that are required to compute a mixture critical point are contained in subroutine CUBDAT, This scheme allows the

user to simply specify a set of identification numbers for the mixture components, and the

number of moles, of each component in order to compute a critical point. However, to determine the critical point of a mixture containing a component (or components) for which the

pure component data are not contained in CUBDAT, one must modify the software to include the

required data. The pure component data which are required for each mixture component are: a critical temperature (K) , critical pressure (bar), and the acentric factor. The arrays in CUBDAT should be dimensioned as follows: PARAMETER (N=# of pure components whose data are stored in CUBDAT, NX=N*(N+1 )/2)

DOUBLE PRECISION PC(N) , TTC(N) , W(N) , XI(NX)

In addition to the pure component data, any binary interaction coefficients (Cij) which

are not 1.0 are set in the subroutines SRKBIP (for the SRK EOS) or PRBIP (for the PR EOS). Just follow the example of the values already set there.

Storage for all arrays which may be variably dimensioned, depending on the number of

mixture components, is set in subroutine CUBCRPT or CUBCR3. The proper dimensions are the

following:

PARA^4ETER (N=# of mixture components, N1=N+1, NX=N* N1/2)

INTEGER KPVT(N)

DOUBLE PRECISION AIJ(NX), AIJT(NX), AIT(N), ALI(N), ATIJ(NX), AW(N), BI(N), BEI(N),

Q(NX), QA(N1,N), TCA(N), VCA(N), Z(N), ZK(N1), ZZ(N)

6.1 Verification of the Software

The accuracy of the software was verified by a twofold procedure:

1) the results of a binary critical point calculation were compared to results obtained

from an independent code which solved the traditional critical point equations of

Gibbs;

2) the software was checked for internal consistency by comparing selected values of

derivatives which are computed using the formulas given in this note, with the values

obtained using different formulas.

Calculation of binary critical points was done for the nCi5-C02 mixture using the Peng

and Robinson equation of state. A detailed consideration of the phase diagram for nCi5-C02 as predicted by the PR EOS, along with formulas for calculating critical points, can be found in Hong [16].

Figure 6.1 shows the result of a complete tracing of the critical line using a binary

interaction coefficient <; = 0.919, and the following pure component data:

nCi6 CO2

Tc(K) 717.0 30^.21

Pc(MPa) l.iJ2 7.38M

w 0.746 0.225

28
40

30 -

S.

(D

D

(0

CO

20 -

10 -

Temperature, K

Figure 6.1 Calculated critical locus for the hexadecane + carbon dioxide system

29
The gas constant was set to 0.00831'J3'< MPa'L/(mol -K) . A benchmark result for x(nC-|5) = 0.99

follows:

Tc = 716.701292254 K Vc = 1 .27759526809 L/raol Pc = 1 .i46980319713 MPa

The independent calculations both produced the above result to the accuracy presented.

The following consistency checks were performed on multi-component mixtures containing up

to 12 components:

1) Values of the cubic form computed using eq (4.10) were compared to values obtained

from the "numerical" version, eq (4.11);

2) Values of dT/dV computed using eq (4,18) were compared to values obtained from a difference formula;

3) Values of dZ/dV obtained by solving the system of eq (4.24) were compared to values

obtained from a difference formula;

4) Values of dC/dV computed using eq (4.20) were compared to values obtained from a difference formula.

Generally, 6 to 8 significant figures of agreement were obtained from these comparisons.

7. Summary of FORTRAN Subprograms

1 ) CUBDAT computes the temperature independent part of the pure component parameters and

mixture cross terms. The calling sequence is

CALL CUBDAT (NTYPE, NC, MXC , G, TCA, VCA, AIJ, AW, HI, R).

On entry, NTYPE specifies type of TCC- EOS. Current options are: 1 for an SRK EOS, 2 for a PR EOS. NC # of mixture components.

MXC(I) identification # of each mixture component.

On return,

G(I) G(1) = 6i, G(2) = 62 TCA(I) critical temperatures (K) of pure mixture components

VCA(I) critical molar volumes (cm3/mol) of pure mixture components

AIJ(IP) upper triangular part of the symmetric matrix of a^^j packed by columns, i.e.,

IP = i + j(j-1)/2. In unpacked notation,

AlJd.I) = fiaR^T^i/Pci

AIJ(I,J) = Cij[AIJ(I,I) » AIJ(J,J)]l/2

AW(I) = m(a)i)

Bid) = % RTci/Pci

R gas constant = 0.00831434 MPa'L/(mol -K)

Externals,

BIP: 3RKBIP, PRBIP

FORTRAN: DSQRT

30
Usage: CUBDAT is only called once for any set of critical point or stability limit

calculations on a given mixture since it computes quantities which are independent of (T, V,

rii).

2) CUBMX computes the temperature and composition dependent coefficients. The calling

sequence is CALL CUBMX (NC, T, XN, XNT, TCA, AIJ, AW, BI, ATIJ, ALI, BEI, A,B).

On entry,

NC # of mixture components (c)

T temperature (T) in kelvins

XN(I) mole numbers (n-| , n2, ..., n,-,)

c XNT N = y n.

i = 1 ' TCA(I) Tci

AIJ(IP) temperature independent part of ay (=a!.)

AW(I) m (wi) BI bi

On return,

ATIJ(IP) aj^j in packed storage of upper triangular part

ALKI) a-- = TT- I n. a. . 1 Na .^^ J ij

BEKI) 6i = bi/b

1 a = — ? r 2, n.n .a. . N^ i,j=1 ' J ^J

1 ? B b = ^ I n.b.

Externals,

FORTRAN: DSQRT

Usage: CUBDAT must be called before CUBMX.

3) CUBMXT computes the temperature derivatives of the EOS coefficients. The calling

sequence is CALL CUBMXT (NC, T, XN , XNT, TCA, AIJ, AW, ATIJ, AT, AIT, AIJT).

On entry,

NC c T T

31
XN(I) rii XNT N TCA(I) Tci

AIJ(IP) a'ij

AW(I) m(a)i)

ATIJ(IP) aij

On return,

c da..

d (aa.

)

c da..

^"^^) -If-^ =

Jf .1 "j -dT

da.. a.. , da.. . da..

dT 2 a . . dT a . . dT

11 JJ

where

da.. -.1/2 -m.

11 t Ii fi ( T ) 1 1 r 1 = a! . {1 +m. [1 - dT

^^ ^ T . (T T .)^^^

ci ci

upper triangular part of symmetric matrix is stored in packed form.

Externals,

FORTRAN : DSQRT

Usage: CUBDAT and CUBMX must be called before CUBMXT.

4) FVOL computes the functions of the dimensionless volume (V/Nb) and their 1st order

derivatives. The calling sequence is

CALL FVOL (V, XNT, B, G, FV)

On entry,

V system volume (V) in cm3 XNT N

B b G(I) G(l) = 6i, G(2) = 62

On return,

FV(I) FV(1) = Fi, FV(2)=F2, FV(3)=F3, FV(i|)=Fij, FV(5)=F5, FV(6)=F5, FV(7)=Fki,

FV(8)=Fk2, FV(9)=Fk3, FV(10)=FKij, FV(11)=Fk5, FV(12)=Fk6. FV(13)=F9

Externals,

FORTRAN: DLOG

32
Usage: CUF3DAT must be called before FVOL. The argument B may be supplied either by a call to CUBMX or simply by summing:

b = - I n.b.

N . , 1 1 1 = 1

5) CUBP is a double precision function which returns the pressure. The calling sequence

is

P s CUBP (T, VM, G, A, B, R)

On entry,

T T (K) VM V/N (cm3/mol)

G(I) 6i , 62 A a

B b

R gas constant

On return,

CUBP system pressure (P) in bar

Usage: CUBDAT and CUBMX must be called before CUBP.

6) CUBDPV is a double precision function v;hich returns TTT? j-r m

DPV = CUBDPV (T, VM, G, A, B, R)

The argument list is identical to that for CUBP.

7) CUBF computes fugacity. The calling sequence is CALL CUBF (I, T, V, XN , XNT, G, A, B, ALI, BEI, F, R)

On entry,

I index for mixture component i T T

V V XN(I) "i XNT N G(I) 61 , 62 A a B b ALI(I) ai BEI(I) 6i R gas constant R

On return.

F RT in fi

33
Externals,

FORTRAN: DLOG

Usage: CUBDAT and CUBMX must be called before CUBF

8) QIJ computes matrix element Qij. The calling sequence is CALL QIJ (I, J, T, FV, XN, XNT, ATIJ, A, B, ALl, BEl, Q, R)

On entry,

I index for mixture component i

J index for mixture component j T T FV(I) Fi, FKi

XN(I) n^ XNT N ATIJ(IP) aj^j (packed storage)

A a

B b ALI(I) ai BEI(I) 6i

R gas constant

On return,

8 In f. Q RT —^—-

J

Usage: CUBDAT, CUBMX, and FVOL must be called before QIJ. This routine is called by PAKQ to assemble Q.

9) PAKQ assembles the upper triangular part of the symmetric Q matrix into packed form,

the calling sequence is CALL PAKQ (NC, T, FV, XN , XNT, Q, ATIJ, A, B, ALI, BEI, R)

On entry,

NC # of mixture components

T T FV(I) Fi and Fki

XIJ(I) "i XNT N ATIJ(IP) aij

A a

B b ALIO ai 3EI(I) 6i

R gas constant

3>i
On return,

Q(IP) packed ^

Externals,

CUBPROP: QIJ

Usage: CUBDAT, CUBMX, and FVOL must be called before PAKQ. PAKQ is called by DETQ which

computes A, and by DZK which computes Zy^.

10) DETQ assembles g by calling PAKQ, then computes its determinant using the LINPACK

routines DSPFA and DSPDI. A flag can be set to indicate that the null vector is to be

computed, in which case DSPCO is called instead of DSPFA. The null vector is normalized to unit length. The calling sequence is

CALL DETQ (NC, T, FV , XN , XNT, Q, KPVT, Z, INL, RCOND, DTQ, ATIJ, A, B, ALI, BEI, R)

On entry.

NC // of mixture components

T T FV(I) Fi and FKi

XN(I) "i XNT N INL null vector computation

ATIJ(IP)

A

BALI(l)

BEI(I)

R

INL = for no null vector computation

INL = 1 for null vector computation

When INL is set to 1 , it is assumed that the £ matrix has already been

determined to be singular.

'ij

gas constant

On return,

Q(IP)

KPVT(I)

Z(I)

RCOND

DTQ

information from the factorization of £ by DSPFA or DSPCO

array used by LINPACK routines to store pivot information

null vector components Z^^. Z(I) is computed by LINPACK routine DSPCO, and

contains an approximate null vector when £ is nearly singular, i.e., when the

reciprocal condition number is on the order of machine roundoff.

reciprocal condition number computed by DSPCO

A=det I a I computed by DSPDI

Externals,

CUBPROP: PAKQ

LINPACK: DSPCO, DSPFA, DSPDI

35
FORTRAN: I DINT, DSQRT

Usage: CUBDAT, CUBMX, and FVOL must be called before DETQ. This routine is called by QNT to find T on the stability limit, then compute !_.

11) QNT performs the quasi-Newton Iteration on temperature to solve A(T) s det |£| = for a given volume. The calling sequence is CALL QNT (NC, T, FV, XN , XNT, RDT, TOLT, ITMX, IPR, IT, IFT, Q, KPVT, Z, TCA, AIJ, AW, ATIJ,

EI, ALI, BEI, R)

On entry,

NC

TFV(I)

XN(I) XNT

RDT TOLT

ITMX IPR

TCA(I)

AIJ(IP)

AW(I)

Bid)

R

# of mixture components

initial guess of TCTq)

N6T/T where 6T is the increment used to compute the derivative 3A(TQ)/tJT

iteration convergence tolerance

I Ti + i - Ti |/Ti < TOLT

maximum number of iterations

print flag:

for no printing

1 for printing of iteration progress on unit 6 Tci

m(wi

)

bigas constant

On return,

T

ITIFT

temperature on stability limit (A=0)

number of iterations to convergence

iteration flag:

successful iteration

1 iteration failed

Q(IP) factored £ (from DSPFA)

KPVT(I) pivots (LINPACK routines)

Z(I) not used

ATIJ(IP) aj^j at T on stability limit

ALI(I) tti at T on stability limit

BEKI) Si

Externals,

CUBPROP: CUBMX, DETQ

FORTRAN: DABS, DSIGN

36
Usage: CUBDAT and FVOL must be called before QNT. CUBMX and DETQ are called from within the

iteration loop on temperature in QNT. This routine is called by STABLIM.

12) STABLIM locates the temperature on the stability limit when the volume is fixed.

Then the null vector is computed. On the 1st call to STABLIM the direction of Z^ is determined by requiring Z-j > 0. On subsequent calls to STABLIM, for a given critical point or

T stability limit calculation, the direction of Z_. is determined by requiring Z._ Z. > 0.

The calling sequence is

CALL STABLIM (NC, T, FV, XN , XNT, RDT, TOLT, ITMX, IPR, IT, IFT, IFC, Q, KPVT, Z, ZZ, TCA,

AIJ, AW, ATIJ, BI, ALI, BEI, A, B, R)

On entry,

NC

TFV(I)

XN(I) XNT

RDT TOLT

ITMX

IPR

IFC

TCA (I)

AIJ(IP)

AW(I)

BI(I)

R

# of mixture components

initial guess of T on stability limit

N6T/T, where 6T is the increment used to compute the derivative 3A(Tq)/8T

iteration convergence tolerance, i.e.,

I Ti + i - Til /Ti < TOLT

maximum number of iterations

iteration print flag:

for no printing

1 for printing of iteration progress on unit 6 first call flag

1 first call; set Zj > T not first call; require Z. , Z. > ^ -1-1 -1

Tci

m((i)j^

)

bigas constant

On return,

TITIFT

Q(IP)

KPVT(I)

Z(I)

ZZ(I)

ATIJ(IP)

ALI(I)

temperature on stability limit

# of iterations for convergence

iteration convergence flag--not currently used--if iteration fails in QNT,

execution is stopped in STABLIM.

factored g (from DSPCO)

pivot information (from DSPCO)

li

Z^j^-I (from previous iteration on volume)

37
BEKI) Bi A a B b

Externals,

CUBPROP: QNT, CUBMX, DETQ

FORTRAN: DABS

Usage: CUBDAT and FVOL must be called before STABLIM.

13) ZSUM computes the summations involving the null vector components which are used to compute the quadratic, cubic, and quartic forms. The calling sequence is CALL ZSUM (NC, XN, XNT, Z, ATIJ, ALI, BEI, A, AIT, AIJT, ZS)

On entry,

NC # of mixture components (c)

XN(I) ni XNT N Z(I) null vector components (Z-| , Z2 1q)

ATIJ(IP) aj^j (packed storage)

ALI(I) tti BEI(I) 6i A a AIT(I) d(aai)/dT

AIJT(IP) d(aij)/dT (packed storage)

On return.

2

c Z. n. ZS(I) ZS(1) = I ^^ wher- " ^

..,.^, e y ^

1=1 '1

c Z.3

ZS(2) = I -ir

c Z.'* ZS(9) = I ^

c zs(3) = z = I z.

i=1 ^

c ZS(4) =6=1 B.Z

i=1

38
ZS(5) = a = I a Z

i = 1

- 1 % ZS(6) = a = - y a, , Z.Z.

c d (aa.

ZS(7) = a, = I —Ht^Z.

1=1

ZS(8) =l^= l^ Z.Z

i.J = l

Usage: CUBDAT, FVOL, STABLIM, and CUBMXT must be called before ZSUM.

1^) QFORM computes the quadratic form. The calling sequence is CALL QFORM (XNT, R, T, A, B, ZS, FV, QF)

On entry, XNT N

R gas constant

T temperature on stability limit

A a

B b ZS(I) sums involving Z FV(I) Fi and F^^

On return,

QF zT g Z Usage: CUBDAT, FVOL, STABLIM, CUBMXT, and ZSUM must be called before QFORM.

15) CFORM computes the cubic form. The calling sequence is

CALL CFORM (XNT, R, T, A, B, ZS, FV, CF)

On entry,

Same as QFORM

On return,

? a^A

OF I. I ^, Z.Z.Z,

. .

, , 3n. 3n . 3n, i J k

Usage: same as QFORM

16) DFORM computes the quartic form. The calling sequence is

39
CALL DFORM (XNT, R, T, A, B, ZS, FV, DF)

On entry, same as QFORM.

On return,

c 4 DF I i:—, ^/ . - Z.Z.Z, Z^

. .

, -, , 9n. 3n .dn, dn, i j k 1 i,j,k,l=l 1 J k 1

Usage: same as QFORM.

17) SLDTK computes dT/dK along the stability limit. The calling sequence is

CALL SLDTK (R, T, A, B, AT, ZS, FV, DTK)

On entry,

R gas constant

T T

A a

B b AT da/dT

ZS(I) sums involving Z FV(I) Fi and Fki

On return, DTK dT/dK

Usage: CUBDAT, FVOL, STABLIM, CUBMXT, and ZSUM are called before SLDTK.

18) QKMAT assembles the upper triangular part of the symmetric matrix d£/dK in packed

form. The calling sequence is CALL QKMAT (NC, T, FV, XN , XNT, QK , ATIJ, A, B, ALI, BEI, AT, AIT, AIJT, DTK, R)

On entry,

NC # of mixture components

T temperature on stability limit

FV(1)

XN(I) n^

Xt>lT N ATIJ(IP) a^j

A a B b ALI(I) ai BEKI) h

AT da/dT

AIT(I) d(aai)/dT

AIJT(IP) d Sij/dT

DTK dT/dK

i»0
R gas constant

On return,

QK(IP) d£/dK = 8g/3K + 0£/8T) (dT/dK)

Usage: CUBDAT, FVOL, STABLIM, CUBMXT, ZSUM, and SLDTK must be called before QKMAT

19) DZK computes dZ^/dK along the stability limit. In assembling the rectangular matrix,

the bottom row (Z^"^) is scaled so that the largest component of Z is the same magnitude as the

largest component of Q. The rectangular system is solved using a simple Gaussian elimination

algorithm with partial pivoting, i.e, row exchanges only. Thus, at some point in the

elimination sequence, a zero pivot is encountered (because £ is singular) and the bottom row will be exchanged with the row containing that pivot. The elimination sequence is completed,

and the back substitution proceeds as if the system were square. The calling sequence is CALL DZK (NC, T, FV, XN, XNT, Q, ATIJ, A, B, ALI, BEI, AT, AIT, AIJT, DTK, NCI, MQA, QA, Z, ZK, R)

On entry,

NC

TFV(I)

XK(I) XNT ATIJ(IP)

A

BALI(I)

BEI(I)

ATAIT(I)

AIJT(IP) DTK

NCI MQA

IJ

# of mixture components (c)

temperature on stability limit

F

n

N

a

a

b

da/dT

d(aui)/dT

d a^j/dT

dT/dK

c + 1 1st dimension of QA array

QA(I,J) work space used for rectangular matrix. Must be dimensioned at least QA (NCI,

NC)

Z(I) Zi

R gas constant

On return,

Q(IP)

ZK(I)

Externals,

CUBPROP

LAS:

PAKQ, QKMAT

MSAX, RMSOL

41
FORTRAN: DABS

Usage: CUBDAT, FVOL, STABLIM, CUBMXT, ZSUM, and SLDTK must be called before DZK.

20) ZKSUM computes the sums involving Zj^j^ . The calling sequence is

CALL ZKSUM (NC, XN , XNT, Z, ZK, ATIJ , A, ALI, BEI, ZS, DTK, ZKS)

On entry,

NC # of mixture components

XN(I) ni XNT N Z(I) Zi

ZK(I) Zki

ATIJ(IP) aij

A a ALI(I) tti BEI(I) 6i ZS(I) sums involving Zj^ DTK dT/dK

On return,

c Z^ Z n. ZKS(I) ZKS(l) = I -^—- where y = —

i=1 y. i N

1

ZKS(2) = Z^ . I Z^.

1=1

ZKS(3) = B, - l^ B.Z^.

c 8(aa ) ZKS(^) = (aa)^=- I [aa. Z^^. . Z. -3^ -

1=1

ZKS(5) = (aa)j^

Usage: CUBDAT, FVOL, STABLIM, CUBMXT, ZSUM, SLDTK, and DZK must be called before ZKSUM.

21) DCK computes the derivative of the cubic form along the stability limit. The callinf

sequence is CALL DCK (XNT, R, T, DTK, A, B, AT, ZS, FV, ZKS, CK)

42
On entry,

XNT N

R gas constant

T temperature on stability limit DTK dT/dK

A a

B b AT da/dl

ZS(I) sums involving Z ZKS(l) sums involving Z^

On return,

CK dC/dK

Usage: CUBDAT, FVOL, STABLIM, CUBMXT, ZSUM, SLDTK, DZK, and ZKSUM must be called before DCK.

22) MSAX computes the product of a symmetric matrix and a vector. The calling sequence

is

CALL MSAX (N, A, X, B)

On entry,

N order of the matrix

A(IP) upper triangular part of A stored in packed form, i.e., A( J*( J-1 )/2+I) = A(I,J)

X(I) vector to be multiplied by A

On return,

8(1) result A X = B

23) SUMVC sums the components of a vector. The calling sequence is

CALL SUMVC (N,X,S)

On entry,

N size of vector

X(I) vector to be summed, x

On return,

N

S y X.

i = l '

24) DOT computes the dot product of two vectors. The calling sequence

CALL DOT (N,X,Y,S)

On entry,

N size of vectors

>i3
On entry,

NA

M

NAd.j)

B(I)

On return,

B(I)

X(I),Y(I) vectors to be dotted; X. X

On return,

N

S y X.Y.

i = l ' '

25) RMSOL solves the rectangular system A 2^ ^

M.' "here A is an M x n matrix. It is

assumed that M = N+1 and that a unique solution exists, i.e., the rank of A is N. The matrix

is factored by Gaussian elimination with partial pivoting. The calling sequence is

CALL RMSOL (NA, M, N, A, B)

1st dimension of the array A(I,J)

# of rows in ^ (assumed to equal N+1)

# of columns in A matrix to be factored

right hand side B

solution X

8. References

[I] Courant, R. ; John, F. Introduction to calculus and analysis, Vol. 2 New York: John Wiley

& Sons; 197^1. 187. [2] Dongarra, J.; Bunch, J. R. ; Moler, C. B. ; Stewart, G. W. LIMPACK user's guide.

Philadelphia, PA: SIAM Publications; 1978,

[3] Gibbs, J. W. On the equilibrium of heterogeneous substances (Oct. 1876-May 1877),

collected works, Vol. 1. New Haven CT: Yale Univ. Press; 1928. 55. [4] Heidemann, R. A.; Khalil, A. M. The calculation of critical points. AIChE J. 26: 769-779; 1980.

[5] Hildebrand, F. B. Methods of applied mathematics, 2nd ed. Englewood Cliffs, NJ

:

Prentice-Hall Inc.; 1965. 52. [6] Michelsen, M. L. Calculation of phase envelopes and critical points for multicomponent

mixtures. Fluid Phase Equilibria M: 1-10; 1980.

[7] Michelsen, M. L.; Heidemann, R. A. Calculations of critical points from cubic

two-constant equations of state. AIChE J. 27: 521-523; 1981.

[8] Modell, M. ; Reid, R. C. Thermodynamics and its applications. Englewood Cliffs, NJ

:

Prentice-Hall Inc.; 197^.

[9] Peng, D.-Y.; Robinson, D. B. A new two-constant equation of state. Ind. Eng. Chem.

Fundam. 15: 59-6iJ; 1976.

[10] Peng, D.-Y.; Robinson, D. B. A rigorous method for predicting the critical properties of

multicomponent systems from an equation of state. AIChE J. 23: 137-1^^1; 1977.

[II] Peters, G. ; Wilkinson, J. H. Inverse iteration, ill-conditioned equations and Newton's

method, SIAM Review 21: 339-360; 1979,
[12] Prausnitz, J. M. Molecular thermodynamics of fluid-phase equilibria. Englewood Cliffs,

NJ: Prentice-Hall Inc.; 1969.

[13] Soave, G. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem

Eng. Sci. 27: 1197-1203; 1972.

[14] Tisza, L. Generalized thermodynamics. Cambridge, MA: The M.I.T. Press; 1977.

[15] Weinhold, F. Thermodynamics and geometry. Physics Today. 23; March 1976.

[16] Hong, G.T.; Binary phase diagrams from a cubic equation of state. Ph.D Thesis. MIT;

1981 .

^5
NBS-lUA (REV. 2-8C) U.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO. NBS/TN-1313

2. Performing Organ. Report No, 3. Publ ication Date

March 1988

4. TITLE AND SUBTITLE

ON THE CALCULATION OF CRITICAL POINTS BY THE METHOD OF HEIDEMANN AND KHALIL

5. AUTHOR(S)

Brian E. Eaton

6. PERFORMING ORGANIZATION (tf joint or other thon NBS, see in struct ion sj

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

7. Contract/Grant No.

8. Type of Report & Period Covered

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

10, SUPPLEMENTARY NOTES

^ Document describes a computer program; SF-185, FlPS Software Summary, is attached.

11. ABSTRACT (A 2Q0-word or less factual summary of most significant information. If document includes a significant

bi bliography or literature survey, mention it here)

The formulation of critical point criteria by Heidemann and Khalil is analyzed

and contrasted with the original formulation of Gibbs. An extension to the

solution technique originally used by Heidemann and Khalil, and later improved

by Michelsen and Heidemann, is presented along with its detailed implementation

for a general two-constant cubic equation of state. Finally, FORTRAN software

developed for these computations is carefully documented.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolons)

algorithm; critical locus; Helmholtz free energy; mixtures; Peng-Robinson equation

of state; thermodynamics

13. AVAILABILITY

|_Xj Unlimited

I

I For Official Distribution. Do Not Release to NTIS

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

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

14. NO. OF

PRINTED PAGES

52

15. Price

USCOMM-DC 6043-P80

w US GOVERNMENT PRINTING OFFICE: 1988—575001.85.119 REGION NO. 8


NBS_ Technical 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 aaive.

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 physicisLs,

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

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—F/PS and NBSIR 's^rom 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. 1 127), 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 spxjnsors (both government and non-government). In general, initial distribution is handled hy 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