Ind. Eng. Chem. Fundam. 1980, 79, 1-10
1
REVIEW The Computation of Chemical Equilibria in Complex Systems Willlam R. Smith” Department of Mathematics and Department of Physiology and Biophysics, Dalhousie University, Halifax, Nova Scotia B3H 4H8, Canada
The problem of computing equilibrium compositions of systems with large numbers of species and elements is
reviewed. Numerical algorithms are discussed in the contexts of ideal, nonideal, and multiphase systems. Important aspects of equilibrium problems not discussed extensively elsewhere include numerical singularities incurred by some algorithms, the dependence of equilibrium compositions on the problem parameters, nonideal system computations, and the computation of special types of equilibria Considered recently in the literature. We conclude that numerical algorithms based on the classical chemical approach using stoichiometric coefficients are especially well suited for the treatment of multiphase equilibrium problems, such as those arising in metallurgical applications.
1. Introduction The computation of equilibrium compositions of complex chemical systems is an important problem in chemical engineering. By “complex” we mean systems consisting of many species (up to about 100) and elements (up to about 10). Prior to 1947, solution methods were oriented toward hand or graphical techniques, which have been reviewed by Kobe and Leland (I). The first general-purpose algorithm was devised in 1947 by Brinkley (2-7), and the development of the subject since then has been intimately linked with the development of digital computers. In spite of the problem’s importance, it is surprising that most chemistry and chemical engineering thermodynamics textbooks discuss only relatively simple systems for which calculations are readily performed by hand. That of Balzhizer et al. (8) is a notable exception. I t is the purpose of this paper to review developments in algorithms to solve this problem up to the present time. Zeleznik and Gordon (9) reviewed this field in 1968, and Van Zeggeren and Storey (IO) published a book which reviews the literature prior to 1970. Klein (11)has discussed algorithms for gaseous systems developed up to 1971, as did Holub and Vonka (12). These reviews primarily consider systems consisting of only a single (usually gaseous) ideal solution phase, although Holub and Vonka (12) have made some attempt to treat gaseous nonideality. Since the early part of this decade, interest in chemical equilibrium computations has tended to develop in mainly two directions. One is an attempt to consider nonideality properly and the other is the efficient numerical treatment of problems involving more than one phase. Recent interest in nonideality (13-15) has been primarily theoretical rather than computational, although it is clear that further developments can be expected for such problems. The late development of interest in multiphase problems (16-20)
* Department of Mathematics and Statistics, University of Guelph, Guelph, Ontario, Canada. 00 19-7874/80/ 1019-0001$01.OO/O
has arisen from metallurgical and geochemical contexts. Early interest in chemical equilibrium computations was motivated by rocketry problems in the late 1950’s and 1960’s. The calculation of the specific impulse of a rocket fuel requires the ability to calculate equilibrium compositions relevant to high temperatures and moderate pressures, where the system is usually gaseous and exhibits ideal solution behavior. Typical metallurgical and geochemical systems often consist of one multispecies phase and potentially many pure condensed phases. We will refer to such problems as “gas-solids systems” in this article. Many algorithms developed primarily for ideal gaseous systems can be very inefficient and sometimes fail completely for such problems, as we shall see in this article. In this review, we first concisely formulate the mathematical problem in several equivalent ways and discuss its general properties in section 2. In section 3 we discuss algorithms for ideal systems in the light of the formulations of section 2. Nonideal systems are treated in section 4, and in section 5 we discuss the dependence of chemical equilibrium compositions on the parameters of the problem. In section 6 we discuss numerical singularities that can arise using the algorithms of section 3, and in section 7 we consider some special topics, including some special types of equilibrium problems that have been considered recently in the literature (20,21). We summarize and make some general conclusions in section 8. A more extensive elaboration of the material of this review appears in a forthcoming article (22).
2. Mathematical Properties of the Chemical Equilibrium Problem 2.1. Problem Formulations. We assume here that we are given a temperature, T , a pressure, P, a list of N chemical species names along with their formula vectors ai, where aji is the number of atoms of element j per molecule of species i, the standard chemical potentials of the species, pi*(T,P),and the phase designation of each species; we treat chemical substances with the same formula vector but with different phase designations as dis0 1980 American Chemical Society
2
Ind. Eng. Chem. Fundam., Vol. 19, No. 1, 1980
tinct chemical species. (In this article, all vectors are column vectors and a superscript T denotes a transposed quantity.) For example, H,O(g) and H20(1)are distinct chemical species. By this convention we subsume all problems involving phase equilibria. We also assume we are given the number of moles of chemically inert species present in each phase, and the total number of gram-atoms of each of the M elements of the system, bi, which are usually positive bounded real numbers (bi corresponding to net charge may be zero, and bi may be nonpositive for some special types of problems discussed in section 7.2). We illustrate all the following problem formulations in the case of a simple example in the Appendix to this article. The first statement of the chemical equilibrium problem is Formulation P1 N
min G = Cnipi i=l
such that N
&..n. 11 1 = bI. i=l
j = 1, 2, ..., C
and
ni I0 i = 1, 2, ..., N In the above, N is the number of species, C is the rank of the matrix A = (aji) (usually C = M). P1 is thus a constrained nonlinear minimization problem, in which the Gibbs function G is minimized subject to the elementabundance and nonnegativity Constraints. The formula matrix A contains the formula vectors as columns, and 1.1 is the vector of species chemical potentials. Throughout this article pi is the usual chemical potential divided by RT (where R is the gas constant and T the absolute temperature). A second formulation results using the Kuhn-Tucker necessary conditions for the minimization formulation P1 (23). These yield the conditions Formulation P2 c pi - C U j i 7 r j = 0 (ni > 0) (la) j=1
c pi - C U j i 7 r j I0 j=1
(ni = 0 )
(Ib)
i = 1, 2, ...,N
j = 1, 2,
..., C
where a is a vector of C Lagrange multipliers. We assume in the above that pi is well defined a t ni = 0. Equations l a and ICof formulation P2 are a set of nonlinear equations in the N + C unknowns n and T . Equations lb, which usually refer to pure condensed species (e.g., solids), state that the addition of a small nonzero amount of the species will lead to an increase in the Gibbs function of the system. Two further formulations result by removing the elemental abundance constraints of eq 1 to yield an unconstrained minimization problem. We call these stoichiometric formulations of the equilibrium problem and formulations P 1 and P2 above are called nonstoichiometric formulations.
Since eq ICare a system of C linear algebraic equations in N unknowns, we may write their general solution as R
n = n* + Cultl ]=1
(2)
where n* is any particular solution of the elemental abundance equations, are real parameters, v, are linearly independent nonzero solutions of the corresponding homogeneous equations for eq IC, and hence satisfy Avl=O
(3)
and R is given by R = N - rank(A) = N - C
(4)
Equations ICand 2-4 are the essence of what is referred to as “chemical stoichiometry” (22, 24, 25). The vectors uI in eq 3 are called stoichiometric vectors and the parameters (, in eq 2 are called reaction extents. A stoichiometric equation in chemistry is essentially a short-hand way of writing eq 3 by replacing the columns of A by the corresponding species names. For our purposes, we regard eq 2 as a linear transformation from the unknown quantities n to a new set of unknowns 5. The chemical equilibrium problem may now be expressed as Formulation P3 N
min G =
Cn,(5hL,(5)
1=1
such that
n, = n,* +
R
Cu1&2 0
1=1
In formulation P 3 the Gibbs function G, considered as a function of the N-C reaction extents 4,is minimized. The element-abundance constraints have been incorporated into the stoichiometry and only the nonnegativity constraints remain. Our final formulation of the chemical equilibrium problem results from using the necessary conditions for formulation P3. The key to incorporating the nonnegativity constraints in a simple way involves forming a stoichiometric matrix u in a special form (22),which results from expressing R formula vectors (columns of A) as linear combinations of a set of C linearly independent formula vectors corresponding to species numbered from 1 to C. Then we can recast the nonnegativity constraints as constraints on the reaction extent variables. The resulting Kuhn-Tucker necessary conditions for formulation P 3 are then Formulation P4 aG
-a t -j - i=l Cujipi = 0 I0
bj+,
(nj+,> 0)
= 0)
(5b)
j = 1, 2, ..., R Equations 5b again assume (cf. eq l b ) that 1.1 is well behaved a t n’ = 0. Equations 5a are the classical form of the equi1ib:ilum conditions given in most thermodynamics texts. The second set of conditions reflects the nonnegativity constraints on the mole numbers, and they are not usually seen in textbook discussions of the problem, in spite of their practical importance. Finally, we note a t this point that all of the above formulations are independent of any assumptions as to the
Ind. Eng. Chem. Fundam., Vol. 19, No. 1, 1980
ideality or otherwise of the phases-that is, the form of the chemical potential vector p. To perform a computation in any specific case, however, we must be given p, usually in the form of some chemical potential model in algebraic form. For example, if the system consists of a single ideal gaseous phase, then eq 5 become the so-called equilibrium-constant form of the equilibrium conditions seen in most thermodynamics texts. 2.2. Existence and Uniqueness of Solutions. Existence of solutions to the chemical equilibrium problem is a relatively straightforward matter. Provided the constraint set (denoted by 3) of formulation P 1 is nonempty and p is continuous in 3, the existence of a solution is a consequence of the Weierstrass theorem in analysis (26). Uniqueness is a more delicate matter and depends on the particular form of the chemical potentials. For ideal systems, p is given by pi =
+
p,*(T,P) In xi
multispecies phase and refer to the literature for the extensions required for other situations. Later in this paper, we discuss difficulties that can arise using these algorithms for multiphase problems (section 6). The necessary and sufficient conditions are the set of nonlinear equations in (n,r) (cf. section 2.1). ”
ti
+ In (ni/ii)= j=l Cajirj
pi*
3. Algorithms for Ideal Systems 3.1. Introduction. Chemical equilibrium algorithms can be classified in many different ways (9-12). Here, we will introduce a classification scheme, based on the problem formulations of section 2.1, which emphasizes the essential similarities and differences between various algorithms in a (hopefully) more enlightening way than earlier classifications. For example the previous designation of some algorithms as “equilibrium constant methods” (9-12) has been the source of much unnecessary confusion. The use of equilibrium constants is solely a notational device and reflects no intrinsic computational properties, as we shall see below. 3.2. Nonstoichiometric Algorithms. The algorithms of Brinkley (2-7), Huff et al. (29), and White et al. (30), as well as many others (9-12) are nonstoichiometric algorithms. They are motivated by problem formulations P 1 and P2 of section 2.1. The last of the above three has been referred to as “the free energy minimization method” (31) and the first as an “equilibrium constant method” (9-12) although we shall see below that they differ in only minor ways from each other and from the second algorithm above. For brevity we consider here only a single ideal
(74
i = 1, 2, ..., N
j = 1, 2, ..., C
where the total number of moles is N
fi=Cni+z
(6)
where the standard chemical potential pi*(T,P) is a function of T and P only and xi is the mole fraction of species i in its phase. Rather complete results have been established for this case. For problems consisting of a single phase, the solution is unique (27,28).For gassolids problems, the solution is not unique only if the columns of (b,A2)are linearly dependent, where A2 consists of the columns of A correspondingto the solids (27). For systems consisting of the more than one multispecies phase, the solution need not be unique (28). In multiphase problems, the mole fractions of the species in the various phases are uniquely determined, however (28). A simple example illustrating nonuniqueness is the system consisting of (H20(1),H2(g),Oz(g))with bH:bo= 2:l at any (P,)‘2 on the vapor pressure line of HzO(l). There are an infinite number of solutions to this problem, since the relative amounts of the liquid and gaseous phases can take on any value. We will discuss equilibrium computations for these and related types of systems in section 6. Othmer (13),Ceram and Scriven ( 1 4 ) ,and Heidemann (15) considered nonideal systems and found that, when several phases are present in the problem formulation, the Gibbs function can possess local minima. This type of mathematical nonuniqueness is a well-known phenomenon in phase equilibrium calculations, and it is not surprising that it also occurs when chemical reactions are possible.
3
i=l
and z is the number of “inert” (nonreacting) species. (In terms of the A matrix for the problem, an inert species is one whose formula vector is linearly independent from those of all the other species). Such species may thus be removed from the formula matrix and considered only in eq 7c. Equations 7a are linear in the logarithms of ni and ii and eq 7b and 7c are linear in niand fi themselves. The algorithms of Brinkley (2-7), Huff et al. (29),and White et al. (30) all use the Newton-Raphson method (32)to solve eq 5 and differ only in their use of either In ni or ni as independent variables. Using ni as independent variables yields the method of White et al. (30) which was originally derived using a minimization point of view. Expanding eq 7 about a solution estimate (no,r,fio) yields
c 6nj = njo(Caij6ri+ u - p?)
(8)
1=1
j = 1, 2, ..., N C
N
N
i=l
m=l
m=l
Cari(C aimajmnmO) + bpu = C ajmnmOpmO + bj - bp (9a) j = 1, 2, .,., C
c
N
b i h i - zu =
i=l
nm0pmo m=l
(9b)
where u is a new variable defined by N
u =
C6nj/iio
j=l
(10)
Equations 9 are a set of C + 1 linear equations in (6r,u), the solution of which is substituted in eq 8 to determine the mole number changes 6n. The next estimate of n is obtained from n = no + X6n (11) where X is a relaxation parameter, chosen to speed convergence and to ensure that n satisfies the nonnegativity constraints. We will refer henceforth to the algorithm employing eq 8-11 as the RAND algorithm. The terms bj - bp in eq 9a were not present in the original formulation of the RAND algorithm (30),being due to Zeleznik and Gordon (33) and Levine (34). The treatment of inerts in the RAND algorithm is due to Apse (35). The algorithms of Brinkley (2-7) and of Huff et al. (29) use logarithmic variables and write eq 7a in the form
4
Ind. Eng. Chem. Fundam., Vol. 19, No. 1, 1980 C
For a single ideal phase, formulation P4 yields the nonlinear equations aG _ - C vjj{Fi* + In (ni(O/fi(t))l= 0
Substitution of this in eq 7b yields N
C
Caij j=l
exp( C akjrj - pj* k=l
+ In i i ) = bi
a tj
(13)
j = 1, 2,
R
exp( C a k j r k - pj*) = 1 - z / f i
j=l
k=l
The Newton-Raphson method for eq 13 yields the Brinkley algorithm C
N
i=l
m=l
E6ai(C aimajmnmO) + bjO6 In ii
= bj - bjO
(14)
j = 1, 2, ..., C C Cbi”6Ti i=l
...,R
n = no + Cvjtj
C
j=l
(15)
where
j = 1, 2, ..., C N
i=i
S t j = (aG/a
tj)o/(a2G/a
j = 1, 2,
N
-z
= A’ -
Villars (46,47)proposed solving eq 15 by iteratively solving each equation in turn. Cruise (48)proposed considering solving all equations approximately on each iteration using the Newton Raphson method
i=l
Equations 14 are a set of C + 1 linear equations in (6?r, 6 In E ) , from the solution of which the new estimates of n are obtained by eq 12. Brinkley usually expressed the elemental abundances in terms of a set of component species rather than the elements (resulting in the replacement of aimby a symbol vim and bi by q iin eq 11-14, but this is not a significant change from the above description. The incorporation of inert species in eq 14 is new. Note that the coefficient matrices of eq 9 and 14 are identical. The right-hand sides differ as a result of the use of different independent variables in each case. We should thus expect the computational properties of the two algorithms to be very similar, as is indeed the case (36). The method of Huff et al. (29)is closely related to the Brinkley algorithm. The differences are that the elements are always included as species, and it is not assumed that (no,fio)obeys eq 12 at each iteration. The Brinkley and the Huff et al. algorithms differ in an essential way only in the latter aspect. In Brinkley’s earlier papers (2-4), he used equilibrium constant notation in eq 12. However, this was a notational convenience only, due probably to the fact that he intended his work to be used by chemists, who were familiar with such notation. Unfortunately, this fact has been obscured in the literature. This is probably partly due to the mathematicians’ and chemists’ lack of familiarity with each others’ notation, and the tendency to regard equilibriumconstant notation in some way as being “old-fashioned” in comparison with the “modern free-energy minimization method”. Partly as a result of this, minor modifications of Brinkley’s algorithm have been “rediscovered”by White (36) and by Vonka and Holub (38). The above nonstoichiometric algorithms have been extended to the case of more than one phase by several workers (32,39-45). We discuss some of the difficulties that these algorithms can encounter for such problems in section 6 of this article. 3.3. Stoichiometric Algorithms. These methods, based on problem formulations P 3 and P4 of section 2, have two general advantages over the nonstoichiometric algorithms that use sets of linear equations on each iteration, such as those described in section 3.2. These advantages are realized in the case of real systems (see section 4) and for multiphase problems and for others in which singularities can occur in the linear equations of the nonstoichiometric algorithms (see section 6).
tj90
(17)
..., R
then adjusting the mole numbers by eq 16. He also incorporated the choice of an optimal stoichiometric matrix following Browne et al. (49). Hutchinson (50),Stone (51), and Bos and Meershoek (52) suggested using the full Newton-Raphson method on eq 15, yielding
Sg = (aG/a 5)o-1(a2G/a t2)o (18) This entails the solution of a set of R linear equations in as many unknowns on each iteration. Since R is usually
+
1, the number of linear very large compared to C equations occurring in the nonstoichiometric methods discussed in section 3.2, this approach is not especially worthwhile. Naphtali (53-55) suggested the use of a first-order gradient method (23)on G(0. Since first-order algorithms have rather slow convergence properties, this approach has not been used in practice to any great extent. Smith (22,56)and Smith and Missen (57) formulated the Cruise algorithm as a minimization method (formulation P3) and incorporated a relaxation parameter. This yielded significant computational improvements, and the algorithm has been used successfullyfor a number of years. 4. Algorithms for Real Systems We assume that the chemical potential is written in the general form (19) pi = pi* (T,P) + In y i x i where xi is the mole fraction of species i and yi(x) is the activity coefficient of species i, depending on the composition of the phase. We remark here that when the nonideality follows the Lewis-Randall fugacity rule (58),the phase is an ideal solution and the methods of section 3 can be used. Rather few general algorithms exist in the literature for real systems. This is partly due to the fact that generally accurate activity coefficient models have yet to be developed. Most existing nonideal system algorithms have been developed in the context of specific chemical potential models. We describe three basic general approaches to solving the equilibrium problem, each of which may be used in the context of any of the four formulations discussed in section 2. The first considers the problem as a perturbation of an ideal system problem, and the second attempts to solve the problem directly. We call these indirect and direct algorithms, respectively. The third class of methods uses an intermediate approach and we call such algorithms semi-direct. We emphasize that in sections 4.1-4.3 (see also section 2.2) the question of uniqueness of equilibria is not con-
Ind. Eng. Chem. Fundam., Vol. 19, No. 1, 1980 5
sidered. The discussion essentially assumes unique solutions exist for the problems considered. 4.1. Indirect Algorithms. Brinkley (2)first suggested the basis of approach, in conjunction with his ideal system algorithm. It has been used by Boynton (59),Cowperthwaite and Zwisler (60),Fickett (61, 62), and Holub and Vonka (12). It is most suitable for systems in which nonideality is not severe, such as some gaseous systems. We indicate possible extensions for more severely nonideal systems. Equation 19 is written as (20) pi = pi*(T,P)+ In yi + In x i y i is a function of the (unknown) equilibrium composition x*,which in turn depends on T and P. The first two terms of eq 20 are combined to give a chemical potential in ideal solution form (21) p i = p i o ( x * ( ~ , ~ )+, In ~ ,x ~ i ) Equation 21 is used iteratively in an ideal system algorithm. First, ideality is assumed, all yi are set to unity, and the ideal system equilibrium composition is computed. This composition is used to compute the next approximation to y, which is used in eq 21 to compute another ideal system equilibrium composition as the next approximation to the actual system equilibrium composition. The iterations proceed in this way until two successive ideal system equilibria differ by a sufficiently small amount. Folkman and Shapiro (63) have given sufficient conditions for convergence of this approach. We can attempt to ensure convergence of the above method by the incorporation of a parameter a and modifying eq 20 to p i = pi*(T,P) a In y i In x i (22) We then pass a through a sequence of values (a,) starting from zero and terminating at unity. At each a, we may use the iterative approach discussed above, using the solution from the previous a value as the initial estimate of the solution a t a,. Another possibility, suggested by Zeleznik and Gordon ( 9 ) , is to integrate the differential equations relating the equilibrium solution of the problem to a. These equations are
+
+
(23) Equations 23 are integrated from a = 0 to a = 1. We discuss the computation of (ani/apj*)in section 5. The terms (a In y j / a n k )are calculated from the specific chemical potential model used. 4.2. Direct Algorithms. These algorithms approach formulations Pl-P4 of section 2 directly, incorporating the nonideality from the outset of the iterations. For example, a nonideal system version of the RAND algorithm of section 3.2 uses the Newton-Raphson method on eq 4 to yield the linear correction equations in (6n,6r) C
N
C
E(a pi/anp)6nl- j=1 Caji6xj= Ea.,a.O 11 - p? 1=1 j=l J
1 = 1, 2,
(24)
..., N
McGee and Heller (64),Michels and Schneidermann (65), and Boynton (59)considered this approach in the context of gaseous systems obeying the virial equation of state. Zeleznik and Gordon (66)considered a variation of this approach for plasma systems. When the Newton-Raphson method is used on the stoichiometric formulation P4 of section 2 there result the R linear equations in t R
N
N
N
j = 1, 2,
...,R
The number of linear equations, R , ( R = N - C) to be solved on each iteration for this approach is always less than for the corresponding number N + C for eq 24, in contrast to the ideal system case (section 3). Boynton (59) has considered a method related to eq 25. For either eq 24 or 25, a general computer program can be constructed which is essentially independent of any particular chemical potential model. A subprogram which calculates p and its compositional derivatives can be appended to the main program, the latter consisting of the implementation of eq 24 or 25. 4.3. Semi-Direct Algorithms. By this term we mean algorithms which use the ideal solution approximation to the compositional derivatives of p but which use the nonideal values for p themselves. These algorithms thus accurately calculate the gradient of the Gibbs free energy (the p vector) but determine its curvature (the matrix a2G/an,anJ= a p l / a n l ) only approximately. Eriksson and Rosen (42)and Eriksson (45)have used the RAND algorithm in this way. Sanderson and Chien (67) proposed a stoichiometric algorithm of this type and used it on some very simple multiphase equilibrium problems. 5. The Effects of Problem Parameter Changes Here we consider the problem of determining the effects on the equilibrium solution of changing the problem parameter vector p = (T,P,b,p*,z)T from some initial value PO. The marginal rates of change of n with respect to the parameters define the Jacobian or sensitivity matrix (an,/ap,). The determination of this matrix is called “sensitivity analysis”. Some uses of the sensitivity matrix are discussed later. Portions of the matrix have been discussed in the case of ideal systems by Hochstim and Adams (68),Neumann (69, 70),Gordon and Zeleznik (71), Storey (72),Bigelow and Shapiro (73), and Van Zeggeren and Storey (10). Pings (74-77) considered the problem in the special case of a single stoichiometric reaction. Finally, Smith (78) discussed the general problem of determining all the elements of the sensitivity matrix, which we briefly describe here. Dinkel and Lakshmanan (79) have also recently discussed this problem. Sensitivity analysis can be considered either in the nonstoichiometric or the stoichiometric formulations of section 2.1. For the nonstoichiometric formulation, differentiation of eq 1 with respect to a parameter Pk yields (we assume, for simplicity, that all n, > 0) N
C
C(ap,/anJo(anl/apd - Ca,,(a x,/apd I-1
1=1
i = 1, 2, ..., N
C a l j 6 n l= bj - C a j l n : 1=1
1=1
j = 1, 2,
..., C
These are N + C linear equations, in contrast to the C + 1 equations of the corresponding ideal system algorithm.
(a p,/a~dO (26)
N
C
=-
N
C (anl/apda,l
1=1
j = 1, 2,
= (ab,/apk)
..., C
6
Ind. Eng. Chem. Fundam., Val. 19, No. 1, 1980
+
Equations 26 are a set of N C linear equations in the derivatives (ani/dpk) and (apj/apk), and are generally applicable to nonideal, as well as ideal systems. For ideal systems, the number of equations can be reduced in general to C 4 (78), where C$ is the number of phases. As a specific example for a single ideal phase problem, consider the parameter pk to be the number of moles of inert species. Then we have
+
c
C(a ~
j=l
N
+ (bi/iio) ( a i i / a z ) = 0
j / a ~ C )~ i p j l n t
1=1
(27)
i = 1, 2, ..., C
c (ani/az) = niO(C U j i b 7 r j / d Z ) j=l
The C
+ (l/iiO)(aii/az))
(28)
+
1 equations of eq 27 are solved for ( a r / a z ) , (aiilaz) and the required elements of the sensitivity matrix
are obtained from eq 28. All quantities in eq 26-28 with superscript 0 are evaluated at the given value po of the parameter vector. For the stoichiometric formulation P4, differentiation of eq 5 with respect to a parameter p k yields the R linear equations (again assuming all ni > 0) R
N N
N
T=
7" + X(H* -
N
N
i=l
i=l
c nihi*(7",P))/(C(ani/aT)ohi*(~,P) +
niO(ahi*/aT)oJ(32) where X is a relaxation parameter, 7" is the current estimate of T, and subscript 0 denotes equilibrium a t ( P , P ) . The sensitivity coefficients (ani/a!I? appearing in eq 32 are computed by the methods discussed previously. It is clear that any non-(T,P) problem can be handled in the above manner. Provided one has available an efficient (T,P) equilibrium algorithm and can calculate sensitivity coefficients with respect to T and/or P, a general computer program can be constructed to solve such problems which requires only a short main program implementing the equations analogous to eq 32 and using the (T,P) algorithm in a subprogram. In the literature this approach has not been used, but special algorithms have been devised, along the lines of those discussed in section 3, for solving problems of this type. The approach discussed here is probably more efficient. This is especially true if only one of the variables T and P must be determined by means of a single nonlinear equation analogous to eq 32. In such cases, simple numerical algorithms exist for the solution of one equation in one unknown which are guaranteed to converge to the solution (23). Another use of the sensitivity coefficients arises when one wishes to optimize some function of the equilibrium composition with respect to a parameter set. If this objective function is f(n) the necessary conditions for the optimum are N
( a f / a p j ) + C (af/ani)(ani/apj)= 0
m = 1, 2, ,.., R
i=l
The sensitivity coefficients are then determined from R
ani/apk =
1=1
vli(a [1/aPk)
(30)
Again, eq 29 and 30 are applicable both to nonideal and to ideal systems. Pings (74-77) considered the above equations in the special case R = 1for an ideal single-phase system. Lu et al. (80,811 extended some of this work to multiple-reaction systems. We note that, for nonideal systems, eq 29 and 30 are usually preferable to eq 26-28, since the number of linear equations that must be solved is always less. The above type of analysis can also be readily carried out to yield the second-order sensitivity coefficients (a2ni/apjapk)(22). The sensitivity matrix can be used in many ways, a few of which we will discuss here. One example is the class of all equilibrium problems that do not specify T and P as the thermodynamic constraint functions. Such "non-(T,P) problems" are readily solved using the (T,P) algorithms discussed in sections 3 and 4 in conjunction with the sensitivity coefficients. For example, comider the (ideal-system) adiabatic equilibrium at fixed enthalpy H* and pressure P*, for which we require, in addition to the equilibrium conditions, satisfaction of the equation N
nih,(T,P*) - H* = 0 i=l
(31)
where T is the unknown adiabatic reaction temperature and hi* is the molar enthalpy of pure species i. We may solve this problem by treating nias a function of the unknown T and use the Newton-Raphson method on eq 31, treating it as one equation in the single unknown, T. This yields the iteration equation
(33)
where (pi)denotes the required parameter set. The sensitivity coefficients (ani/apj)thus also appear in eq 33. The sensitivity coefficients with respect to pi* can be used to approximately compute the effects of errors 6pi* in the thermodynamic data of the problem (82). If the error-bounds on pi* are 16pi*I then the errors 6ni are bounded approximately by N
ISniI 5 CI(ani/a pi*)/ I f i ~ j * I 1=1
(34)
Alternatively, the sensitivity coefficients (ani/apj*) can be used to approximately determine the standard deviations of the errors 6n? in terms of the standard deviations 8pi*2 by means of N
6ni2=
C(ani/a ~ j *6$) ~
j= 1
(35)
In eq 35 we have assumed that all the errors in p* are independent. The above applications of sensitivity analysis are meant to be merely illustrative, and it is clear that many problems involving chemical equilibria can be fruitfully approached in this way. 6. Numerical Singularities in Equilibrium Computations Here we are concerned with situations in which the algorithms discussed in section 3 fail and the associated reasons for this failure. This is a crucially important problem for the scientist who wishes to employ these algorithms to compute equilibria. The three nonstoichiometric algorithms discussed in section 3.3 can fail for some problems in which more than one phase is present, but also at times for single-phase problems as well. The difficulty is associated with the
Ind. Eng. Chem. Fundam., Vol. 19, No. 1, 1980 7
singularity of the coefficient matrix of eq 9 and 14, which must be solved on each iteration of the methods. Some of the difficulties discussed below have been mentioned briefly in the literature by various authors, including Barnhard and Hawkins (83),Oliver et al. (31),and Gordon and McBride (84). We will analyze the source of the singularities in a general way which has not appeared in the literature. We will suggest ways to circumvent such difficulties, but we will see that some types of problems can be practically unsolvable by certain algorithms. We will also see that certain stoichiometric algorithms, like some of those discussed in section 3.3, are ideally suited to solve these types of problems. We will confine our discussion to the RAND algorithm (30),eq 8-11 of section 3.2, although is is clear that similar remarks apply to the algorithms of Brinkley (2-7) and of Huff et al. (29). We consider gas-solids problems (see section l),since this type of system occurs frequently in metallurgical applications (16-20) and illustrates the difficulties involved. The coefficient matrix of eq 9 for these problems is (31, 39-45, 56)
’=h:
)
A , A A I T b, A, Q
(36)
where AI and Az are the portions of the formula matrix referring to the gaseous and solid species, respectively. A = diag(nl, n2,..., nh),where h is the number of gaseous species; bl is the elemental abundance vector for the gaseous species, nl is the vector of gaseous species mole numbers bl = Alnl (37) Q is the 4 X 4 null matrix (4 being the total number of phases) except for QI1 = -z (the amount of inerts). 8 is always singular for problems consisting of solids alone, since all entries in P with subscript 1vanish. Thus, these algorithms cannot be used to compute equilibria among solids. 8 is also singular whenever the submatrix Q is too large, which occurs when 4 2 C + 1. For problems consisting of only a gas phase it is readily shown that 8 is never singular (in principle, but see below), since for that case 8 can be expressed as the product of a nonsingular matrix and its transpose. In general, one can show (56) that 8 can be singular only if the columns of the matrix (Az,b)are linearly dependent. Chemically, this means that b for the system can be distributed among the solids alone. The way out of the difficulty is to remove one or more solids from the species list (and hence the calculation), since the equilibrium amounts of these solids are zero. However, it is usually difficult to decide a priori which solids to remove. Even when the columns of (A2,b) are linearly independent, the nonnegativity con8traints for the solids must be taken into account when using the algorithms of section 3. If a solid persists in taking on negative values during the course of the iterations, it should be removed from the calculation. For cases in which 8 becomes singular as discussed previously, one or more solids must be removed. In the following, we will consider gas-solids problems in the context of three “classes of difficulty”, and discuss the use of the RAND algorithm on each type. Type I: Pure Gaseous Systems. Although is in principle never singular, as discussed above, finite computer word lengths can sometimes result in numerical singularity. Gordon and McBride (84) discussed the
problem consisting of (02,H2,H20)when (bH,bo)= (2,l) with T and P such that the only species present in a significant amount is H20. This yields a coefficient matrix 4 p =
+ +
E l
el
2
+ +
b)
E2
e3
(38)
ci are small numbers, and depending on their magnitude relative to the computer’s rounding errors, 8 can appear to be singular. Gordon and McBride (84) suggest alleviating such difficulties by changing the elemental abundance vector to (2 + r, 1 + r ) , where r is a suitably small number. Type 11: Columns of (Az,b)Linearly Independent. 8 cannot be singular in this case, and the only point of concern is the efficient treatment of the nonnegativity constraints on the solids. Type IIA consists of problems in which it is possible to distribute the b vector entirely among the gaseous species (i.e., the set W = (nl:Alnl= b; n j > 0) is nonempty). This case has been discussed in the literature and its treatment is straightforward. One first computes equilibrium assuming all solids have zero amounts (they are not considered in 8). Then the quantities (cf. formulation P2 of section 2)
c
AGi = p i - X a j i r j J=1
(39)
are computed for each solid. The solid with the smallest negative value of AGi is inserted into the calculation procedure and equilibrium recomputed. Tests of AGi for the remaining species are made a t the new equilibrium, and a solid is re-inserted again according to the test. This procedure continues until an equilibrium composition is computed a t which AGi 2 0 for the remaining solids. Equation 39, which arises from the non-negativity constraint conditions in formulation P2 of section 2, can be seen chemically to test whether or not the inclusion of a solid will marginally lower the system’s Gibb’s free energy. Type IIB problems are those for which the set W defined above is empty. Thus, one or more solids must be present in nonvanishing amounts at equilibrium. The difficulty is that there is often no a priori way of deciding which solids to include in a calculation. It may be feasible in some cases to insert and remove solids in the course of a calculation according to the criterion of eq 39, but this is tedious at best. For Type I1 problems numerical singularities like those occurring for Type I problems may also arise, when components of b are linearly dependent on one or more columns of A2. For example (31),in the system consisting of A1203(s)and several gaseous species composed of the elements 0 and H, when ( b ~ l , b o = ) (2,3), we can have
\o
3
2 0
o/
where x is relatively large, b~ is the elemental abundance of hydrogen, and ci are small positive numbers. The second and third columns (rows) of 8 can give numerical singularity. These singularities can be avoided in the same way as suggested by Gordon and McBride (84) for Type I problems.
8
Ind. Eng. Chem. Fundam., Vol. 19, No. 1, 1980
Type 111: Columns of (A2,b)Linearly Dependent. Here, j3 may become singular in the course of a calculation. Difficulties associated with a singular j3 are connected with the fact that the solution to the equilibrium problem need not be unique (see section 2.2). Type IIIA problems are analogous to Type IIA above. One can attempt to treat them using the RAND algorithm in the same way, provided z # 0. When t = 0, if any subset of the set of all pure condensed species is included in the calculation having formula vectors linearly dependent on b, this immediately results in j3 being singular through the linear dependence of the columns of (b1,A2). A striking simple example of this is the system consisting of (H,(g),02(g),H20(1))with ( b ~ , b o=) (2,l). If H20(1)is included in the calculation, j3 is singular. In general, if one computes the gaseous equilibrium composition, as for problems of Type IIA, and the test for the presence of solids according to eq 39 gives a positive test for the presence of a solid with a formula vector which is a multiple of b, /3 will become singular when that solid is included in the calculation. Type IIIB problems are those for which the set W defined above is empty, and thus at least one solid must be included in the calculation procedure at the outset. The difficulty is in determining a priori which solids to include. One cannot simply insert all of them and hope to be able to remove them later in the course of the calculation, since insertion of all solids automatically makes j3 singular. Type IIIB problems arise frequently in metallurgical applications. One can circumvent the difficulty by evaluating equilibria at a large number of cases and comparing them, but this becomes practically infeasible when the number of possible solids is large. Finally, Type IIIC problems are those for which the columns of A2 are linearly dependent. j3 is singular if more than C - 1 solids are present in a calculation, and sometimes even when fewer are present. Boll (85) has suggested treating multiphase problems using the Brinkley algorithm by inserting and deleting phases in the course of the calculation according to whether or not C
T , = C exp(-pi* 9
+ jCaji7rj) -1 =1
(41)
is positive or negative, respectively, when denotes summation over all species in the phase. For a pure condensed species, this test is equivalent to that using eq 39. However, this procedure will not handle all the difficult cases, especially the numerical singularities of Types I and
11. Bigelow (28) has discussed a means of avoiding singularities for Type I11 problems. He inserts fictitious inert species into each solid phase and solves a sequence of equilibrium problems in which the amounts of these fictitious inerts are gradually reduced to near-zero amounts. The solutions to this sequence of problems can be shown to be unique (and j3 is always nonsingular) and approaches the correct solution in the absence of the inerts. This technique, although it may work, seems very inefficient, since it requires the solution of several intermediate equilibrium problems. Madeley and Toguri (18,19)have proposed solving such problems in a different way, by using a first-order nonstoichiometric algorithm due to Storey and Van Zeggeren (86) in the early stages of the calculation to eliminate the appropriate solids. Ma and Shipman (87) have suggested a related approach. Again, both these approaches may work, but they seem rather inefficient, since they essen-
tially require two different algorithms in a single computation. We remark in passing here that first-order algorithms such as those suggested above were not discussed in section 3 and are not generally used in practice because their speed of convergence is very slow. The most efficient way of solving any of the above types of gas-solids problems is by using the Smith and Missen (57) modification of the Villars-Cruise (46-48) algorithm discussed in section 3.3. No numerical singularities can occur, and the nonnegativity constraints on the solids are readily treated using the test quantities (cf. problem formulation P4 of section 2)
If the solid amount is zero and AGj < 0, the solid is reinserted into the calculation. When the amount of a solid species becomes zero in the course of the calculation, it is kept at this level as long as AG, 1 0. Solids are easily inserted and removed in the course of the calculation, and even problems consisting entirely of solids are readily handled in this manner. No singularities occur since no linear equations must be solved on each iteration, and the nonnegativity constraints on solids are easily handled. Provided sufficient care is used, it should even be possible to test for situations in which the entire gas phase is absent a t equilibrium. For example, one simple way of handling the nonnegativity constraints on the amounts of species in a multispecies phase in the context of a stoichiometric algorithm is to replace the condition ni 1 0 by ni 1 t, where t is a suitable small positive number. Another possibility is to use an analogue of eq 41. Stoichiometric algorithms are especially well suited to implementation of these approaches. We thus recommend the use of this stoichiometric algorithm in general for gassolids problems over the RAND algorithm and other nonstoichiometric algorithms discussed in section 3.2. At the very least, one must be aware when using these nonstoichiometric algorithms that numerical difficulties can occur in the course of a calculation. 7. Some Special Topics 7.1. Initial Equilibrium Estimates. For the stoichiometric algorithms it is necessary, and for the nonstoichiometric algorithms it is usually desirable, that the initial estimate noof the equilibrium composition satisfy
Ano = b
(43)
n/’ > 0 (multispecies phase species) Often, this can be readily achieved by using the given amounts of a set of “reactant” species and setting the remaining species amounts to some arbitrarily small positive amounts. A more systematic way is to use the linear programming approach of Smith and Missen (57) or some of the related techniques of Clasen (88). These are readily implemented in a computer program, and the automatic generation of an initial estimate no satisfying eq 43 is readily accomplished. 7.2. Restricted Equilibrium Problems. Here we discuss equilibrium problems in which compositional constraints are imposed in addition to the elemental abundance constraints. These additional constraints may be stated either explicitly, by means of the specification of certain linear combinations of the species mole numbers, or implicitly, by restricting the number of stoichiometric equations to be less than the maximum possible number, N - C. Restricted equilibrium problems have been dis-
Ind. Eng. Chem. Fundam., Vol. 19, No. 1, 1980 9
cussed in the literature by Kaskan and Schott (89),Schott (go), Bjornbom (21), and Smith (20, 82). When the additional compositional restrictions are explicit, and may be expressed in the form
Bn = d
(44)
one may simply define a new effective system formula matrix A and elemental abundance vector b by means of A* =
b* =
(t) (:)
(45)
This yields the usual form of (unrestricted) chemical equilibrium problem for which the effective elemental abundance constraints are A*n = b* (46) Any of the algorithms discussed in section 3 may then be used for the solution of this modified problem. The examples recently discussed in the literature by Bjornbom (21) and by Smith (20,82) are all of the explicit type, although Bjornbom originally formulated them in implicit form. An example is the system consisting of (C6H5CH3,H2,CHI, C6H6),where equal amounts of CHI and C6H6are specified at equilibrium. A possible modified formula matrix and elemental abundance vector are thus A*=
f 0”
i)
We explored the use of sensitivity analysis for the chemical equilibrium problem in section 5, which we believe is an extremely useful tool that will find many future uses in problems involving chemical equilibria. Finally, we discussed some special topics in section 7, including the computation of equilibria for some special types of problems that have been of recent interest in the literature. Acknowledgment The financial assistance of the National Research Council of Canada is gratefully acknowledged as are useful discussions with Professor R. W. Missen. Appendix We illustrate formulations Pl-P4 of section 1of the text for an elementary example. We consider the system at some specified T and P consisting of the species {CH4(g), H20(g),C02(g),CO(g), C(s))with element abundances of carbon, hydrogen and oxygen given by bc, bH,and bo, respectively. Thus there are N = 5 species and M = 3 elements: numbering the species in the order given, the unknown composition is the vector nT = (nl,n2,n3,n4,n5). The system formula matrix is o2 l 0 l 0, 1 2 1
A = ( 4l 0
) 0
(Al)
which is of full rank, since C = Rank(A) = 3 = M. Problem formulation P1 is thus 5
min G = C nipi
0 0 1 - 1
i=l
(47)
such that
+ n3 + n4 + n5 = bc 4nl + 2n2 = bH n2 + 2n3 + n4 = bo
nl For problems in which implicit compositional restrictions are imposed by restricting the stoichiometrically assessible states of the system, the determination of a suitable A* and b* was first discussed by Kaskan and Schott (89)and by Schott (90) and later by Smith (20,82). 8. Summary and Conclusions The chemical equilibrium problem is mathematically a form of nonlinear programming problem, and we have discussed its mathematical and numerical aspects in this review. Many algorithms can in principle be used for its numerical solution, but all can be classed as either stoichiometric and nonstoichiometric, depending on the way in which the elemental abundance constraints are utilized in the algorithm (see section 3). For problems involving more than one phase, the nonnegativity constraints on the species amounts become of crucial importance, and numerical difficulties can beset the unwary (see section 6). Stoichiometric algorithms, once regarded as being “older and limited in scope” (31), are especially well suited for multiphase problems and will probably enjoy a resurgence of popularity in the future. In this respect, it is of interest to note that in the nonlinear programming literature, problems of the same general form as the chemical equilibrium problem (nonlinear objective function and linear equality constraints) are only recently beginning to be treated using what are essentially (in the chemical equilibrium problem context) stoichiometric algorithms (91). We have indicated possible directions in which research on algorithms for nonideal systems will develop in section 4. The efficient implementation of the species nonnegativity constraints will also be important here.
and ni 2 0. For formulation P2, we have 111 = a1 47r2
+ 112 = 2a2 + 7r3 p 3 = al + 2a3
+ a3 115 = a1 (115 > 0 ) 115 > a1 (n5 = 0 ) 114
=
a1
(A3b) Equations A3c are also part of formulation P1 and are identical with eq A2 above. Equation A3b and the last of eq A3a is essentially a test for the presence or absence of C(s) at equilibrium. For problem formulations P3 and P4 we first form a set of stoichiometric equations for the system (25). Since C = 3, R = 5 - 3 = 2 (see eq 4 of the text). A convenient set of stoichiometric vectors is
corresponding respectively to the “reactions” forming one mole of CO(g) and C(s) from the species CH,(g), H20(g), and C02(g) 7 4 CH4(g) + 3/4 CO2(g) = CO(g) + 7 2 H2O(g) (A5)
10
Ind. Eng. Chem. Fundam., Vol. 19, No. 1, 1980
f/z CH4k) + '/z CO&)
= C ( S )+ H2O(g)
If the starting composition is denoted by (nl*, n2*, n3*, n4*, n5*)problem formulation P3 is 5
min G = C ni(l)pi(U i=l
(A6)
such that nl = nl* - 1/4t1 - l/zt21 0
n2 = n2* + Y2t1+ t21 0 n3 = n3* - 7 4 t 1 - Xt2 2 0
n4 = n4* + t1 1 0 n5 = n5* + t21 0 Finally formulation P4 is the set of nonlinear equations
aG 1 AG2 = - = - Lpl + p2 - - p 3 at2 2 2 1 0 (n5 = 0)
+ p5 = 0 (n5> 0)
Literature Cited (1) Kobe, K. A., Leland, T. W., Univeristy of Texas Bureau of Engineering Research, Special Publication No. 26, University of Texas, Austin, Texas, 1954. (2) Brinkley, S.R., Jr., J. Chem. fhys., 15, 107, (1947). (3) Kandiner, H. J., Brinkley, S. R., Jr., Ind. Eng. Chem., 42, 850, 1526 (1950). (4) Brinkley, S.R., Jr., Ind. Eng. Chem., 43, 2471 (1951). (5) Brinkley, S. R., Jr., "High Speed Aerodynamics and Jet PropulsionCombustion Processes", 8. Lewis, Ed., Vd. 11, p 64, Princeton University Press, 1956. (6) Brinkley, S.R., Jr., "Proceedings of the First Conference, Western States Section, The Combustion Institute, on Kinetics Equlllbria and Performance of High Temperature Systems", G. s. Bahn and E. E. Zukoski, Ed., Butterworths, London, 1960, p 74. (7) Brinkley, S. R., Jr., paper presented at the 16th Canadlan Chemical Engineering Conference, Windsor, Ont., 1966. (8) Balzhiser, R. E., Samuels, M. R., Eliassen, J. D., "Chemical Engineering Thermodynamics", Prentice-Hall, Englewood Cllffs, N.J., 1972. (9) Zeleznik, F. J., Gordon, S.,Ind. Eng. Chem., 60, 27 (1968). (IO) Van Zeggeren, F., Storey, S. H., "The Computation of Chemical Equilibria", Cambrldge University Press, 1970. (11) Klein, M., In "Physical Chemistry, An Advanced Treatise", Eyring, H., Henderson, D., and Jost, W., Ed., Academic Press, New York, N.Y., Vol. I,1971. (12) Holub, R., Vonka, P. "The Chemical Equilibrium of Gaseous Systems", Academia, Prague, 1975. (13) Othmer, H. G., Chem. Eng. Sci., 31, 993 (1976). (14) Ceram, H. S., Scriven, L. E., Chem. Eng. Sci., 31, 163 (1976). (15) Heidemann, R. A., Chem. Eng. Sci., 33, 1517 (1978). (16) Samuels, M. R., Ind. €ng. Chem. Fundam., 10, 643 (1971). (17) Helgeson, H. C., Brown, T. H., Nlgrini, A., Jones, T. A,, Geochim. Cosmochim. Acta, 34, 569 (1970). (18) Madeley, W. G., Toguri, J., Ind. Eng. Chem. Fundam., 12, 261 (1973). (19) Madeley, W. G., Toguri, J., Can. Metall. Quart., 12, 71 (1973). (20) Smith, W. R., Ind. Eng. Chem. Fundam., 15, 227 (1976). (21) Bjornbom, P., Ind. Eng. Chem. Fundam., 14, 102 (1975). (22) Smith, W. R., in "Theoretical Chemistry, Advances and Perspectives", Vol. 5, Academic Press, New York, N.Y., 1979. (23) Wa!sh, G. R., "Methods of Optimization", Wliey, London, 1975. (24) Smlth, W. R., (Can. Math. Soc.) Appi. Math Notes, 3, 33 (1977). (25) Smith, W. R., Missen, R. W., Chem. Eng. Educ., 13, 26 (1979). (26) Hadley, "Nonllnear and Dynamic Programming" Addison-Wesley, Reading, Mass., 1964. (27) Hancock, H. J., Motzkin, T. S., ref 6, p 82. (28) Bigeiow, J. H., TR70-3, Operations Research House, Stanford University, 1970.
(29) Huff, V. N., Gordon, S.,Morrell, V. E., NACA Rept. 1037 (1951). (30) White, W. B., Johnson, S.M., Dantzig, G. B., J. Chem. fhys., 28, 751 (1958). (31) Oliver, R. C., Stephanou, S. E., Baier, R. W., Chem. Eng., 69, 121 (1962). (32) Ortega, J. M., Rheinboldt, W. C., "Iterative Solution of Nonlinear Equations in Several Variables", Academic Press, London and New York, 1970. (33) Zeleznik, F. J., Gordon, S., NASA TN D-1454 (1962). (34) Levlne, H. B., J. Chem. fhys., 36, 3049 (1962). (35) Apse, J. I., B.A.Sc. Thesis, University of Toronto, 1965. (36) Zeleznik, F. J., Gordon, S., NASA TND-473 (1960). (37) White, W. B., J. Chem. fhys., 46, 4171 (1967). (36) Vonka, P., Holub, R., Collect. Czech. Chem. Commun., 36, 2446 (1971). (39) Boynton, F. P., J . Chem. fhys., 32, 1860 (1960). (40) Boll, R. H., J. Chem. fhys., 34, 1108 (1961). (41) Raju, B. N., Krishnaswaml, C. S.,Indlan J. Techno/., 9, 99 (1966). (42) Eriksson, G., Rosen, E., Chem. Scripta, 4, 193 (1973). (43) Kubert, B. R., Stephanou, S. E., ref 6, p 166. (44) Core, T. C., Saunders, S. G., McKitttick, D. S.,"Proceedings of the Second Conference, Westem States Section, The Combustion Institute, on Kinetics Equilibria and Performance of High Temperature Systems", G. S. Bahn, Ed., Gordon and Breach, New York, N.Y., 1963, p 243. (45) Eriksson, G., Acta Chem. Scand., 25, 2651 (1971). (46) Villars, D. S., J. fhys. Chem., 63, 521 (1959). (47) Vlllars, D. S., ref 6, p 141. (48) Crulse, D. R., J. Phys. Chem., 68, 3797 (1964). (49) Browne, H. N., Williams, M. M.. Cruise, D. R., NOTS-Tech. Publ. 2434 (1960). Hutchinson, P., Chem. Eng. Sci., 17, 703 (1962). Stone, G. E., J. Chem. Educ., 43, 241 (1966). Bos, M., Meershoek, H. Q. J., Anal. Chem. Acta, 61, 185 (1972). NaphtaB, L. M.,J. Chem. Phys., 31, 263 (1959). Naphtali, L. M., ref 6, p 181. Naphtali, L. M., Ind. Eng. Chem., 53, 387 (1961). Smith, W. R., M.A.Sc. Thesis, University of Toronto, 1966. Smith, W. R., Missen, R. W., Can. J . Chem. Eng., 48, 269 (1968). Lewis, G. N., Randall, M., "Thermodynamics", McGraw-Hill, New York, N.Y., 1923. Boynton, F. P., ref 44, p 187. Cowperthwaite, M., Zwisler, W. H., SRI Publication No. 2106, Stanford Research Institute, Menio Park, Calif., 1973. Fickett, W., Phys. Fluids, 6, 997 (1963). Fickett, W., Los Alamos Report LA-6250 (1976). Folkman, J., Shapiro, N. Z., SIAM J. Appl. Math., 16. 993 (1968). McGee, H. A., Heller, G., J . Am. Rocket Soc., 32, 203 (1962). Michels, H. H., Schneldermann, S.B., ref 44, p 205. Zeleznik, F. J., Gordon, S.,Can. J. fhys., 44, 877 (1966). Sanderson, R. V., Chien, H. H. Y., Ind. Eng. Chem. frocessDes. Dev., 12, 81 (1973). Hochstlm, A. R., Adams, B., "Progress in International Research in Thermodynamic and Transport Properties, 2nd Symposium on Thermophysical Properties", Princeton University Press, 1962. Neumann, K. K., in "Progress in International Research in Thermodynamic and Transport Properties, 2nd Symposium on Thermophysical Properties", Princeton University Press, 1962. Neumann, K. K., Ber. Bunsenges. fhys. Chem., 67, 373 (1963). Gordon, S.,Zeleznik, F. J., ARS J., 32, 1195 (1962). Storey. S. H., Can. J . Chem. Eng., 43, 168 (1965). Bigeiow, J. H., Shapiro, N. Z., RAND Corp. Report P-4626 (1970). Pings, C. J., Chem. Eng. Sci., 16, 181 (1961). Pings, C. J., Chem. Eng. Sci., 18, 671 (1963). Pings, C. J., Am. Inst. Chem. Eng. J., 10, 934 (1964). Plngs, C. J., Chem. Eng. Sci., 21, 693 (1966). Smith, W. R., Can. J. Chem. Eng., 47, 95 (1969). Dinkel, J. J., Lakshmanan, R., Comp. Chem. Eng., 1, 41 (1977). Lu, C. S.,Ph.D. Thesis, California Institute of Technology, 1967. Larsen, A. H., Lu, C. S., Pings, C. J., Chem. Eng. Sci., 23, 269 (1968). Smith, W. R., Ind. Eng. Chem. Fundam., 17. 69 (1978). Barnhard, P., Hawklns, A. W., ref 44, p 235. Gordon, S., McBride, B., NASA SP-273 (1971). Boil, R. H., J. Chem. fhys., 34, 1108 (1961). Storey, S.H.,Van Zeggeren, F., Can. J. Chem. Eng.. 4Q, 59 (1970). Ma, Y. H., Shipman, C. W., AIChE J., 18, 299 (1972). Clasen, R. J., Rand Corp. Memo 4345-PR, AD-609 904. Kaskan, W. E., Schott, G. L., Combust. Flame, 6, 73 (1962). Schott, J. Chem. fhys., 40, 2065 (1964). Chllds, B., Maron, M. J., Appl. Math. Comp., 2, 273 (1976).
Received f o r review M a r c h 7, 1979 Accepted August 21, 1979
