APPLIED THERMODYNAMICSSYMPOSIUM
CALCULAT ON OF COMPLEX CHEMICAL EQUILIBRIA The calculation of complex equilibria is an interplay o f thermodynamic funda m entals and n urn e rical analysis. This review begins
by
separately examining the pertinent thermodynamic
principles and the algorithms for solving systems of nonlinear equations and is followed
by a historical review of some of the
more sign iflcan t developments
umerical calculations of chemical equilibria have N concerned chemists and chemical engineers for a great many years; however, only rarely was it necessary
~
for them to consider more than a single chemical reaction. Thus, to many individuals, chemical equilibrium was characterized by a single equilibrium constant, and equilibrium compositions could be calculated almost trivially. In a few situations, it was necessary to consider simultaneous equilibria involving perhaps two or three reactions. A relatively familiar case of such a situation is the successive ionization of polybasic acids such as the dibasic carbonic acid, HzC03, or the tribasic phosphoric acid, H3P04. As people began to study chemical processes at more extreme temperature and pressure conditions, it soon became apparent that they could no longer consider a small number of simultaneous equilibria. For example, to calculate the flame temperature for the combustion of a hydrocarbon in air it might be necessary to consider as many as 20 or more chemical reactions. As the number of reactions increased, so did the mathematical difficulties. No longer could the simultaneous equilibrium constant relations be solved in closed form, even approximately. It became necessary to use either a trial and
error, or an iterative approach to obtain solutions of the system of simultaneous equations. Several different approaches have been used to obtain solutions of the system of simultaneous equations describing chemical equilibrium. Some of the calculational methods were designed for specific problems and often took advantage of some special characteristic of the particular problem to facilitate its solution. Other methods were intended to be multipurpose schemes that could, at least in principle, be applied to any chemical equilibrium problem. Prior to the advent of computers the special purpose schemes were the most widely used. However, as computers became more generally available the tendency was to abandon the special methods in favor of the multipurpose schemes. At the present time a predominant proportion of the computations of chemical equilibria are done with the multipurpose schemes. Accordingly, we will emphasize these methods in our review and will go into a considerable amount of detail. On the other hand, only a cursory examination will be given to the special purpose schemes. We shall begin by briefly reviewing some of the thermodynamic principles and deriving the equations that are used in the computation of chemical VOL. 6 0
NO. 6
JUNE 1968
27
equilibria. This will be followed by a presentation of some of the calculating techniques that are available for solving the nonlinear chemical equilibrium equations. Having established the necessary background, we will then be able to survey the literature and point out some of the advantages and disadvantages of the various equilibrium computation schemes that have been proposed. Finally we shall complete our paper by discussing the application of equilibrium computations to typical problems. 1. T H E R M O D Y N A M I C FUNDAMENTALS Let us begin our considerations of chemical equilibria by giving a brief review of the pertinent thermodynamic fundamentals. This section will serve a dual purpose; it will act as a repository for some necessary formulas and also will introduce our notation. Those readers interested in a more detailed exposition of the subject should consult one of the several textbooks listed in the references (74, 26, 29, 30, 43, 59, 64, 77, 86). The excellent textbooks by Callen (Id),Kirkwood and Oppenheim (43),and Wilson (86) are particularly appropriate. Fundamental Relations
Every thermodynamic system can be completely characterized by any one of its fundamental relations; that is, all the thermodynamic information about a system can be derived from it. The particular fundamental relation that is used is immaterial and is selected on the basis of convenience. For the independent variables, temperature, T , and pressure, P, the appropriate fundamental relation is the one expressing the Gibbs free energy, G, in terms of T , P, and the composition 1-ariables. All other thermodynamic quantities can be obtained from G and are expressed in these same variables. It should be noted that we will consider oiily pressure-volume work in our discussion. The function, G, is sometimes called a thermodynamic potential. If one wished to work instead with pressure and entropy, S, as the independent variables, then the fundamental relation would express the enthalpy, H , in terms of P, S, and the composition variables. Since the temperature and pressure are convenient experimental variables, the Gibbs function is probably most familiar and will serve as the vehicle for our discussion. ]Ye will assume that several phases can exist in our system and that each chemical species is potentially present in each phase. We can then write our fundamental relation as
G
=
G ( T , P, ndQ)
where nia (z' = 1, 2 , . . .m;a = 1, 2, . . .p) represents the number of moles of species i in the phase a. The Greek index on the mole numbers will label the phase while the Latin index will label the chemical species. Now the Gibbs free energy is an extensive property, proportional to the amount of material in the system. This implies that Equation 1 must be a homogeneous function of degree one in n,":
G ( T , P, An:) 28
=
X G ( I ; P, n,")
INDUSTRIAL A N D ENGINEERING CHEMISTRY
From Euler's theorem on homogeneous functions we have immediately that (3)
i,"
where the chemical potential
p;
is defined to be (4)
and is homogeneous of degree zero in n,'; that is, p i a is an intensive quantity. A comparison of the differential of G, for fixed T and P, calculated from Equation 1 with that calculated from Equation 3 gives the well-known Gibbs-Duhem relationship (43): nZadpia = 0
(5)
Thermodynamic Equilibrium
When the Gibbs function is used to describe the thermodynamic system, the condition for thermodynamic equilibrium has been stated in the following manner by Callen (14): ". . . the equilibrium state niinimizes the Gibbs function over the manifold of states of constant temperature and pressure . . .." Thus at equilibrium, G is an extremum and therefore the variation in G, produced by the independent variations, must vanish. Not all variations in the mole numbers nia are independent. The variations an," must, as a minimum, satisfy the requirement that the total mass of each element is constant regardless of how the element is distributed among the different chemical species in the system. This might be considered the least constrained problem. Additional types of constraints are possible. For example, an analogous but somewhat more constrained problem is obtained by requiring that the mass of each chemical species remain constant. Such a constraint is appropriate for the study of the distribution of chemical species among several phases. Pt'e will examine in some detail the minimum constraint of conservation of the element. T o formulate these constraints mathematically let Z,(i = 1, 2, . . . I ) be the symbol for the I elements that make up the thermodynamic system. These are distributed among the various species Y,"(z = 1, 2 , . . .m; a = 1, 2 , . , .p). I n this notation the chemical foriiiula of a species can be written in the form 1
i=l
The formula numbers Q , ~represent the number of atoms of the ith element in the j t h species. If the system contains b,o(i = 1, 2 , . . . I ) gram-atoms of elementi?,, then the conservation of the elements can be written as P
"
ai& a=l j=1
- bio
=
0
i = 1, 2 ,
. . .I
(7)
If ionization is considered, then the conservation of charge can be expressed in the same form as Equation 7 where charge is assumed to be the ( I f 1)st "element." The value of b,+10 is then zero and u ~ + ~is, ,the charge on
the j t h species. Equation 7 represents 1 constraints on the variations Snp. These are supplemented by the nonI) negativity constraints n p 2 0. There are only (mp independent variations because of the constraints (Equation 7 ) . However, all of the 8nJa can be considered independently variable if we introduce 1 additional variables At, the Lagrangian multipliers. When we use Lagrangian multipliers X, (i = 1, 2, , , .Z), the condition for equilibrium becomes
-
8s
=
c (d + c Aiai,) 1,"
6n301
2
+
Chemical Reactions
Since, in the above formulation, the an," could be regarded as independent variations only at the expense of introducing 1 new variables, A,, it would seem that this formulation is perhaps not the most economical. I n chemical thermodynamics it has long been the practice to use independent variations corresponding to relatively simple processes that can be written in the form of chemical reactions. Thus, for example, we might write the reactions
A;, Yp
=
r
0
=
1, 2, . . . (mp
- I)
(11)
j,"
where
Treating the variations of n; and X i as independent not only gives PUla
+c
Xiail =
z
0
(9 1
but also Equation 7 , since the use of the X t permitted us to regard all of the n p as independently variable. We could, of course, have used any other set of independent variations expressible as a linear combination of the Snp. The only effect of this change would be to replace Equation 9 by a new set of equations formed from a linear combination of the equations in Equation 9. I n any case, Equation 7 and either Equation 9 or some linear combination of the equations in Equation 9 represent a set of ( I mp) nonlinear equations that must be solved to determine the equilibrium compositions and the Lagrangian multipliers. However, Equation 9 is interesting from still another point of view. Notice that the term involving the Lagrangian multipliers X i is independent of the phase index a. Thus, it contains the conditions for phase equilibria; that is, the chemical potential of a species is the same in all phases. This is precisely the condition that leads to the well-known phase rule of thermodynamics [Kirkwood and Oppenheim (43)l. Finally let us point out that although the X i were introduced as a mathematical device, they can be given a physical interpretation. Substituting Equation 9 into Equation 3 and using Equation 7, we immediately obtain
+
*
G = -
Xdb,O
Because of the one-to-one correspondence between the n; and the YJa,the 1 relations (Equation 7 ) among the n," implies that we have I relations among the Y,". Equivalently, this means that the matrix Ajarmust be of rank (mp - 1). T h e AT, are called the stoichiometric coefficients for the rth reaction and are chosen so that the elements are conserved in each of the reactions represented by Equation 11. This implies that the A,", must satisfy the relations
atjAj"7 = 0 Conventionally, the Aja7are taken to be positive for the products and negative for the reactants. A concrete example of Equation 11, for which the condition in Equation 12 is satisfied, is the reaction l/zcOz(g)
+ H2(g) - '/zC(s)
- HzO(g)
= 0
I t should be pointed out that although the reactions of Equation 11 may be significant kinetically, they are superfluous thermodynamically since all the compositions are already determined by Equations 7 and 9. I n a sense, these chemical reactions violate the spirit of thermodynamics since they specify a mechanism, or path, by which reactants are converted to products. As is well known, classical thermodynamics deals only with state functions, and for any thermodynamic process it is necessary only to know the two end points of the process and not the detailed path connecting them. The reactions in Equation 11 should most properly be regarded only as an artifice for introducing (mp - I ) independent variations, 8t7, one for each reaction. I n order to relate the variations Snt" to the St,, we proceed in the following manner. The change in the number of moles of Y; produced by the rth reaction can be written
i
Thus - A , represents the contribution of the element 2, to the Gibbs free energy of the system.
AUTHORS Frank J . Zeleznik and Sanford Gordon are Research Scientists for the National Aeronautics and Space Administration, Lewis Research Center, Cleveland, 0 hio. The authors extend their thanks and appreciation to Miss Holly Smith and Mrs. Bonnie McBride for their assistance and cooperation during the preparation of this paper.
= A;&
(%?7
(13)
The quantity t7is called the extent of reaction for the rth reaction and is defined by Equation 13. The total change in Y," is then obtained by summing Equation 13 over all the reactions
Snp =
c A;& r
(14)
I t is apparent from Equation 14 that we can replace the variations 8nJ" in Equation 8 with independent variations SET. As was pointed out, thermodynamically the chemical reactions are quite arbitrary and hence so, too, VOL. 6 0
NO. 6 J U N E 1 9 6 8
29
are the A;, with only the requirement that its rank be (m,b - 1). Kow substituting Equation 14 into Equation 8 and rearranging we obtain
have the independent reactions of Equation 11 written in the form 1
Y; -
c,v,j
=
0
(13)
k= 1
Using Equation 12, we find that the conditions of chemical equilibrium take the form
and also we again obtain Equation 7. Because of the use of independent variations the A, have disappeared from the problem. As a matter of fact, we could have obtained this result directly from a variation of G without introducing the hi if we had regarded G as a function of t;?. Alternately, the direct elimination of X i from Equation 9 using Equation 12 also gives Equation 16. These equations represent the equilibrium constant expressions for the chemical reactions in Equation 11. I n this equilibrium constant formulation of chemical equilibrium, the compositions are obtained by the simultaneous solutions of the set of (mp)nonlinear equations represented by Equations 7 and 16. I n passing we might point out that the expression on the left-hand side of Equation 16 is the negative of de Donders affinity for the rth reaction (86). Since the two formulations of the equations governing chemical equilibrium (Equations 7, 9 and Equations 7, 16) have a common origin, their content is the same. T h e only difference at this point is that the former set involves (mp I ) equations while the latter involves only (mp) equations. This is somewhat compensated by the fact that in the second formulation one must determine a set of stoichiometric coefficients Ajar. In a small system where the number of constituents (mp)is relatively small, the determination of the Ajarusually presents no problem; however, when (mj)is large, it is important to have a more systematic procedure for determining the stoichiometric coefficients. Brinkley (8) has discussed this problem in terms of the linear independence of a set of vectors constructed from the a i j . With each species Y j awe can associate a column vector yj”
+
where the T on the column vector indicates the transpose. Not all of these \-ectors will be linearly independent and, in fact, only c 5 1 will be linearly independent. Generally, the equality holds and only in special situations does the inequality apply. JYe will henceforth assume c = 1. Designating the linearl>-independent vectors by ci (i = 1, . . ., I), we can write
T h e phase index does not appear on v k j since y% = y j p and therefore the viCj must be independent of phase. Corresponding to the relationship in Equation 18, we 30
INDUSTRIAL AND ENGINEERING C H E M I S T R Y
where C, represents the chemical formula of the species corresponding to the independent vectors c,. The Ch are called components and need not all be from the same phase. The vjk play the role of stoichiometric coefficients. TYhen the indices a, j take the values corresponding to the components, the Equations 18 and 19 are identically zero. The choice of components is not unique and, from a theoretical point of view, all choices are equally good. The chemical reactions of Equation 19 can now be used to obtain a particular set of A,*, by comparison with Equation 11. Additionally the conservation of elements can be written in terms of the v i J . Based on Equation 19, the analog to Equation 6 would be I
The conservation of elements would now be written as a conservation of components
E vX,nja = q k O
ii.
=
1, 2,
. . ,I
(21)
3P
where the qho represents the moles of the component C, that must be introduced into the system to achieve the over-all composition. The numbers qko can be expressed in terms of the bko and the formula numbers, a i j , of the components. The conservation relations in terms of components could be used in the mininiization of free energy. However, no particular advantage is gained since this merely replaces a i j with v i j in Equation 9, causes - X i to be interpreted as the contribution of the ith component to the Gibbs free energy, and replaces bio with qio in Equation 10. The conditions of chemical equilibrium in Equation 16: correspcnding to the reactions of Equation 19, can be written in the dimensionless form
where the p k are the chemical potentials of the coinponents. The formulation of the equations for determining equilibrium compositions assumed that the thermodynamic state was specified by assigning the temperature and pressure. That is, in addition to either the set of Equations 7 and 9 or Equations 7 and 16, we had the pair of trivial equations
T
=
To
P
=
Po
(23)
where To and Po are constants. However, any two thermodynamic state functions could have been used for the same purpose. For example, to calculate the flame temperature for a constant pressure combustion, Equation 23 would be replaced by
O f t e n c h e m i c a l e q u i l i b r i u m p r o b l e m s o c c u r that
do
n o t r e q u i r e the g e n e r a l i t y
of the minimum c o n s t r a i n t p r o b l e m that is discussed
H(T,
P,)n: = HO P = Po
(24)
where H i s the enthalpy and Ho is a constant equal to the enthalpy of the reactants. For an isentropic or constant entropy process, the state may be specified as
S ( T , P,np) P
=
=
So
.
P
n:
= bp
nma = cma
i
=
+ ylAz( 5 b
i= 1 rm-1
rida
- bzo)
a=l
+
P
(rima
A”, a= 1
-c
,.)I
=
1
Po
are used are largely a matter of convenience. Often chemical equilibrium problems occur that do not require the generality of the minimum constraint problem just discussed. Typical of such problems is the class of problems related to the distribution of chemical species among several phases. A simple and familiar example of such solvent extraction problems is the distribution of iodine between water and carbon tetrachloride. An example with physiological importance is the distribution of chemicals between the interior of the red blood cells and the blood plasma (26, 67). Problems of this kind are all characterized by the fact that no “chemical reactions’’ take place; that is, the total amount of each species remains constant. Additionally, we may wish to take into account the existence of actual or imaginary semipermeable membranes. Thus, in our water-carbon tetrachloride-iodine example we can imagine the water phase to be separated from the carbon tetrachloride phase by a membrane that is permeable only to iodine. This is a permissible assumption because of the small mutual solubility of water and carbon tetrachloride. For the purpose of formulating the equilibrium equations let us suppose that each phase is surrounded by a semipermeable membrane that is permeable to all species except the mth species. Then our constraining equations are
a=l
6 [G
(25)
I n general, we require any two equations involving T , P, and nza to specify the state. T h e particular ones that
.“
ber of components or independent species. Each of the first m - 1 species is now a component while the mth species furnishes p additional components. The condition for equilibrium now becomes
1 , 2, CY
..., m
= 1, 2,
-1
(26)
. . .p
where bi0 now represents the total moles of the ith species represents the introduced into the system, and “,c number of moles of the mth species introduced into the a t h phase. Additional nontransferable species could be introduced without any essential complication. The chief difference between the reaction equilibria discussed previously and phase equilibria is the larger num-
Therefore we obtain l.4:
+ A< +
=
0
i=1,2,
. . .m
a = 1, 2,
. . .p . . .p
- 1
(27) a = 1, 2, 0 and in addition we also obtain the constraining Equations 26. Equations 27, together with Equations 26, can be used to calculate the equilibrium compositions. It should be noted that here, too, the Lagrangian multipliers could be eliminated to reduce the number of equations and obtain the equations corresponding to the “equilibrium constant’’ formulation of the problem, Pma
Ama =
Thermodynamic Derivatives
Principally we have relied on the axioms of thermodynamics, and our manipulations to this point have been quite general. There is still one more topic we can discuss without introducing simplifying assumptions. I n addition to the equilibrium compositions and the thermodynamic functions, such as enthalpy and entropy, one often would like to have the first derivatives of these functions. Fortunately it turns out that all first derivatives can be expressed in terms of three of them. For example, using the so-called Bridgman tables as tabulated by Glasstone (30))all the derivatives can be expressed in terms of C p = bH/bT, b V / b T , and dV/bP. I n comparison with the determination of equilibrium compositions, the evaluation of these three derivatives is almost trivial once the equilibrium compositions have been determined. The reason for this is that the calculation of these derivatives involves the solution of a set of linear equations whereas equilibrium compositions can be obtained only by solving a set of nonlinear equations. The enthalpy and volume of the system can be written as
VOL. 6 0
NO. 6 J U N E 1 9 6 8
31
(29) where H," dH/dn,* and VZa= b T7/dn," are the partial molar enthalpies and volumes, respectively, and where H a = c H p n p and V a = V z a n p are the enthalpy
This may be regarded as a variant of the Gibbs-Duhem relation of Equation 5. Multiplying the first members of both Equations 33 and 34 by nja and summing over a , j give
1
i
and volume of the a t h phase. Equation 29 can be regarded as the equation of state for the system since it expresses the 'i olume in terms of the temperature, pressure, and composition. The differentiation of these expressions gives C,
=
Cdanp
dn a
+
H,"
$
Alternately these relations could have been obtained directly from Equation 10 by differentiation.
(30) Form of Gibbs Free Energy
1,a
where we used the fact that H," and V," are homogeneous of degree zero in nJP. I n each of these three expressions, the second term represents the effect of composition changes produced by changes in the system temperature or pressure. To evaluate the three derivatives, Equations 30, 31, and 32, one must know dn,*/bT and dn,"/dP. These can be obtained by sol\-iiig the set of linear equations obtained by differentiation from either the set of Equations 7 and 9 or Equations 7 and 16, or perhaps the set of Equations 26 and 27. Thus, from the set of Equations 7 and 9, we obtain for the temperature derivatives
Some pragmatic considerations permit us to make a considerable simplification, not only in Equations 33 and 34 that are used for calculating derivatives, but also in all thermodynamic relations. T o accomplish this simplification we will make our first, albeit very minor, assumption. We shall assume that each phase contributes additively to the Gibbs free energy or, equivalently, that vie are neglecting all interactions among the phases
G a ( T ,P,n;)
G =
(37)
a
This is a good assumption so long as the phases are not finely dispersed in one another. Henceforth all our calculations will be within the framework of the assumption (Equation 37). This assumption makes the matrix, pIaRP, the direct sun1 of submatrices corresponding to each of the phases:
(33) Some of our subsequent discussions ill be facilitated if we have some idea of the structure of the chemical potential. Therefore, we will briefly look at the functional form of the chemical potential. A much more detailed presentation is available in the text by Kirkwood and Oppenheim (43). The chemical potential for the ith species in the a t h phase can generally be written in the form
where u pi;P = p k P a -
b2G/bnk' an,*
and where
S;
E
bS/bnja
is the partial molar entropy. Similarly the derivatives with respect to pressure are obtained by solving
T h e Lagrangian multipliers X i were given a physical interpretation by Equation 10. An analogous interpretation can be given to their temperature and pressure derivatives. Recalling that pi" is homogeneous of degree zero in n j P , we have from Euler's theorem on homogeneous functions
(35) 32
INDUSTRIAL AND ENGINEERING C H E M I S T R Y
where p," represents the contribution from ideal behavior and A p i a represents the contribution from nonideality. The quantity Apta is usually called the excess chemical potential and often must be obtaincd from experimental measurements (26, 64, 77). I n some cases it can be estimated theoretically as, for example, in the case of nonideal gases obeying a virial equation of state (55) or in the case of Debye-Huckel plasmas (90). The ideal chemical potential, pia, is characterized by an extremely simple concentration dependence, the logarithm of the mole fraction
p,"(np, T , P )
= Op;(T,
P)
+ RT In (n,"/na)
(40)
I n this equation Opza is an arbitrarily selected reference value for the chemical potential of the ith species in the octh
phase; na represents the total number of moles in the a t h phase and is defined by
T h e reference values of the chemical potentials for the gaseous phase (a = 1) are generally taken to be the chemical potentials of the pure species considered as ideal gases at a temperature T and a pressure P: Optl(T,
P)
= *ptl(T)
+ RT In P
(42)
Here *ptl( T ) is the ideal gas chemical potential of the ith species at zero pressure. This is a convenient choice since *p: can be calculated from spectroscopic constants by evaluating the canonical partition function of statistical thermodynamics. Calculations of this type have been performed for many species, and the results are available in convenient tabulations (77). Computer programs (53) are also readily available to perform additional calculations as new or revised spectroscopic data become available. The reference value of the chemical potential for condensed phases is sometimes chosen to be the chemical potential of the pure species a t the same temperature and pressure; however, other choices can be made and these are discussed by Kirkwood and Oppenheim (43). T o generate the usual equilibrium constant relations from the chemical potential, we need only apply Equations 39 and 40 to the chemical reactions of Equation 19 by substituting into the equilibrium relations (Equation 22). We will simplify the notation by assuming that we are only interested in the gaseous phase and then suppress the index a. This gives 1
(.dn)
=
j =1
+ 1, . . .m
(43)
(n~/fl)’kJ k= 1
where we’ve assumed the first 1 species were chosen as components. Here K,(T, P) is the usual mole fraction equilibrium constant given by
(44) while AK,(T, P, n;) has a corresponding definition in terms of the deviation from ideality Ap6.
position variables. This aspect of the problem will not be considered and we will assume p; is known in subsequent portions of the paper. Perhaps the problem that causes the greatest practical difficulty is the solution of the nonlinear equations that describe equilibrium. I n general, these equations cannot be solved in closed form but must be solved numerically by the application of some iteration scheme. T h e solution of a system of nonlinear equations is not a problem that is unique to thermodynamics but, on the contrary, it is a problem that occurs in all areas of science and engineering. As a result, the literature of numerical analysis contains a number of techniques that can be employed to effect a solution. We shall describe a few of these methods; however, we shall not go into their derivation nor shall we give any detailed discussion of the numerical aspects of these methods. O u r primary concern is to present briefly some of the methods that have been successfully used for equilibrium computations in the past and a few others that potentially might be used for this purpose in the future. The thermodynamic notation of the preceding section is not well suited to a discussion of methods for solving systems of nonlinear equations. I n this section we shall use the more concise notation of matrix algebra. Thus, for example, the set of Equations 7 and 9 will be written as
f = f(x) = 0
(46)
where f and x are real, N-dimensional column vectors and x represents the N independent variables. I n our case, x would stand for the composition variables and possibly the Lagrangian multipliers, the temperature, and pressure. Broadly speaking, the techniques for solving nonlinear equations have been divided into two categories by Householder (39). These categories are (1) the functional iteration methods and (2) the descent methods. The functional iteration reduces the problem of solving a system of nonlinear equations to the problem of solving a n infinite sequence of linear equations. The descent methods reduce the problem to solving an infinite sequence of single, generally nonlinear, equations, We first turn our attention to the functional iterations. Functional Iterations
Let x* represent the solution of Equation 46; that is, f(x*) = 0. Then the functional iteration is characterized by a vector function, g(x), such that
g(x*) = x*
j =1
II.
+ 1, . . .m
(45)
NUMERICAL SOLUTION OF NONLINEAR EQUATIONS
I n principle, the problem of chemical equilibrium has now been solved. I n practice, a great deal remains undone. For example, we must still obtain, from experiment or theory, the dependence of the chemical potential, p;, on the temperature, pressure, and com-
(47)
I n terms of this function and some initial estimate, xo, for x* we can define a sequence of vectors xk by the formula
xk+l - g(xJ (48) If xo is sufficiently close to x*, this sequence will converge to x* (39). Although convergence is thus guaranteed, it is often difficult to obtain a n initial estimate that is sufficiently close. This difficulty is common to all functional iterations. VOL. 6 0
NO. 6
JUNE 1 9 6 8
33
Thus in thermodynamics here is a very direct connection between minimiza-
tion and the so/ution of a system of equations
T h e various functional iterations differ from one another by the choice of g(x). Perhaps the most widely used iteration of this type is the Newton or NewtonRaphson iteration. I n this iteration the function g(x) is given by the formula
approximate Jacobian can be evaluated with the aid of N auxiliary points xki(i = 1, 2 , . . . N ) in the vicinity of xk by the formula
J(x3
0
(52)
Af.[AxxI-’
where the matrices Af, and Ax, are the matrices whose columns are formed from the column vectors f and x. Here J(x) is the Jacobian matrix-that formed from the partial derivatives o f f
is, the matrix
Af, = [f(xni)
- f(~k),f(xkz) - f(xd, . ., f(XXN)
af
J = -
ax
AX, = (Xkl
Xlc2
- Xlc,
(50)
from which it is apparent that the iteration is obtained by truncating the Taylor expansion of f(x) after the first derivative. Although the iteration will converge for a sufficiently good estimate, the method has some practical difficulties. Not only must one obtain a sufficiently good estimate, which in itself might present a considerable problem, but also at each stage of the iteration the Jacobian J must be nonsingular. Additionally, the repeated evaluation of the Jacobian and its inverse will require a considerable amount of computation time if the number of equations, N , is large. I t is possible to somewhat compensate for this last difficulty by employing a modified form of Newton’s iteration (39, 49). I n this form of the iteration, the Jacobian is not re-evaluated for each iteration but is fixed at its initial value.
. ., X~GN- XX)
g(x) = x - (Ax)(Af)-’f(x)
As a practical matter this form of the iteration should probably only be used in the later stages of a calculation. Newton’s method presupposes that the functional form of f(x) is known since one must be able to calculate the Jacobian. When the functions f(x) are so complex that it becomes impractical to differentiate them or when expressions for the functions are not known, then Kewton’s iteration cannot be used. T h e latter situation can arise in thermodynamics when the chemical potential is only available in tabular form. I n such a case, it is possible to use a generalized secant method (3, 88) which is a n example of the so-called multipoint iterative methods whose convergence has been discussed by Tornheim (73). T h e multipoint iterative methods are characterized by the fact that no derivatives are required in the evaluation of g(x). The generalized secant method is, in effect, Newton’s method with a n approximate Jacobian. If xk is the current estimate for x* then an I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
(54)
where the indices on the matrices Ax and Af have been suppressed. This method has the advantage that no derivatives of f(x) need be calculated; however, this means that when the iteration converges, we have only an approximate J(x*). Since the Jacobian is needed to calculate the thermodynamic first derivatives by means of Equations 30 to 34, this means that only approximate values of these derivatives can be obtained. The last method of the functional iteration type is the method of successive substitution (37). It is based on being able to write f(x) in the form
f(x) = x
- g(x)
Thus, the function g(x) is defined to be
g(x) = x
34
I
(53)
Thus the generalized secant method is defined by the function
T h e Equations 48 and 49 can be combined into
J ( ~ d ( ~ k +-i XJ = -f(xd
- Xk,
- f(XX)l
- f(x)
(55)
Thus, in relation to Newton’s method, the Jacobian is here taken to be the unit matrix. Descent M e t h o d s
The solution of thermodynamic equilibrium problems can be regarded either as the solution of a system of nonlinear equations or, as the direct minimization of the Gibbs free energy, G. I n thermodynamics the nonlinear equations themselves arise from the constrained minimization of G. Thus in thermodynamics there is a very direct connection between minimization and the solution of a system of equations, I n fact, the solution of any system of nonlinear equations can be regarded as the location of the minimum of some function +(x). If Equation 46 arises from the minimization of then
+
a+ f(x) = -
ax
and the connection is immediately apparent. When the functions f(x) are obtained in some other way, then we can always define a suitable function by
+
+(x)
.,
E
1 2
- fTMf
(57)
where f Tin the transpose o f f and M is any real, positive definite matrix. Because M is positive definite, the function + ( x ) , as defined by Equation 57, is always positive and vanishes only when f itself vanishes. This, in fact, is the definition of a positive definite matrix. The simplest choice for M would be the unit matrix, in which case +(x) is proportional to the sum of the squares of the components of f and is a measure of the error. I n the following discussion +(x) will represent the objective function to be minimized and could be either the Gibbs function regarded as a function of the ,$? or the equivalent of Equation 57. A review of minimization techniques has been given by Spang (69). We shall confine our attention to only a few of the many possibilities. All of the methods that we consider will generate a sequence, x k , of approximations to the minimum, x*, by the formula xk+l
=
xk
$.
(58)
Auk
T h e vector ukis arbitrary except that it must not lie in the surface + ( x ) = constant at the point, x k . I t specifies a direction of descent while the parameter X determines the size of the step. The methods under consideration differ in the choice of U k ; however, they all essentially determine X by the requirement that the one-variable function
(59) be a minimum. This is equivalent to determining X by solving the equation
.
where the gradient b+/bx is the column vector whose components are the partial derivatives of with respect to the components of x. I n principle, this gives the optimum choice of X. I n practice, however, it is often advantageous to use a n approximate value for X rather than to devote a n excessive amount of time to a search for the optimum value. Thus, Householder (39) sug-
*
gests determining a n approximate X by one step of Newton's method applied to Equation 6 0 with an initial estimate of X = 0. Spang (69) suggests some other possibilities. One should also recognize the possibility that
+
an approximate value of X might be obtained empirically based on experience with related systems of equations. There are many possibilities for the choice of a direction vector u k . For example, all of the functional iterations discussed previously could be used to supply a direction. Thus, from Newton's method, Equation 49, we could take
u = -J-'f
(61)
to obtain the descent version of Newton's method. Taking u from the method of successive substitutions in Equation 55 gives
u = -f I n general, a large number of choices of this type can be written in the form
(63)
u = -Bf
where B is some matrix. This form even encompasses the so-called gradient methods where u is related to b+/bx; for either f is already the gradient, Equation 56, or else from Equation 57, b+/bx = JTMf. I n particular, the method of steepest descent is given by the choice
u =
a+ --
(64)
dX
Perhaps the simplest choice of a descent direction is to take u along one of the coordinate lines. If one selects the j t h coordinate line, it then follows from Equation 6 0 that X is determined so as to make thejth component of the gradient equal to zero by changing only the j t h coordinate. This method of operation is known as a univariate or relaxation method. The choice of coordinate lines can be made in a number of ways. For example, they might be taken in sequence or perhaps they might be chosen so as to reduce the largest component of the gradient to zero. As the final technique for solving systems of nonlinear equations, we will briefly consider a method that is basically a descent method but incorporates some of the features of the multipoint iterations. I t possesses the advantage, like the multipoint iterations, of not requiring the calculation of a Jacobian but unlike these methods it gives the correct Jacobian when the iteration converges. Because of these two properties, it should have a considerable utility in thermodynamic calculations. This type of iteration was first proposed by Davidon (27) and was subsequently refined by Fletcher and Powell (25). Somewhat related algorithms were later suggested by Barnes (2) and Broyden (72). A numerical comparison of some of these methods was made by Rosen (62). (Note added in proof: For a recent theoretical comparison of these methods see Zeleznik, F. J., (i = 1, 2 , . . . I ) and A In nor(a = 1,2, . . . p ) as independent variables and calculate the constituent mole numbers from the components and estimates of the equilibrium constants based on the current compositions; and (3) interpret the logarithmic corrections A In x ( ~ as ) (x@*l) - dk))/x@), To obtain the Huff variant we must: (1) retain the temperature as a variable; (2) select atomic species as VOL 60
NO. 6 J U N E 1 9 6 8
39
components ( v i I = ai?,42 = b k 0 ) ; (3) treat both components and constituents as independent variables during the iteration; and (4) use all corrections in the logarithmic form rather than the linear form. Items 3 and 4 of the Huff version are quite important from a practical standpoint. For example, the use of logarithmic corrections automatically ensures the fact that all variables will remain positive if initially chosen to be positive and thus will satisfy the non-negativity constraints in the mole numbers, ng". The disadvantage of treating only the component moles as independent variables is that one needs rather good estimates for the components to obtain a convergent iteration. This problem has prompted some authors to devise schemes for obtaining initial estimates for the Newton-Raphson iteration (56, 80, 84). An alternate approach to this problem is to select as components only those species present in major amounts. This is not completely satisfactory because it creates the problem of changing components during the course of the iteration. This "optimum component" procedure has been used ( 7 7, 79, 76); however, the problem can be handled in a less complicated manner by treating all species as independent during the iteration. This often produces a convergent iteration even with relatively poor estimates. This latter method, while less susceptible to divergence, has been known to diverge. The divergent cases can be handled bb- using a descent Newton iteration. The basic disadvantage of the working Equations 82 or of their predecessors (Equation 79) is the asymmetric treatment of the species. From these equations it is obvious that the components are singled out for preferential treatment. I n the case of Equations 82 this is more obvious than true, since we demonstrated (97) that the equations could be put in a symmetric form by a simple transformation that eliminates those linear combination terms in Equations 82 that are componat-dependent. Introducing a variable u d by the definition A In n$
= ug
+ ( H i l / R T )A In T p d l + A ln nl i
=
1, 2,
,
,
.1
(84)
and substituting this into Equations 82 give as the symmetric equations determining u,,A In n" and A In T
c Vk:n"i.+T
k = 1, 2,
. . .I
">I
nu
-
nia
+ G"
CL
= 1, 2,
,
2
40
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
. .p
(Ho - H ) / R T
+
H;n;pj a31
"I
R2T2 ( 8 5 )
The substitution of Equation 84 into the first member of Equation 79, neglecting r, gives the corrections to all species in the form A In nJa
=
- pL/RT
+
I i=l
uplj
+
(H,"/RT) A In T
j
=
1, 2,
,
..VI; 01 =
1, 2,
+ A In n" , ,
(86)
.p
This form (Equations 85 and 86) of the Brinkley-Huff iteration is computationally superior to the original formulation because all species are given an equivalent treatment. Free Energy Minimization Methods
Prior to 1958 all equilibrium computations were carried out using the equilibrium constant formulation of the governinq equations. I n 1958 ll'hite, Johnson, and Dantzig (82) suggested that equilibrium compositions be calculated by "free energy minimization." Their procedure soon captured the fancy of some of the people making thermodynamic calculations, and it became the basis for a number of computer programs (5, 76, 48, 50, 83). The world of equilibrium computations was then divided into two camps, the free energ)- minimizers and the reactionary equilibrium constant formulators. It was not long before extra\Tagant claims and counterclaims of guaranteed convergence were heard from proponents of each method. So heated became the controversy that when a panel discussion was arranged to discuss equilibrium computations in 1959 ( I ) , it was necessary to divide the panel into a free energy panel and an equilibrium constant panel. To see what it was that pro\,oked such a controversy let us examine the method of free energy minimization. '4s we did in the discussion of Brinkley-Huff iteration we will deviate somewhat from the original presentation. This will enable us to incorporate some subsequent developments and it will also put the method in proper perspective. LVe \\ill formulate the method for the case where the thermodynamic state is specified by pressure and enthalpy since this will facilitate a comparison with our development of the Brinkley-Huff method ; however, it should be clear that any two thermodynamic variables will suffice. T o obtain the method, we need only apply the descent Newton-Raphson method to the Equations 7 and 24 and the dimensionless form of Equation 9 obtained by dividing the equation by R T . The function # to be minimized may be the Gibbs free energy G to which has been added one half the sum of the squares of Equations 7 and 41. It is no restriction to require that the estimates for the Lagrangian multipliers, rt - A X , / R T ,be taken as zero for each iteration since the multipliers occur linearly in Equation 9. Again, as before, we will use the na as independent variables and thus we will add Equation 41 to our set of equations and, in
b In n,"
addition, we will use, except for the nl, the logarithmic form of the variables in the expansion. With these conditions, we obtain as the Newton-Raphson equations:
b In T
, . . .p;
j = 1 , 2, . . . m
a = 1,2
a k j n p A In n p
=
bko -
k = 1,2,
bk
k = 1,2,,
..I
b In na - 0 a = 1, 2 , .
Cnp--n"-b In n? (H,"/RT) A In T = -p,"/RT
- 0
b In T
..p
(88)
T h e equations that determine the derivatives with respect to P are
. . .I
",I
a.i
n p A In n p
- nu A In nu
= n"
-
i
In nu -b b In P
n;
. . .p;
a = 1, 2,
axjnp
I
*
sa
The basic similarity of these equations and the corresponding Equations 79 of the Brinkley-Huff method is now clear; the only difference occurs between the first member of each set. As a matter of fact, if Equations 79 are written with atomic species as components and if we eliminate the nt from Equation 87 by using the first 1 equations of the first member of Equation 87 (a = 1, j = 1, 2 . . . I ) , we then obtain exactly Equations 79 But this should have = 1, 2, . . . I ) . because a t j = been obvious from the outset since it is quite immaterial whether the Lagrangian multipliers are eliminated prior to linearization or subsequent to linearization. I t should also be abundantly clear that neither method can have any inherent advantage; however, this is not meant to imply that "free energy minimization'' has no computational advantage. The benefits that exist are associated with the fact that no species are singled out for special handling. This, for example, makes it possible to incorporate the effect of the matrix r into the iteration while still retaining a small set of working equations analogous to Equations 85. The corresponding development in the case of Equations 79 is considerably more awkward. Although it is desirable to use r in the iteration, it is not absolutely necessary to do so in order to obtain a convergent iteration. However, once the equilibrium compositions have been calculated and one wishes to calculate the thermodynamic derivatives, it then becomes necessary to use r to obtain the correct values for these derivatives. The equations that must be solved in order to obtain the thermodynamic derivatives C,, bV/b T , and bV/bP from Equations 30 to 32 are very similar to the iteration Equations 87. This similarity enables us to give a simultaneous discussion of the procedure for taking account of nonideal effects. The equations for the derivatives with respect to T are given by the expressions
1
eau-
i=l
ant b In T
a = 1 , 2 , . . . p; j
b In n"
- b- H,"/RT In T = l . . .m
a,]
j = 1, 2,
'0
3 In n." b In P
b In n f p
- y n a - -bin^ n a p =
- PVj"/RT
=
b In nu b In P
,
. .m
k = 1, 2,
...I
a = 1, 2,
...p (89)
0
Apart from the fact that the coefficient matrix of Equations 87 has one more row and column than either Equation 88 or 89, the three sets of equations are identical. T o illustrate the reduction to a smaller set of working equations, we shall use the more concise matrix notation M v = w to symbolize any one of the three sets of Equations 87, 88, or 89. Further, we shall assume that the matrix M and the column vectors v and w have been partitioned so that the set of equations can be written in the form
The column vector VI is associated with the variables nja while v2 is associated with the remaining variables. If the square submatrix M11 = I r is assumed to be nonsingular, then we can write
+
This expression can now be used to eliminate v1 from the second member of Equation 90 to give:
N v ~= y
(92)
where N
y
M2z E
- MzlMI1-lM1z (93)
wz - Mz1M11-'wl
The matrix N and the vector y both contain the inverse of M11 which must either be calculated numerically or else must be known as a closed form analytical expression. Numerical inversion of M11 would defeat our avowed purpose of providing a smaller set of working equations analogous to Equations 85, and a closed form expression of M11-1 is not possible under all circumstances. However, when the contribution of r to Mll can be regarded as a perturbation, then M11-l may be replaced by its iterative expansion VOL. 6 0
NO. 6
JUNE 1 9 6 8
41
- r
where r0 G
I
The substitution of Equation 94 into Equation 93 gives as an alternate expression for N and y . m
m
N
=
k=O
N(k)
y
1
n
=
k=0
y(k)
(95)
Explicit expressions for the matrices Nck) and the vectors y ( Q can be readily constructed for the iteration Equations 87. These can be written in a compact form by first introducing the notation AuP = [Apj"k"]
where the matrices N(&) and the column vectors y(Ic)are defined by
nu =
h"611,I
(100)
and then defining the sequence of symmetric matrices
@)Acr = nu
The index k effectively gives the order to which the perturbation matrix r appears in the various terms. Therefore, we are now in a position to apply conventional perturbation theory to the solution of the reduced set of Equations 92. Writing
(97) and equating terms in Equation 92 with the same order of perturbation, we obtain the usual hierarchy of equations N(0)v2(O)= yCO? N(0)v2(k)= y(k) -
k
N(l') vz("-J*) j= 1
k = 1,2,, . . (98)
These expressions can be used to calculate v2 to the desired order of perturbation. T h e expression for VI can then be obtained from v2 by using Equation 91.
while the column vector y(")has the form
42
INDUSTRIAL AND ENGINEERING CHEMISTRY
... - .... With this notation we have Equation 102 for N(") The expressions for N(")and y(")(Equation 103) were considerably simplified by using the properties in Equation 81. The zero order approximation of Equations 98 and 99 is now easily seen to be identical to the modified Brinkley-Huff Equations 85 and 86 when the components of the Brinkley-Huff method are atomic species. Equivalently, the identity could be established by writing the free energy equations in terms of components. The original presentation of the free energy method differed considerably from our discussion. White et al. (82) derived their equations only for the ideal gas phase and minimized the quadratic approximation to the Gibbs free energy for an assigned temperature and pressure subject to the mass balance constraint (Equation 7). Their equations may be obtained from our zero order approximation by taking p = 1, deleting the
The c o n c l u s i o n in these studies of ideal systems is t h a t
if a s o l u t i o n exists it will be u n i q u e
1
I
last row and column of N(O) and striking out the.last element of y(O),setting bi equal to b,O, and interpreting all correction linearly rather than logarithmically. The implementation of their version of the free energy method requires that the estimates always satisfy the mass balance constraints. I n principle this is no problem; however, it does cause unnecessary difficulties in practice. These arise because round-off errors can cause compositions to violate mass balance even though the initial estimates satisfied mass balance. Both Levine (47) and we (97) indicated how to modify their equations within the framework of their derivation. The extension of the White method to multiphases and nonideal systems was done by a number of people (6, 7, 45, 67), all using the quadratic approximation to the Gibbs free energy. We developed the perturbation approach for incorporating the effects of nonideality on both the iteration and the calculation of thermodynamic derivatives (90). For the sake of completeness, we will now give the formulas for calculating the derivatives C,/R, b In V / b In T and b In V/d In P. I n all cases the matrices N@)are the same as Equation 102 with the last row and column deleted. The vectors y@), however, differ from Equation 103. For the temperature derivatives y ( n )has the form
while for the pressure derivatives it becomes
i=1,2
.,
)
'
"
I; r = 1 , 2 , . . . p
I n terms of the solutions to the hierarchy of Equations 98 and 99 the expressions, Equations 30 to 32 for the thermodynamic derivatives become m
/
I n addition to these derivatives with respect to T and P, we could also calculate derivatives with respect to the independent parameters bio. These would be significant in the application of thermodynamics to fluid dynamical problems when one assumes that the flow field is characterized by local thermodynamic equilibrium. The resulting formulas are similar to Equations 106 to 108 and therefore we will not reproduce them here. The free energy method we have just discussed was shown to be a descent Newton method with the objective function chosen to be essentially the Gibbs free energy. This choice of the objective function could easily be supplanted by any other suitable test function, and thus it seems inappropriate to label the method by a particular choice of test function. There is another method, due to Naphtali (57, 58) that perhaps seems more worthy of the name. He regards the Gibbs free energy as a function of the extent of reaction variables & and then uses the method of steepest descent to minimize G for assigned values of temperature and pressure. The components of the gradient vector are just the equilibrium constant expressions, as is apparent from Equation 15.
r
Because the t Tare essentially used as the independent variables, the estimates for the composition variables must satisfy the mass balance constraints. This causes the same problems here as in White's original version of the descent Newton calculation. This calculation suffers from the same drawbacks that any equilibrium constant method must endure; that is, one must write appropriate chemical reactions and one must make a judicious choice of components. Snow (68) has suggested two changes to improve the method. The first VOL. 6 0
NO. 6 J U N E 1 9 6 8
43
of these is to use a weighting matrix B (Equation 63), that is diagonal with the rth diagonal element being chosen as the smallest mole number in the 7th reaction. Second, he suggests that “parallel reactions” be written which, in effect, means making a more appropriate choice of components. Dobbins (22) has used what could be considered a univariate version of Naphtali’s iteration. The direction of descent is along the coordinate line ET that has the largest component of the gradient. I n addition, he uses atomic species as components and uses an approximate form for the Gibbs free energy during the calculations in order to reduce the calculating time. Storey and Van Zeggeren (70) used a slightly different approach in the calculation of equilibrium compositions of a n ideal gaseous phase a t an assigned temperature and pressure. Instead of working with the moles nt (the phase index has been suppressed because only the gas phase is considered) they essentially introduce logarithmic variables by defining variables zi by the relationships n1 = nioezc (110)
Table 1.
T h e nto plays the role of the current estimate of composition. They obtain a set of equations analogous to Equations 7 and 3 by minimizing the Gibbs free energy subject to the constraints that z is a unit vector and that mass-balance is preserved. These equations are then solved for z by one step of a Newton-Raphson iteration with z = 0 as an initial estimate. To facilitate the computations of z some terms in the Newton-Raphson equations are neglected. T h e length of the step in direction z is then obtained to minimize the Gibbs free energy. We’ve gone into considerable detail in the exposition of methods for calculating equilibria in the presence of chemical reactions. Concomitantly we’ve slighted the calculation of phase equilibria and ignored the problems associated with obtaining solutions to Equations 26 and 27. However, when thermodynamic data are available, the same kinds of calculational methods can be applied to the computation of phase equilibria as are used in the calculation of reaction equilibria. The only real difference is that in phase equilibria each species is a component whereas in reaction equilibria some species are components and others are constituents.
Thermodynamic Properties for an isentropic Process Assuming
Reactants: 1 ,000 50.00 3636 381 0.06384
10.080 5.000
20 000 2.500 2823 5500 0.05936 --1176.2
4-0.5 0
3,8451
50.000 1 .ooo 261 3 12492 0.05827 -- 1467.2 3.8451
17.088
17.162
-1.02213 1.50594 2.0931 9 1.14153 1.11682 1221.2 0.08415 0.12399
- 1 .01990
-1 -01822 1.43516 1.93619 1 .13762 1.11727 11 89.2 0.08583 0.12097
1 .00000 1 .00000
- 1 ,00080
- 1 00000 1 .OOOOQ
0.72823 1.I9136 1 .I9136 1261.3 0. I6062 8.16062
0.16046 0.116046
0.16039 0.16039
0.01 280 0.00001 0.06097 0.84682 0.00493 0.04169 0.02479
0.01107 0.00001 0.06362 0.86072 0.0041 7 0.03723 0.02318
0.00981 0.00001 0.05948 0.87129 0.00362 0.03390 0.02190
40 000
‘I 667 2728
1 .a50
16.846
7907 0.05886 -1308.4 3.8451 16.990
-1.02359 1 .53093 2.14609 1.14321 1.11686 1232.9 0.08376 0.12527
3.8451 15.665
-1 .OS240 1 .89599 2.72389 1.18918 1.12996 1476.8 0.08830 0.15908
-1.03129 1.65197 2.38266 1 ,15346 1.11 846 1294.7 0.08305 0.13304
-1.02541 1.56118 2.20834 1 ,14544 1.11706 1247.5 0,08340 0.12697
- 1 .00000
- 1 .OOOOO
- 1 .00000
3.84651
9
30.000
25.000 2 * 000 2771 671 8 0.05908 -1249.7 3.8451 16.926
e
2990 2958 0.06028 -936.8 3.8451 16.590
0.0
Hz
a
e
2663 10230 0.0585’2 -
Derivatives (Reacting Mixture)
( 3 In v / b In PIT ( 3 In V / b In TIP Cp, cal/g°K. Y , CPlCV
p l b In @Is Sound vel, a, m/sec In r/a In P)s (Y - 1)/Y
Y S , ( 3 In
(a
Derivutives (Nonreacting Mixture)
(3In “/a in P ) r (3In v / b In T ) p C p , cul/g O K . Y,CPl cv ys, (a In p / a In Q ) s Saund vel, a, m/see
(a In J/d In PIS - 1)/Y (A!
.ooooo
-”
1 .00000 0.76257 1 .19956 1 .19956 1521.6 0. I6636 0.16636
1 .O0000
1
0.74052 1.19298 1 ,19298 1337. ‘I 0.16177 0.16177
0.73294 1 .I9184 1.1918% 0.16095 0.16095
0.73037 1.19155 4 .I9155 1273.5 0.18875 0.16075
0.04020 0.00013 0.1291 5 0.66758 0.01 898 0.10523 0.03873
0.02036 0.00002 0.08921 0.79128 0.00847 0.06028 0.03037
0.01 542 0.00001 0.07649 0.82672 0.00612 0.04828 0.02696
0.01 395 0.00001 0.07235 0.83785 0.00545 0.04461 0.02578
i2a8.6
1.46617 2.00627 1.13918 1.11696 1203 .O 0.08499 0.12218
1 ,00600
0.72474 1.19114 1 .I9114 1242.3
~~~~
L4
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
1398.9
-
1 .00000 1.00000 0.72194
1.19103 1 .I9103 1327.8
_ _ = I -
.4
This concludes our review of calculation methods; however, we still wish to express our opinion on what is the best method of calculation. We feel that any method of calculation that can be made reliable is a good method when equilibrium compositions are the only concern. If, however, thermodynamic derivatives must also be calculated then the scale shifts strongly to the descent Newton method. I n particular, we niean the iteration equations in the form of Equation 87 or the reduced form of Equations 98 and 99. The reason for this is that the iteration and the calculation of derivatives are cast in the same form. With the other iteration methods, it is necessary to write one program to calculate compositions and another to calculate the derivatives since these can only be calculated from Equations 88 and 89 or their counterparts in the Brinkley-Huff calculation. The Brinkley-Huff equations for calculating derivatives are identical to Equations 88 and 89 in the ideal case; however, when nonideality must be taken into account, then Equations 88 and 89 are preferred because they lead, in a relatively simple manner, to reduced working equations. I t might be argued that
Equilibrium Composition 100.000 0.5000 2459 2321 9 0.05754 -1669.4 3.8451 17.379
-1,01336 1.33930 1.70665 1.13455 1.11960 1147.6 0.08973 0.11859
400.000 0.1250 21 50 79592 0.05639 -2031 .O 3.8451 17.735
500.000 0.1000 2099 96876 0.05625 -2084.1 3.8451 17.779
- 1 .00566
- 1 .00472
- 1 .00242
1.16449 1.23047 1.13999 1.13357 1069.0 0.10605 0.12280
1.14063 1.15785 1 14286 1 13749 056.7 0 11011 0 12501
1.07810 0.95338 1.15618 1.15339 1018.4 0.12560 0.13508
- 1 .ooooo
- 1 ,00000
1 .ooooo
0.71 266 1.19111 1.191 11 1183.7 0.16045 0.16045
1 .ooooo 0.69040 1.19375 1.19375 1097.0 0.16231 0.16231
00000 1 00000 C 6861 7 1.19459 1.19459 1082.9 0.16290 0.16290
0.00636 0.00000 0.04673 0.90258 0.00221 0.02437 0.01776
0.00181 0.00000 0.02318 0.95548 0.00053 0.00963 0.00937
0.00137 0.00000 0.01989 0.96234 0.00039 0.00789 0.0081 1
-1
1000.00 0.0500 1935 177475 0.05589 -2239.9 3.8451 17.892
- 1 .00000 1 .ooooo 0.67151 1.19819 1.19819 1037.9 0.16541 0.16541
0.00048 0.00000 0.01 110 0.97989 0.00012 0.00374 0.00467
thermodynamic derivatives are relatively unimportant and that when they are required they can easily be approximated by the nonreacting contributions to Equations 30 to 32. Thus the choice of a method, based on its ability to calculate derivatives, seemingly has no significance. Such an argument is fallacious on both counts. First, thermodynamic derivatives are significant for many applications of thermodynamics. Second, we will show shortly, by numerical example, that in many cases the contributions of the reactive part far exceed the nonreactive contributions to the thermodynamic derivatives. Existence of Solutions
Our discussion of the numerical methods that can be used to generate solutions of the thermodynamic equations has completely ignored the possibility that these solutions might not be unique. That is, such solutions could be either global minima or else they might merely be local minima of the Gibbs free energy. Thermodynamically the local minima correspond to the socalled metastable states while global minima correspond to true equilibrium states. Such a question, while important, cannot be given a general answer because the answer depends upon the shape of the Gibbs free energy surface for the system under consideration. T h e Gibbs free energy varies from system to system because it reflects the varying degree of nonideality present in each system. A definitive answer can only be given for ideal systems for which an analytical expression for the Gibbs function is known. Hancock and Motzkin (35) investigated a system composed of an ideal gaseous phase and pure condensed phases while Shapiro and Shapley (66) examined a system in which all phases were ideal. The conclusion in these studies of ideal systems is that if a solution exists it will be unique.
IV. APPLICATIONS The previous sections have reviewed the thermodynamic relationships and numerical techniques which permit the calculation of equilibrium compositions, thermodynamic properties, and thermodynamic first derivatives for mixtures. Theoretically calculated properties have been used in a wide variety of practical applications in chemistry and chemical engineering such as the analysis and design of chemical processing equipment, heat exchangers, steam power plants, engines, turbines, compressors, nozzles, and shock tubes. T h e usefulness of theoretical thermodynamic properties in predicting actual results depends largely on how closely the assumptions used in the calculations approximate physical reality. In some cases significantly different values may be obtained from various assumptions. I n this section a few examples have been selected (a) to illustrate the calculation and use of thermodynamic properties and (b) to point out differences due to several assumptions. These examples include (1) calculation of flame temperatures, (2) calculation of specific heat and other derivatives, (3) isentropic expansion, (4) VOL. 6 0
NO. 6 J U N E 1 9 6 8
45
Table II.
Thermodynamic Properties for an Isentropic Process Assuming Fixed" Reactants:
* j I v
Hz
+ 0.5
0 2
P
, Derivatives (Nonreacting Mixture) ( b In V / b In P)T (a In v / b In TIP Cp, cal/g, O K . 7,CP/CV Sound vel, a, m/sec (Y
-
1)Y
-.l
,00000 1.00000 0.76257 1.19956 1521.6 0.16636
- 1 ,00000
- 1 * 00000
- 1 * 00000
1.00000 0.71553 1.21550 1300.9 0.17729
1.00000 0.69673 1.22261 1353.4 0.18200
1.00000 0.69026 1.22517 1342.3 0.18379
- 1 .ooooo
%
a
1.00000 0.67601 1.23101 1319.1 0.18766
,
Composition flxed at equilibrium composition of first point.
use of derivatives to predict effect of change of an initial condition on some final condition of a thermodynamic process, and (5) use and effect of a convergence control factor A. Numerical data to illustrate the examples are given in Tables I to IV. All the examples selected are for one phase, ideal gas systems. Before proceeding with the examples it would be well to clarify what we mean by a n ideal gas. T h e chemical potential of an ideal gas p t is given by Equation 40. For this gas, the partial molar volume is
V,l
=
RT/P
(111)
Equation 111, together with Equation 29, gives the familiar equation of state
PV =
niRT = nRT
(112)
E
Equation 112 applies to an ideal gas in the sense that no interparticle forces are assumed to exist. There is no restriction, however, as to whether the nZ remain constant or are variable as a result of a change in system temperature and/or pressure. Thermodynamic texts usually assume ideal gases to be only those for which the ni in Equation 112 remain constant or, equivalently, those for which the internal energy, U , satisfies the equation
T h e only criterion necessary for a gaseous mixture to be considered ideal is Equation 112. A gas which, in addition to Equation 112, also obeys Equation 113 is a special case of an ideal gas. We will refer to ideal gas46
- 1 .ooooo
1.00000 0.60402 1.22736 1333.2 0.18524
I N D U S T R I A L A N D ENGINEERING CHEMISTRY
eous mixtures which obey both Equations 112 and 113 (ni constant) as nonreacting; while those mixtures to which Equation 112 applies but Equation 113 does not (ni variable), we will refer to as reacting. Calculation of Flame Temperature
Theoretical flame temperatures are useful for many engineering applications. They indicate the maximum temperature at which thermal energy released during chemical reactions will be available for conversion to other energy forms such as kinetic energy or work, and they also indicate possible problems relating to materials and cooling. Equilibrium compositions and adiabatic flame temperatures for specified reactants may be obtained from the solutions of iteration Equations 85 or the first member of Equation 98. (The second member of Equation 98 is identically zero for the case of ideal gases.) However, the results so obtained depend considerably on which species are assumed to exist in the reacted mixture. For example, for a stoichiometric mixture of Hz(g) and Oz(g) at a temperature of 298.1 5"K, the combustion temperature is calculated to be 4998°K if H20 is assumed to be the only combustion product. If, however, other species such as H, H2, HOz, 0, 0 2 , and O H are also assumed to exist and the combustion pressure is assumed to be 50 atm., the combustion temperature is calculated to be 3636°K. T h e results for this particular calculation are shown in the first column of Table I corresponding to the pressure ratio Po/P = 1. Other parts of this table and Table I1 will be discussed in the later examples. T h e conclusion to be reached from this example is that one should include as possible constituents all species considered likely to be significant for some set of
Composition
50.000 1 .ooo 1826 13606 0.06384 -1313.8 3.8451 15.665
-1
-
.ooooo 1.00000 0.66898 1.23400 1308.2 0.18963
100.000 0.5000 1597 25760 0.06384 1464.1 3.8451 15.665
.ooooo
-1
1.00000 0.64625 1.24425 1274.3 0.19630
400.000 0.1250 1204 901 07 0.06384 - 1709.2 3.8451 15.665
-1
.ooooo 1.00000 0.59814 1.26918 1203.5 0.21209
500.000 0.1000 1148 109950 0.06384 - 1742.5 3.8451 15.665
1000.00 0.0500 986 202700 0.06384 -1836.1 3,8451 15.665
.ooooo
- 1 .ooooo
1.00000 0.59025 1.27376 1191.2 0.21492
1.00000 0.56603 1.28886 1150.5 0.22412
-1
n
conditions. Unfortunately, for complex systems it is difficult, if not impossible, to predict a priori which species are significant for all possible conditions. For example, in the C, H, 0, N, F, C1 system, we might initially consider the 69 gaseous species indicated in Table 111. Depending on the various atom ratios which might be selected and the temperature and pressure defining the thermodynamic state, some of these species will sometimes be significant and other times insignificant. This is the type of situation for which the general computer programs discussed in Section I11 are particularly well suited, in that a priori decisions do not have to be made as to which species to include or exclude. For example, consider the reaction of perchlorylfluoride (C103F) with unsymmetrical dimethyl hydrazine (C2HsN2). The results for the stoichiometric reaction of CzHsN2 2 C103F are given in Table I11 for a pressure of 1 atm. and for several temperatures from 500" to 5000°K. Of the 69 species considered, 22 have a mole fraction greater than 0.000005 for at least some of the specified conditions. At the two lowest temperatures selected (500' and 1000°K.), only the five most stable species have a mole fraction greater than 0.000005.
+
*
Calculation of Specific Heat and Other Derivatives
Thermodynamic derivatives have many applications in chemical processes. For example, specific heat, which is one of the most commonly used thermodynamic derivatives, appears in applications such as heat transfer calculations, isentropic or isenthalpic relationships, or shock wave parameter calculations. As will be made apparent by the discussion and numerical examples that follow, the assumption of reacting or nonreacting mix-
tures may affect considerably the computed values of thermodynamic derivatives. In addition, as we will show in later sections, relationships involving these derivatives which are often valid for nonreacting mixtures may not be valid for reacting mixtures. A knowledge of specific heat and two other first derivatives such as (bV/bT), and (bP'/bP), permit one to obtain relationships among all thermodynamic first derivatives by use of tables such as devised by Bridgman (30). We will first give some discussion relating to the evaluation of these three derivatives and then give some numerical examples to illustrate differences in values between reacting and nonreacting mixtures. Expressions for evaluating C p , bV/bT and bV/bP are given in Equations 30 to 32. Values of dn,"/bP which are needed in these equations may be obtained by solution of Equations 88 and 89. Alternatively, expressions for C, and the volume derivatives in logarithmic form, b In V / b In T and b In V / b In P,are also given in Equations 106 to 108. Since we are dealing only with ideal gases in our examples, we will suppress the index a and we will need just the zero order derivatives in Equations 106 to 108; Le., b?r,(O)/dIn T , ( b In n/d In T)co),b?r,(O)/bIn P, and ( b In n / d In P)(O). These zero order derivatives may be obtained from the solution of the first member of Equation 98 where N(O) is given by Equation 102 and y(O)by Equations 104 and 105. Prior to solving Equations 88 and 89 (or the first member of Equation 98), a n explicit equation of state is needed to obtain the Vta. For the case of gases only (a = l ) , Michels (55) has evaluated the V,l for a nonideal equation of state. For gases behaving ideally, Equation 112 applies. For the case of a nonreacting gas obeying Equations 112 and 113, Equations 30 to 32 become rn
(b In V / b In T ) , = 1
(1 14)
( b In V / b In P)T = -1 The logarithmic form of the volume derivatives was selected to indicate percentage changes between reacting and nonreacting mixtures in the numerical examples which follow. T h e data in Tables I, 111, and I V illustrate the differences in magnitude between derivatives for reacting and nonreacting mixtures. The data show that for the same temperature, pressure, and equilibrium compositions two sets of values exist for these derivatives. T h e first set of values is for reacting mixtures and is obtained from the contribution of both terms in Equations 30 to 32, whereas the second set is for nonreacting mixtures V O L 6 0 NO. 6 J U N E 1 9 6 8
47
Table 111.
Thermodynamic Equilibrium Properties at Assigned Temperatures Reactants:
CzHBNz
1 0 1I t 1 J
Derivatives (Reacting Mixture)
(a In V / b In TIP Cp, cal/g OK. Yr CPICV ys, (3In P/b In
- 1 .02530 1.34669 1.10524 1.28644 1 ,25470 1911.0 0.16952 0.22266
@)e
Sound vel, a, m/sec ( b In rlb In PI, (Y
-
I)/?
-
1 0\)(310 i)[J
I
'0
0
I
l " > ) 1,i I,
i
0.00006 0.00001 0.10760 0.0001 3 0.10615 0,00002 0.00000 0.09216 0.4021 5 0.00000 0.001 63 0.01563 0.00000 0.00395 0.00004 0.00726 0.00004 0.0033 1 0.04858 0.20666 0.00373 0.00088
48
INDUSTRIAL
A N D
ENGINEERING
-1.05416 1.77045 1.88076 1.26221 1.19737 1721.1 0.12369 0.20774
- 1 ,08726 2.29605 2.67465 1.26735 1.16565 1506.7 0.09989 0.21096
-1.09330 2.56637 2.91 166 1 .24477 1.13846 1258.8 0.08378 0.19664
.
- 1 .04423 I .a9696 1 .82590
1.17358 1.12387 1051.5 0.081 42 0,14791
- 1 .01341 1.33361 0.96688 1 .14666 1.13148 913.2 0.09719 0.12790
, , OOf'OO
-! OflOii~ 1 .nocic~ii 0 " ,IO I O 9
, > t /I
k
)/ut:
I
~
(3In v / b In P)T
,out
:r
j/:
u o/oo
+ 2CIOsF 1100
i115;o
CHEMISTRY
1
*
w/
J
I .Ad/ 3 1311;. # O*J?/(,O 0. ?>/ 6 0
0.00000 0.00000 0.11358 0.00055 0.10887 0.00004 0.00000 0.06447 0.36653 0.00000 0.00521 0,04967 0.00000 0.01 136 0,00049 0.00207 0.00003 0.00539 0.05332 0.20361 0.01142 0,00339
0.00000 0* 00000 0.12583 0.00304 0.11001 0.00009 0.00001 0.02399 0.28119 0 * 00001 0.01876 0.10489 0.00000 0.03145 0.00647 0.00043 0.00002 0.00871 0.05986 0.17623 0.03482 0.01420
0. 0. 0.14034 0.01744 0.09875 0.00015 0.00002 0.00548 0.13749 0.00001 0.05884 0.15230 0.00001 0.05476 0.05438 0.00006 0.00001 0.01066 0.07353 0.09442 0.06420 0.03715
0. 0. 0.11861 0.07275 0.06989 0.00013 0.00005 0.00106 0.02944 0.00000 0.12124 0.19030 0.00001 0.03510 0.15524 0.00000 0.00000 0.00744 0.09196 0,02260 0.04380 0,04037
0. 0. 0.04733 0.16550 0.03228 0.00005 0.00009 0.00012 0.00250 0.00000 0.18033 0.21 271 0.00000 0.00998 0.21237 0.00000 0.00000 0.00245 0.1051 9 0.00183 0.01 108 0.01619
and Pressure 1 .ooo 2000 5607 0.0341 6 -1235.8 2.2436 29.271
1 .ooo 1500 din1 0.03397 - 1435.0 2.1299 29.440
1000 2787 0.03396 -1599.8 1.9967 29.446
500 1393 0.03396 - 1747.5 1.7935 29.446
-1.00197 1.05954 0.48478 1.18611 1.18378 820.1 0.14838 0.15691
- 1.00008
- 1 .ooooo
1 .ooooo 0.31329 1 ,27457 1 ,27457 599.9 0.21542 0.21542
-1 * 00000
1.00286 0.34854 1.24188 1 .24178 725.3 0.19423 0.19477
- 1 .ooooo
- 1 * 00000
- 1 .ooooo
- 1 00000
1 .ooooo 0.35891 1 ,23329 1.23329 837.0 0.18916 0.1891 6
1.00000 0.34079 1 .24700 1 .24700 726.8 0.19808 0.19808
1 .ooooo 0.31325 1.27460 1.27460 599.9 0.21544 0.21544
1 .ooooo 0.27615 1 .32343 1.32343 432.3 0.24439 0.24439
0. 0. 0.00648 0.21 442 0.00679 0.00000 0.00008 0.00000 0.00006 0. 0.21394 0.22090 0.00000 0.001 49 0.22241 0. 0. 0.00029 0.11031 0.00003 0.00093 0.001 86
0. 0. 0.00021 0.22197 0.00036 0.00000 0.00004 0.00000 0.00000 0. 0.22174 0.22217 0. 0.00008 0.22230 0. 0. 0.00001 0.1 1108 0.00000 0 00001 0.00003
0. 0. 0.00000
0. 0. 0.00000 0.22222 0. 0. 0. 0. 0. 0. 0.22222 0.22222 0. 0.00000 0.22222 0. 0. 0. 0.1 1 1 1 1 0. 0. 0.
1
.ooo
1
.ooo
1 * 00000 0.2761 4 1 .32344 1.32344 432.3 0.24439 0.24439
s
CF4(g)
0
0.22222 0.00000 0. 0.00000 0. 0. 0. 0.22222 0.22222 0. 0.00000 0.22222 0. 0. 0. 0.1 1 1 1 1 0. 0. 0.
.
CHz(g)
Czlg)
CzHz(g)
ClOztgl "3(d
FCN(d NOW
and includes the first term only, Equation 114. T h e differences in these two sets may be large or small depending on the contribution of the second term. Thus in Table 111, at T = 50OoK., the two sets ofderivatives are equal, but at T = 3500'K. the values for C, are 2.91166 and 0.38301 cal./gram"K. for reacting and nonreacting mixtures, respectively. Table I V is given to further illustrate that these differences may be dramatically large. T h e data in Table I V are for a stoichiometric Hz-02 mixture at a pressure of 0.001 atm. and at several temperatures from 1600' to 3200°K. T h e largest differences occurred at T = 2600°K., where values for C, are 23.61233 and 0.76922 cal./gram"K. for reacting and nonreacting mixtures, respectively. Thus, in this instance, the contribution of the second term in Equation 30, due to composition changes produced by a change in temperature, is many times greater than the contribution due to the heat capacities of the equilibrium species at the given temperature. Similar large differences due to the assumption of reacting or nonreacting mixtures are also found in calculated values of heat conductivity. Svehla (72) gives a n extensive tabulation of thermodynamic and transport properties for the Hz-02 system. Large differences, although not so dramatic as the one just illustrated for specific heat, are also found in the volume derivatives in Table IV. The existence of two sets of derivatives may lead one to ask which set is "correct." Unfortunately, thermodynamics cannot settle the question inasmuch as the correct answer is determined by the kinetics of the reactions involved. T h e two sets of values represent the extreme limits of infinitely fast reactions and infinitely slow reactions. If a process is such that the important reactions involved have time to reach equilibrium or near-equilibrium, then the reacting values are probably more nearly correct; contrarily, for slow reactions, the values of the nonreacting mixture are probably preferable. For intermediate situations, neither assumption may be satisfactory, and calculations using rate constants should be made (93). Isentropic Processes
Isentropic processes represent a good approximation to many actual processes. They are used in thermodynamic cycle analyses involving expansions or compressions, in various flow processes such as flow through a nozzle, or in the calculation of parameters such as the velocity of sound. I n analyses involving isentropic processes, one often would like to predict conditions at the end of a process from a knowledge of conditions at the beginning of the process. For example, knowing the pressure and temperature or pressure and volume at the beginning of a process, one might like to know the \
VOL. 6 0
NO. 6
JUNE 1 9 6 8
49
Table IV.
Thermodynamic Equilibrium Properties at Assigned Temperatures Reactants:
Hz
4- 0.5
0 2
3800 51 883333 0.16639 1 1 980.2 9.7841 6.010
0.0010 3600 491 02623 0.16622 1 1794.9 9.7340 6.016
0.0010 3400 46248747 0.16576 'I 1578 2 9.6720 6.033
O.OOI0 3200 43177742 0.16443 11262.8 9.5761 6 * 082
0.0010 3000 39407729 0.16008 10609.0 9.3641 6.247
0.0010 2800 33469673 0.14567 8819.2 8.7431 6.865
- 1 .00080
-1.001 82
- 1 .00451
1.02832 0.97781 1.55416 1.551 34 2778.2 0.34738 0.35657
1 .07456 1 .23705 1.44113 1.43466 2592.9 0.2861 5 0.30610
-1.01242 1.21838 2.08583 1.29820 1.28227 2368.5 0.19087 0.22970
-. l a 03756
1.01 184 0.89076 1.61 221 1.61092 2910.1 0.37559 0.37973
-1.11226 3.27879 14.2351 2 1.24464 1.11901 1948.1 0.06668 0.19655
0.0010
0.0010
2600 23048788 0.11178 4793.9 7.2470 8.946
Derivatives (Reacting Mixture)
( b In V / b In p ) ~ ( b In v/b In nP CP, cal/g
OK.
Y , CPICV Y S , (a In P/a In Q)s Sound vel, a, rn/sec
( b In r/b (Y
-
In
P)s
l)/Y
1.70649 5.06184 1.21416 1.17022 2161 .6 0.10725 0.17639
Derivatives (Nonreacting Mixture)
- 1 ,00000 - 1 ,00000 - 1 .00000 1 .ooooo 1 .00000 1 .00000
( b In v / b In P)?, (3In v / b In C p , cal/g
OK.
Y , CP/CV
YS,( b In p / b In (2)s Sound vel, a, m/sec
(a In J/b In P)s (Y - l)/Y
- 1 .00000
- 1 * 00000
1 .00000
0.83284 1 .65842 1 .65842 2952.7 0.39702 0.39702
0.831 63 1 ,65889 1.65889 2872.9 0.39719 0.39719
0.83027 1 .65771 1 .65771 2787.2 0.39676 0.39676
0.82817 1.65170 1.65170 2688.2 0.39456 0.39456
0.66615 0.00037 0 * 00000 0.33305 0.00032 0.0001 1
0.66552 0.00083 0.00000 0.33266 0.00071 0.00028
0.66389 0.00204 0.00001 0.33158 0.00172 0.00077
0.65922 0.00556 0.00010 0.32817 0.00460 0.00235
.
.- 1 00000
- 1 ,20831 5.62061 23.61 233 1.32619 1 ,09756 1628.6 0.05288 0.24596
- 1 .ooooo
0.3864Y
1.00000 0.80727 1 ,55907 1 .55907 2299.4 0.35859 0.35859
1 .ooooo 0.76922 1 ,40604 1 .40604 1843.3 0.28878 0.28878
0.64436 0.01681 0.00104 0.31645 0.01342 0.00792
0.59238 0.05292 0.01234 0.27683 0.03920 0.02634
0.43546 0.1 2992 0.10713 0.18197 0.08359 0.06194
1 .O00OO 0.82309 'I .62997 I . 62997 2551.2 0 .3 8 649
Additional products which were considered b u t whose mole fractions wore less than 0.000005 for all assigned condi'rions. HOdg)
temperature or volume corresponding to an assigned pressure at the end of the process. Relationships between temperature, pressure, and volume for an isentropic process involving a nonreacting ideal gas are found in any elementary thermodynamic textbook. These relationships include the following :
T / To = (P/Po)(Y-l)'Y
(115)
and
PVY = POVOY
(1 16)
CdC,
(117)
where 7 =
I n Equations 115 and 116, y is assumed to be constant or a n average value between the initial and final states of the process. Another common expression involving 7, which is also based on an isentropic relationship, is for the velocity of sound 50
INDUSTRIAL
A N D
ENGINEERING
CHEMISTRY
a =
KT
Here n must be taken as the moles per unit mass. Equations 115, 116, and 118 usually give reasonably accurate results for nonreacting ideal gas mixtures. Unfortunately, because of their ubiquitousness in the literature, these equations are often mistakenly used for reacting ideal gas mixtures for which they do not apply. M7e will first present relationships that are similar in appearance to those of Equations 115, 116, and 118 and which do apply for reacting as well as nonreacting mixtures. We will then illustrate these relationships with some numerical examples. T o obtain these relationships, we will need to evaluate the two isentropic derivatives (b In T / b In P ) , and (b In P / b In V), in terms of C,, (b In V / b In T ) p ,and ( b In V / b In P)T:
and Pressure 0.0010
0.0010
2200 11445200 0.06340 -1140.3 4.0114 15.773
2000 9509320 0,05794 -1944.0 4.4344 17.259
-1.14661 4.58921 14.5361 6 1 .24623 1 .08689 1306.4 0,04937 0.19758
- 1 .OS673
-1.01697 1.51669 2.50421 1.1 1608 1 .09745 1028.3 0.06974 0.10401
- 1 .00000
- 1 ,00000 - 1 .ooooo
0.0010 2400 15495967 0.07060 003.6 5.6552 12.709
2.54375 6.02867 1.14674 1.08518 1121 . 8 0.05316 0.12797
0.0010 1000 0296523 0.05617 -2295.0 4.2506 17.804
-. l ,
00431 1 .14736 1.22513 1.13562 1.13075 975.0 0.10454 0.11943
1 * 00000 0.72794 1.27356 1.27356 1414.1 0.21 400 0.21 400
1 ,00000 0.70056 1.21927 1.21927 1109.1 0.17904 0.17904
0.16940 0.16940
-1 .ooooo 1.00000 0.65063 1.20406 1 .20406 1006.1 0.16940 0.16940
0.20344 0.16510 0.39388 0.07410 0.08949 0.07399
0.05993 0.11429 0.70544 0.01894 0.05164 0.04975
0.01102 0.05332 0.88892 0.00319 0.01964 0.02311
0.00153 0.01831 0.96649 0.00034 0.00529 0.00804
1 .ooooo 0.67942 1 e 20406 1.20406
1077.1
0.0010 1600 7306900 0.05565 -2409.0 4.1374 17.969
evaluate the velocity of sound according to the relationship
where p is density. I n the case of a nonreacting gas, the relationships in Equations 119 to 121 reduce to
- 1 ,00090 1.03402 0.78509 1.17746 1.17640 933.3 0.14578 0.15072
- 1 .00000 '
.
Cp
1 00000 0.63503 1.21057 1.21057 946.7 0.17394 0.17394
- Cv
Integrations of Equations 119 and 120, assuming represent some mean value over the interval of integration, and using the identity of Equation 123 give
0.00012 0.00457 0.99226 0.00002 0.00090 0.00206
T / T~ = ( p / p o ()a l n T / a w s
(128)
PV75 = POVo's
(129)
For practical applications, only the exponent at the initial point may be available for use in Equations 128 and 129. We may write Equation 124 as
and
where
(121) Equation 120 may be written
V / b In
(127)
(b In T / b In P), and ( b In P/b In V ) , are constant or
a =
y s = - y / ( b In
R
(122)
where
The calculation of (b In P/b In V ) , permits one to
~/T~RRT
(130)
Here, as in Equation 118, n must be on a per-unit mass basis. From the definitions of Equations 125 and 126, we see that Equations 115, 116, and 118 for nonreacting mixtures are just special cases of Equations 128, 129, and 130, respectively. We will use the data in Tables I and I1 to illustrate the use (and misuse) of Equations 115 and 128 and also Equations 118 and 130. Table I presents the data for stoichiometric Hz-Q~at several specified pressures and at constant entropy. These data represent an isentropic process where the data in each column are the equilibrium properties at the temperature and pressure shown. Table I1 presents similar data to that of Table I with the important difference that species composition is assumed fixed (nonreacting) for all the tabulated columns. The fixed composition is taken to be that corresponding to equilibrium composition in the first column of Table I ( T = 3636°K.) and is shown again at the bottom of Table 11. For nonreacting mixtures, the values in each column for the pair y and y s given in Table I are equal as are also the values for the pair ( b In T / b In P)s and ( y - l)/y. I n contrast to this, the corresponding pair of derivatives for reacting mixtures (defined by Equations 119 to 123) vary considerably. Let us now consider the problem of estimating the VOL. 6 0
NO. 6 J U N E 1 9 6 0
51
Beginning of Process
(")
'p
b In Po
. Equation 141 Isentropic process
In T
Equation 140
T,.P,~P
Equation 142
(CP)O/CP
H
In
POVO
V
Equation 144 lsenthalpic process
In T
(s)
rocp
("V)
b In To
pa
Equation 143
(wore
b In T
p
Equation 145
b In V
cpr
S
In V
temperature which would result from an isentropic expansion over a pressure ratio Po/P of 50, assuming that we are starting from the first point in Tables I and 11. For the nonreacting mixture, from Table I1 and using Equation 115, we obtain
T
=
3636(1/50)@.166a6 = 1896°K.
Considering the long extrapolation from 3636"K., the estimate of 1896°K. compares reasonably well with the accurately calculated value of 1826°K. given in Table 11. For the case of a reacting mixture, we will start from the same point as in the previous example but will now use the data of Table I. If Equation 115, which we have just illustrated for a nonreacting gas, were mistakenly used to estimate T for a reacting gas [using the equilibrium value of (7 - l ) / r = 0.159081. the result would be
T
=
3636(1/50)0~'6908 = 1952°K.
This estimate of 1952°K. compares very poorly with the accurately calculated value of 2613°K. given in Table I. If, on the other hand, the correct formula of Equation 128 were used, we obtain
T
=
3636(1/50)0.@88a0 = 2574°K.
This estimate of 2574°K. compares well with the accurately calculated value of 2613°K. We will now look briefly- at the values of velocity of 52
INDUSTRIAL A N D E N G I N E E R I N G CHEMISTRY
sound for stoichiometric Hz-02 in Table IV. Velocity of sound (or Mach number) appears in applications such as the determination of shock wave parameters (28) or detonation velocities (89). Values for nonreacting mixtures are obtained from Equation 118 using 7 , and for reacting mixtures from Equation 130 using ys. For example, at ?' = 2600"K., a = 1843.3 and 1628.6 m./sec. for nonreacting and reacting mixtures, respectively. If Equation 118 were mistakenly used for the reacting gas (using y = 1.32619 instead of ys = 1.09756), one would obtain the incorrect value a = 1790 m./sec. T h e previous illustrations indicate the need to understand the circumstances under which approximate relationships such as Equations 114 and 118 are useful and also when they should carefully be avoided. Effect of Change in Initial Conditions on End Point of Process
We have previously discussed (32) some thermodynamic derivatives which deal with the effect of a change in initial conditions on the end point of a process. The following discussion summarizes the pertinent parts of that presentation and gives some additional derivatives. A s discussed in the first section of this paper, all thermodynamic properties of a system of known composition can be specified uniquely in terms of any two thermodynamic functions, say CY and p, which can be regarded as the coordinates of a two-dimensional space. At any
point (a,p), not only are all the thermodynamic propeities of the system determined, but it is also possible to determine the rate of change of these properties along some curve in (a,0)space. If, for example, 3. is a third thermodynamic function, the derivative (bP/da)# expresses the rate of change of P with respect to a along a curve of constant 9. This partial derivative is the usual thermodynamic first derivative which appears in thermodynamic textbooks and which we have discussed so far in this paper. By a “process” in thermodynamics, we mean that a system originally at some point (ao,00)has moved to a new point (a, 0) where (a, 0) may differ by an infinitesimal or by a finite amount from (ao,PO). An infinitesimal process can be completely characterized by a derivative of the form (bp/ba)$. A finite process can be specified by giving a starting point (ao,PO),a path (say a curve of constant #), and one of the coordinates of the end point, a. For a given path, the only independent variables of the process are the coordinates of the initial point (010, PO) and a coordinate, say a, of the final point. T h e other coordinate of the final point, P, and all other thermodynamic functions are dependent variables in the prbcess. Let cp be any dependent thermodynamic variable associated with the end point of a finite process. T h e rates of change of cp with respect to the independent thermodynamic variables of the process are of two different types. For a change in a, the usual type of derivative (bp/ba)+is needed. However, for a change in one of the coordinates of the initial points, derivatives of the type (bcp/bao)po and (ba/b/30)uo are needed. Expressions for the latter type of derivative will now be obtained. For a process taking place along a curve of constant *o(ao,
P o ) = *(a, P ) = *(a,cp>
the standard thermodynamic first derivatives and can be immediately evaluated. T h e two exceptions can be evaluated for a specified form of a relation which expresses the end point of a process in terms of the initial point Q! = 4%P o ) (135) Two forms of Equation 135 are considered in this paper: a = kl and a = k2a0, where kl and k 2 are constants. For a = k l , Equations 133 and 134 reduce to
(137) For the particular choice a = kzao, Equations 133 and 134 give
T h e right-hand sides of Equations 136 and 138 are identical. Two common processes to which Equations 137 to 139 apply are isentropic and isenthalpic. For an isentropic process, we can make the following associations: $0 = Soand tc, = S. I n addition, we can use temperature and pressure to specify the thermodynamic state and therefore make the following additional associations: 00 = TO,a0 = PO,a = P,and ao/a = Po/P. Equations 137 to 139 then become
(131)
Whatever change in $0 results from a change in the initial point (010, PO) must be equal to the change in J/. This may be expressed in differential form as
(132) Imposing alternately the conditions of constant a0 and constant Po gives the following two desired expressions for the partial derivatives of a function a t the end point of a process with respect to the initial coordinates:
Equations 140 to 142 have been evaluated for cp = In T , H, and In Vdifferentiated with respect to In POand In TOand the results are given in Table V. For a n isenthalpic process we can make the same associations as in the previous case except that now $0 = Ho and rC, = H. For this case, Equations 137 and 139 become
I n Equations 133 and 134, all the derivatives in the right-hand side except ( b a / b P ~and ) ~ ~(ba/bao)ooare VOL. 6 0
NO. 6 J U N E 1 9 6 8
53
Equations 143 to 145 have been evaluated for cp = 111 T , S, and In V differentiated with respect to In PO and In TOand are also given in Table V. Changes in In cp due to a change in initial conditions can be estimated by using the first term in a Taylor's expansion; that is,
(146) Either correction or both may be used in Equation 146. I n Table V, numerous expressions appear which inFor volve groups of terms such as P V / T or PV/TC,. ideal gases, the equation of state is given by Equation 112-Le., P V / T = nR. If H is given in cal./gram and C, and S in cal./gram°K. as in Tables I to IV, then R = 1.98726 cal./mole"K. One numerical example will be given to illustrate use of these derivatives. We will use the data in Table I which are for an isentropic process. The problem may be stated as follows: Starting with data for an initial point of P o = 50 atm. and To = 3636°K. and a final point of P = 1 atm. and T = 2613"K., calculate the temperature at P = 1 atm. for an isentropic expansion from Po = 25 atm. and To = 3636°K. This involves the derivative ( b In T / b In PO)^^,^. Using the expression in Table V for Equation 140 with cp = In T and the data from Table I , we obtain
-; (E : ::
_
_
:)p
- (0.06384)(1.98726)(1.89599) 1.93619
=
-
-
~
- 0.12423
From Equation 146, with In 9 = In T a n d A In To = 0, bln T =
(-) In POT bln
In
T
T
= 4-0.0861 10
+ 0.086110 = 7.95435
= In2613
+
+.
AlnPo = Ta,P
- 0.12423( -0.69315)
of the trivial case for which a closed form solution exists, each computational method uses a set of iterative working equations. As a result, in order to ensure convergence, each iterative method must be concerned with controlling the step size of corrections and/or with obtaining good initial estimates. I n this section we will give an example to illustrate the iterative process which starts with arbitrary initial estimates. For this purpose we have chosen the method and program with which we are most familiar (33, 92) and which were used to obtain all of the numerical results given in the tables. This descent method uses the Kewton-Raphson iterative Equations 85 with atomic species as components (or, for the case of ideal gases, the first part of Equation 98 with N(O) and y(O) given by Equations 102 and 103. This method requires no special initial estimates for composition (an arbitrary estimate is used by the program) and no constraints on mass balance or equilibrium relationships during iteration. The problem selected is the determination of the equilibrium composition of a stoichiometric mixture of a hydrocarbon (CHZ), and elemental oxygen 0 2 at a temperature of 300°K. and a combustion pressure of P = 0.01 atm. We realize that little chemical intuition is required to permit one to assume correctly that the significant gaseous species resulting from this reaction would be COz and H 2 0 . However, the point of this example is not to select the best method of solution, but rather to show that convergence is possible using arbitrary estimates which may be very poor. The course of the iteration will be indicated by two (Equation 57 with M parameters: (1) the function taken as the unit matrix), and (2) the free energy per gram of the mixture. For convenience we will call the error parameter. I n addition, we will discuss the stepsize control parameter X (Equation 58). A brief discussion of the control factor X was given in Section I1 where we pointed out that an optimum X could be obtained by minimizing Further, in Section I11 we indicated that an optimum X could also be determined by the Gibbs free energy to which has been added one half the sum of the squares of Equations 7 and 41. However, for observing the course of iteration, G alone is adequate. We use an empirically determined X such that \ A In nl 5 2 (147)
= 2848°K.
The estimate of 2848'K. compares fairly well with a n accurately calculated value of 2820°K.
A In (nt/n) Iteration a n d Convergence
A discussion was given in Sections I1 and I11 of various computational methods for obtaining chemical equilibrium compositions and the problems of obtaining convergence in iterative calculations. With the exception 54
INDUSTRIAL A N D ENGINEERING CHEMISTRY
5
2 for In ( n i / n )
A In (ni/n) = -9.212
- in
2 -18.5
(148)
(nd/n)
for In ( n i / n ) < -18.5
(149)
For each iteration, the restrictions given by Equations 147 and 148 limit the number of moles n and the mole
:j -1603
~
g
._ Y
-2030
-10
-2400
-12
-2800
-140
5
10
15 20 Ituotion n m b u
25
Figurs 1. Rora of c o w g m e fo* (CH&
-3m
30
+ (3 x/2)0r
fraction of each currently significant species (nJn > -10-8) to an increase 58. Equation 149 limits the mole fraction of each currently insignscant species (nc/n < -10-0) from becoming larger than -104. Restrictions of the type given by Equations 147 to 149 are necessary to control the ovvcorrecting which might otherwise occur when the current estimate is far from the final solution. The error parameter, is composed of three terms which show the combined error in mass balance, number of moles, and chemical potential
+,
5 (& u,,y] RT w1RT
G,>G,>G,>G,
", Control Factor and Error Parameters Table VI. During II lion fc CHI 1.5 I at T : 1°K. ai .~ P 0.01 , 'm - -
(lines of constant Gibbs free er
+
=
=
tQ.tl"
and qvcltlon
C o d fdw, h
mu-
dl0l -
0 1 2
...
7 8 9
0.01118 0.05523 0.36412 0.54797 1 .o 1 .o 1 .o 1 .o 1 .o
10 11 12 13 14 15 16 17 18 19
1 .o 0.58217 0.00021 0.00161 0.01204 0.09254 0.97068 1 .o 1 .o 1 .o
20 21 22 23 24 25
1 .o 1 .o 1 .o 1 .o 1 .o 1 .o
3 4
5 6
16 -
1.o
Ne. blwrnin-
1"s h -
...
01.148 *Hr. 141 0,148 Hc 148
EnOr pommelw, $ (Equdlon 150)
.nWY* 3/%iMG
-1/a I
1.56 1.52 1.36 3.31 6.68 1.48 1.98 2.65 3.74
X 10. X lo8 X 10'
X 1oL X 10' X 10' X 10.
X lo-'
x lo-'
7.07 X 10-'
0,149 0,148 0.148 0.148 20,148 2 0 , 148
Olbbo h.
+1185 -976 -1825 -1511 -1526 -1652 -1680 -1751 1929 -2265 -2476 -2502 -2502 -2503 -2514 -2588 -3181
3.87 x 10-8 5.25 X 10-o 7.22 X lo-" 1 . 0 7 X lo-'' 1.58 X lo-" 9.55 x lo-'$ 1.14 X lo-"
-3014.5 -3014.9 -3015.10 -301 5.15 -3015.173 -3015.179 -3015.1 81
-3083 -3022 -3013
~
Figure 2 is a sketch of free energy with contour lines representing lines of constant Gibbs free energy. Also shown in Figure 2 are a solid line and a dashed line. The solid line is a plot of Equation 152, with mass constrained at the value blo. Starting at point A, for the first iteration, the free energy decreases continuously to the solution at B. The dashed line indicates the free energy which could result from an iteration with mass unconstrained. When we start the first iteration at point C,the free energy reaches a temporary minimum at D, increases to E, and then decreases again until finally the dashed line approaches the line AB as the mass converges to bP. The previous discussion of free energy is not intended to imply an advantage for constraining mass during 56
Figure 2. Sketch of two possible pnth to cowegence. fine AB(al,m alrnr = bP): mass constrained. Curve CDEB: mass unconstrained
+
-
2.10 x 10-1 1.69 X 1.69 X lo-' 1.69X10-* 1.69X 1.68 X lo-' 2.77 X lo-: 1.85 X lo-' 6.27 X 10' 107 3.06
x
'11
INDUSTRIAL AND ENGINEERING CHEMISTRY
iteration. On the contrary, as pointed out in Section 111, with mass constrained the calculations introduce mass balance errors that lead continuously to mass imbalance for which there is no effective removal mechanism. Therefore it is a computational advantage not to impose mass constraints. This may be accomplished by including the terms (b,O - 6,) in Equation 103. The convergence example was given to indicate that convergence in chemical equilibrium calculations is usually possible with a descent Newton-Raphson method in spite of rather poor initial estimates. The fact that the 26 iterations were necessary to obtain convergence simply indicates that the initial estimates were quite poor. Typical problems, that appear more difficult to solve than the example problem, require considerably fewer iterations. For instance, consider the results in Table I11 which are for a problem involving 69 possible species. For the fmt point shown (T = 5000°K.), for which arbitrary estimates of composition were also used, 18 iterations were required. Each of the remaining nine points started with the solution to the previous point for its initial estimate and required an average of 5.4 iterations for convergence.
V.
SUMMARY
The calculation of complex chemical equilibria is an interplay of thermodynamic fundamentals and numerical analysis. In our review we have tried to consider both aspects of the problem without placing undue emphasis on either. Thus we have contrasted the two alternate, but equivalent, formulations of the conditions of chemical equilibrium: first as a minimization of the Gibbs free energy and second as a set of equilibrium constant relations. Next we considered some of the different numerical techniques that could be used in
equilibrium computations and how some of these were implemented by various people. Finally we looked a t some applications of the calculated results. REFERENCES (1) Bahn, G. S., Zukoski, E. E., Eds., “Kinetics, Equilibria and Performance of High Temperature Systems,” Butterworth, Washington, D. C., 1960. (2) Barnes, J. G. P., ComputerJ. 8, 66 (1965). (3) Bittner, L., Z.Angew. Math. Meck. 43, 111 (1963). (4) Boll, R. H., J . Chem. Phyr. 34, 1108 (1961). (5) Bowman, C. M., Foy, R. H., Rept. No. AR-3s-61, Dow Chemical Co., Midland, Mich., August 1961. (6) Boynton, F. B., in “Kinetics, Equilibria and Performance of High Temperature Spems,” G. S . Bahn, Ed., p. 187, Gordon and Breach Science Publishers, New ork, 1963. (7) Boynton, F. P., J . Chem. Phyr., 32, 1880 (1960). (8) Brinkley, S. R., Zbid., 14, 563 (1946). (9) Brinkley, S. R., Zbid., 15, 107 (1947). (10) Brinkley, S. R., in “Kinetics E uilibria and Performance of High Temperature Systems,” G. S. Bahn and E. k. Iukoski, Eds., p. 74, Butterworth, Washington, D. C.. 1960. -(11) Browne, H . N., Williams, M. M., Cruise, D. R., NAVWEPS Rept. 7043, U. S. Naval Ordnance Test Station, June 1960. (12) Broyden, C. G., Math. Comput. 19, 577 (1965). (13) Busch C. W Laderman A J Oppenheim A. K Re t. No. SSL-TN-6 IER-64-{2, Uni;: of Califorka,‘Be;keley (Avaidble fro: D E C as AD-607380j August 1964. (14) Callen, H. B., “Thermodynamics,” Wiley, New York, 1960. (15) Chu, S. T., Jet Prop. 28, 252 (1958). (16) Core, T. C., Saunders, S . G., McKittrick, P. S . , in “Kinetics, Equilibria, and Performance of High Temperature Systems,” G. S. Bahn, Ed., p. 243, Gordon and Breach Science Publishers. New York. 1963. (17) Crampel, B., Barrere S Lemaitre P Jaubert J Rept. No. N. T. 90, Office National d’Etudes et dekdrherches Akrdipatiales,’19’i5. (18) Crisman, P. A., Goldwasser, S. R., Petrozzi, P. J., in“Proceedings ofthe Propellant Conference,” p. 293, Eng. Expt. Station, Ohio State Univ., Columbus, 1960. (19) Cruise, D. R., J. Phyr. Chem. 68, 3797 (1964). (20) Dantzig, G . B.,DeHaven, J. C., Ibid., 36, 2620 (1962). (21) Davidon, W. C., AEC Report ANL-5990 (Rev.), Argonne National Lab., Lemont, Ill., November 1959. (22) Dobbins, T. O., Rept. No. WADC TR-59-757, Wright-Patterson AFB, Ohio, December 1959. (23) Donegan, A. J., Farber, M., Jet Prop. 26, 164 (1956). (24) Erickson, W. D., Kemper, J. T., Allison, D. O., Tech. Note D-3488, NASA, Washington, D. C., August 1966. (25) Fletcher, R., Powell, M . J. D., ComputerJ. 6, 163 (1963/1964). (26) Flugge, S., Ed., “Handbuch Der Physik,” Vol. 12, “Thermodynamics of Gases,” Springer-Verlag. Berlin, 1958. (27) Freudenstein, F., Roth, B., J . Assoc. Comput. Mach. 10, 550 (1963). (28) Gaydon, A. G., Hurle, I. R., “The Shock Tube in High-Temperature Chemical Physics,” Reinhold, New York, 1963. (29) Giles, R., “Mathematical Foundations of Thermodynamics,” Macmillan, New York, 1964. (30) Glasstone, S., “Thermodynamics for Chemists,” Van Nostrand, Princeton, N. J., 1950. (31) Goldwasser, S. R., IND.ENG.CHEM.51, 595 (1959). (32) Gordon, S., Zeleznik, F. J., A R S J . 32, 1195 (1962). (33) Gordon, S., Zeleznik, F. J., Tech. Note D-1737, NASA, Washington, D. C., October 1963. (34) Gordon, S., Zeleznik, F. J., Huff, V. N., Tech. Note D-132, Ibid., October 1959. (35) Hancock, J. H., Motzki?: T S in “Kinetics Equilibria and Performance of High Temperature Systems, G: S:’Bahn and E.’ E. Zukoski, Eds., p. 82, Butterworth, Washington, D. C., 1960. (36) Hanzel, P. C., Rept. No. 30-3, Jet Prop. Lab., Calif. Inst. Tech., Pasadena, April 1959. (37) Hildebrand, F. B., “Introduction t o Numerical Analysis,” McGraw-Hill, New York, 1956. (38) Hilsenrath J. Klein M Sumida D Y in “Thermodynamic and Trans ort Properties of %a;es, Li4uids: and Solids,;’ p:’416, ASME and McGraw-Hill, 8 e w York, 1959. (39) Householder, A. S., “Principles of Numerical Analysis,” McGraw-Hill, New York, 1953. (40) Huff, V. N., Calvert, C. S., Tech. Note 1653, NASA, Washington, D. C., July 1948. (41) Huff, V. N., Gordon, S., Morrell, V., Tech. Rept. 1037, Zbid., 1951. (42) Kandiner, H. J., Brinkley, S. R., IND.ENG.CHEM.42, 850 (1950). ~~
~
c
(43) Kirkwood, J. G., Oppenheim, I., “Chemical Thermodynamics,” McGraw-Hill, New York, 1961. (44) Krieger, F. J., White, W. B.,J. Chem. Phys. 16, 358 (1948). (45) Kubert, B. R., Stephanou, S E. in “Kinetics E uilibria and Performance of High Temperature Systems,” G: S.’Bahn and E. ’E. %ukoski: Eds., p. 166, Butterworth, Washington, D. C., 1960. (46) LeGrives, E., Barrere, S., La Recherche Apronnutique 68, 31 (1959). (47) Levine, H. B., J . Chem. Phyr. 36, 3049(1962). (48) Levy, S. L., Reynolds, 0. A., General Chemical Div., Allied Chemical, Morristown, N. J., May 1960. (49) Lphr, L. R., Rall, L. B., M R C Tech. Summary Rept. No. 464, Univ. of Wisconsin, Madison (available from DDC as AD-600736) March 1964. (50) Marek, J., Holub, R., Collection Czech. Chem. Commun. 29, 1085 (1964). (51) Martin, F. J., Yachter, M., IND.ENC.C H m . 43, 2446 (1951). (52) Martinez, J. S . , Elverum, G. W., Jr., Memo. No. 20-121, Jet Prop. Lab., Calif. Inst. Tech., Pasadena, December 1955. (53) McBride, B. J., Gordon, S . , Tech. Note D-4097, Washington, D. C., 1967. (54) McMahon, D . G., Roba:k, R . , in “Kinetics, Equilibria, and Performance of High Temperature Systems, G. S. Bahn and E. E. Zukoski, Eds., p. 105, Butterworth, Washington, D. C., 1960. (55) Michels, H. H., Schneiderman, S. B., in “Kinetics, Equilibria, and Performance of High Temperature Systems,” G. S. Bahn, Ed., p. 205, Gordon and Breach Science Publisher, New York, 1963. (56) Mingle, J. O., Special Rept. No. 25, Kansas State Univ., August 1962. (57) Naphtali, L. M., IND.ENO.CHEM.53, 387 (1961). (58) Naphtali, L. M., J . Chem. Phys. 31, 263 (1959). (59) Pitzer, K. S., Brewer, L., “Thermodynamics,” 2nd ed., McGraw-Hill, New York, 1961. (60) Potter, R. L., Vanderkulk, W., J . Chem. Phyr. 32, 1304 (1960). (61) Raju, B. N., Krishnaswami, C. S . , Indian J . Tech. 4, 99 (1966). (62) Rosen, E. M., in “Proceedings of2lst National Conference of ACM,” pp. 37-41, Thompson Book, Washington, D. C . , 1966. (63) Sachsel, G. F., Mantis, M. E., Bell, J. C., in “Third Symposium on Combustion, Flame and Explosion Phenomena,” p. 620, Williams and Wilkins, Baltimore, 1949. (64) Sage, B. H., “Thermodynamics of Multicomponent Systems,” Reinhold, New York, 1965. (65) Scully, D. B., Chem. Eng. Sci. 17, 977 (1962). (66) Shapiro, N. Z., Shapley, L. S., J. Soc. Ind. Appl. Math. 13, 353 (1965). (67) Shapiro,, N. Z.,, Shapley, L. S., Rept. No. RM-4464-PR, T h e Rand Corp., Santa Monica (available from DDC as AD-636605), July 1966. (68) Snow, R. H., Rept. No. IITRI-C929-3, I I T Res. Inst., Chicago, September 1963. (69) Spang, H. A., SIAM R m . 4, 343 (1962). (70) Storey, S. H., Van Zeggeren, F., Can. J. Chem. Eng. 42,54 (1964). (71) Stull, D. R., “JANAF Thermochemical Tables,” Clearinghouse, U. S. Dept. Comm., Springfield, Va., August 1965. (72) Svehla, R. A., Spec. Publ. 3011, NASA, Washington, D. C., 1964. (73) Tornheim, L., J . Arsoc. Comput. Mach. 11, 210 (1964). (74) Traustel, S., VDI Zcit. 88, 688 (1944). (75) Turner, L. R., Ann. New York Acnd. Sci. 86, 817 (1960). (76) Vale, H. J., in “Kinetics, Equilibria and Performance of High Temperature Systems,” G. S. Bahn, Ed., p. 271, Gordon and Breach Science Publishers, 1963. (77) Van Ness, H . C., “Classical Thermodynamics of Non-Electrolyte Solutions,” Macmillan, New York, 1964. (78) Villars, D. S., J. Phys. Chem., 63, 521 (1959). (79) Villars, D. S., Rept. No. N O T s T P 2354, U. S. Naval Ordnance Test Station, China Lake, November 1959. (80) Warga, J., J. SOC. Ind. Appl. Math. 11, 594 (1963). Ser. A . 241, 132 (1957). (81) Weinberg, F. J., Proc. Roy. SOC., (82) White, W. B., Johnson, S. M., Dantzig, G. B., J . Chem. Phyr. 28,751 (1958). (83) Wiederkehr, R. R. V., Rept. No. AR-IS-60, Dow Chemical, Midland, June 1960. (84) Wilkins, R. L., in “Proceedings of the Propellant Conference,” p. 315, Eng. Expt. Station, Ohio State Univ., Columbus, 1960. (85) Wilkinson, J. H., “The Algebraic Eigenvalue Problem,” Clarendon, Oxford, England, 1965. (86) Wilson, A. H., “Thermodynamics and Statistical Mechanics,” Cambridge Univ., London, England, 1957. (87) Winternitz, P. F., in “Third Symposium on Combustion, Flame and Explosion Phenomena,” p. 623, Williams and Wilkins, Baltimore, 1949. (88) Wolfe, P., Comm. ACM. 2, 12 (1959). (89) Zeleznik, F. J., Gordon, S., A R S J . 32, 606 (1962). (90) Zeleznik, F. J., Gordon, S., Can. J. Phyr. 44, 877 (1966). (91) Zeleznik, F. J., Gordon, S., Tech. Note D-473, NASA, Washington, September 1960. (92) Zeleznik, F. J., Gordon, S . , Tech. Note D-1454, Ibid.,October 1962. (93) Zupnik, T. F., Nilson, E. N., Sarli, J. J., Rept. No. UACRL-C910096-11 (NASA CR-54042), United Aircraft, East Hartford, September 1964.
VOL. 6 0
NO. 6
JUNE 1968
57
