1 Engineering
Approaches
I
Calculate Complex Chemical Equilibria . . . with this technique, based o n minimizing f r e e energy
A
.)
I
LEONARD M. NAPHTALI Polytechnic Institute of Brooklyn, Brooklyn, N. Y.
TECHNIQUE for computing complex equilibria is presented which is well adapted for machine computation. I t has no limitations to systems in which few components predominate. Nonideal gases or condensed phases are easily included. I n a test case on a 10-component gaseous system, it computed equilibrium composition in 15 seconds on the IBM-704. The problem of chemical equilibrium computation is essentially one of finding a minimum of a function of many variables (normally free energy), subject to constraints, by algebraic manipulation, evaluating partial derivatives with respect to new variables, and equating the partial derivatives to zero. There are three variations on the conventional approach of setting derivatives equal to zero and solving simultaneous equations. The oldest is the conventional equilibrium constant method. The equilibrium equation is derived by considering a small amount of reaction (in which the stoichiometric equation embodies the material balance on the elements) and setting dG = 0. I n the complex case, this results in one nonlinear equation per reaction (see Equation 8). A more recent variation (4) uses the method of Lagrange multipliers to combine a free energy equation with the conservation of the elements constraint. This method also makes use of an algebraic simplification of free energy expression which is dependent on ideal gas expression for chemical potential. The last variation is the conventional (2). This minimization approach method also uses a single expression for free energy of the mixture, algebraically incorporating the material balances, and sets the derivatives equal to zero. In both of these methods (2, 4), simultaneous equations to be solved are fewer in number than the chemical equilibria, and therefore simpler to solve. The present method is based on minimizing free energy by the gradient method. This process is greatly simplified by special relationships which hold for chemical equilibrium. We can describe the change in free energy for chemical react5ons by a AG for each reaction. If we proceed forward in those reactions for which AG is negative
and backward in those for which AG is positive, a decreasing free energy is assured. Then evaluate the composition changes that result and start again. When combined, these steps lead to an algebraically very simple procedure. The result is suited to automatic computation. I t does not require ideal gases; can have any number of phases; and can even be modified to allow computation of equilibrium under constraints other than T a n d P,such as H a n d P for adiabatic flame temperatures, and S and P for isentropic expansions.
Theory The theory leading to computation of the direction of composition change for which AG 5 0 was published previously ( 3 ) . ,An account is given of it here, which -while not simpler-avoids the use of matrix algebra. A table of stoichiometric coefficients is formed in which each column represents a reaction, and each row a component. The entry is the coefficient in the equation, with a minus sign for compounds consumed, and a plus sign for those produced. The entry in row i, columnj, is vi! = (ani)/(dEi) (1) = where ni = moles of component i and extent of reaction j . A composition change is the sum, over all reactions :
C pJ2n2
AG, d&
(9)
I
Since we require that at equilibrium
dG = 0, this condition is met when each AGi = 0. These are the conventional
equilibrium equations. They are nonlinear algebraic equations, when expressed in terms of compositions; and when there are three or more of them they can become extremely difficult to solve. We propose to dodge this difficulty by minimizing free energy by a gradient method. Any composition change for which dG 5 0 is a thermodynamically possible one, and a series of such changes inevitably will terminate at the minimum point, or equilibrium composition. Let us choose d& = -AGidA where X is a parameter which determines how far we g o in the chosen direction. Then we have : dG =
-
{F
(AG,)s)dA
dG 2 0 if dX 2 0 For this direction of change, Equation 2 :
(10)
or
from
Substituting Equation 8 :
Interchanging the order of summation:
and factoring out of the sum overj in which K is fixed :
According t o thermodynamics: G =
dG =
(3)
(dni)/(dX) =
-
T M( F~
vijvhi)
(14)
erwk
(16)
At a fixed T and P: dG =
pz dn,
(4)
1
Then (dn,)/(dX) = -
A table of values of ezk can be formed from Equation 15, and the table vi,. It will have a column and a row for each compound and will be found to be symmetrical-ie., = ent. The table of values of erk depends only on the system being studied and is computed only once. Then by Equation 16, can we compute a direction of decreasing free energy from any composition. From Equation 4 we get:
It is obvious that by definition:
wzl= AG,
k
(dG)/(dX) = (8)
where AGI is the free energy change for reaction j . Hence, substituting Equation 8 into 7:
C
~ , [ ( d n z ) / ( d X ) l (17)
and from Equation 10:
VOL. 53, NO. 5 ' a
MAY 1961
387
W e can then compute the sum of the squares of the AGj’s for all the reactions by Equation 18 to compare with the accuracy of the free energy data to determine whether we are close enough to equilibrium. Near equilibrium each individually should be near zero, within the precision of the data, and hence the sum in Equation 18 should be near zero.
Computing Method and Step Size This theoretical development does not constitute a ‘komputational algorithm unless certain questions of programming strategy are resolved. The methods suggested here are adequate, but may well be improved significantly. If we are not close enough to equilibrium, a change in composition in the indicated direction is made. But what size step should be taken? For a first estimate, take a step so that the greatest decrease in moles is 10%. At all later steps the step size is computed by comparing the current direction with the preceding one. At each step, if any value of n, 5 0 or if the free energy is not lower than previously computed value, then step size is reduced by a factor of 4, and the composition is returned to its previous value. After the first step, the step size is chosen by the following technique. The values (dn,/dX) are regarded as components of a vector (which has geometric meaning for three components, but is an abstract idea in general). The cosine of the angle between two successive directions is given by:
-&(2)‘(2)”
cos
a: =
dF ($y
(I9)
This may be derived by considering the dot product of the unit vectors in each direction. I t is clear that when no more than three components are present, the geometry tells us that the Equation 19 is a cosine, and hence : 15 cosa: 5 1 (20) Inequality (Equation 20) holds generally and is Schwartz’s inequality ( 7 ) . Making use of the geometric analogy, we interpret a value of cos cy 5 0 as implying an acute angle between successive steps. This is undesirable (zigzagging) and results from too large a step size. I n this case the step size is reduced by a factor of 4. If cos a: N 1, the steps are too small, since this indicates that two successive steps are being made in substantially the same direction. Whenever 0 < cos cy 5 1, we multiply step size by the factor F , given by : F = 0.5 1.5 cos4 a: (211
-
+
F ranges from 0.5 when cos cy = 0 to 2.0 when cos CY = 1. Thus the step size is varied, in accordance with the local character of the surface.
388
Computation Experience To test the computation scheme herein described, a general subroutine was written in Fortran I1 for the IBM-704 which uses the ideal gas formula for computing chemical potential. A numerical example from the literature ( 4 ) involving ten components and seven reactions was computed by the present method. 4 n answer was obtained, correct within the precision specified (which was the precision of the data) in 15 seconds. hnother program for the basic IBM 650: written in Fortransit, ran typically for 15 minutes for hydrocarbon-oxygen combustion. The running time quoted only indicates the order of magnitude to be expected, and will depend on the choice of computational rules and their effectiveness, as well as on the size of the error tolerated, the closeness of the original guesses, and probably the choice of reactions to represent the system. Condensed Phases The inclusion of condensed phases is simple. For each condensed phase a new series of species is introduced representing those in the condensed phase. For example, if carbon is supposed to exist in the solid and vapor phases C (s) and C (9) are treated as two different species and the reaction C (s) e C (g) is included with the others. Since the chemical potential of a pure liquid or solid species does not become infinite as its amount becomes zero-i.e., the logarithmic term in the expression for pi is missing-the restriction to positive amounts is not needed to prevent blowup of computations. A solution with negative moles of a condensed phase merely indicates that if supplied the material would have been consumed to the indicated extent. Constraints Other than T and P Minimization of free energy is thermodynamically the correct criterion for equilibrium only at constant T and P. At constant H a n d P (as in an adiabatic flame temperature calculation) entropy should be maximized. At constant S and P (as in a nozzle or turbine expansion) enthalpy should be minimized. We may proceed in an analysis analogous to that for constant T and P. Instead of Equation 4 we have: dH = T d S f VdP
+
pidni i
(22)
At constant H and P we have : dS = - ( l / T )
i
pidni
(23)
And a t constant S a n d P we have: dH = Z pidni (24) Hence the direction given by Equation 17 is one for which
(2)
INDUSTRIAL AND ENGINEERING CHEMISTRY
T.P
IO;
(E)
HIP
2 0 ;
and Thus for a small enough X each function goes toward its equilibrium value, for the corresponding constraints. However, after each step, the value of pi is computed differently for each set of constraints, and the criteria of improvement vary-Le., 6G < 0 for T , P constant; 6s > 0 for €€, P constant; and 6H < 0 for S, Pconstant. The chemical potential pi is normally given as a function of temperature pressure and composition. This is convenient for computing equilibria at constant T and P, but when the constraint is constant H a n d P the chemical potential should be a function of H, P and composition (which is not usually practical) or T should be a function of H, P, and composition. The latter course is most promising. Temperature should be determined by solving the equation : H( T , Pi ni) = H( To, Pa; nit) (25) for T . A rough calculation of the change in temperature can be obtained from: A T - [(Z RiAn)/(ZCp,ni) 1 (?6) and refined by iterating in Equation 27 : Tn+1 Tn - [ ( H - Hb)/(ZCpini)l (27) obtained by Newton’s method from Equation 25. Thus after each increment in composition the temperature computed from Equations 26 and 27 is used in determining p2 to evaluate whether an improvement was obtained, and the next step in composition. In computing equilibrium at constant 5’ and P the chemical potential should ideally depend on S, P? and composition. Usually it is not available in this form and a procedure analogous to that discussed for constant H and P applies, using instead of Equation 25
S ( T , P; ni) = S(Ta, Po:
nia)
= SO
(28)
and instead of Equation 26 : AT
-[(ZSiAni)/T(ZCpni)]
(29)
and instead of Equation 27 :
literature Cited (1) Bodewig, E., “Matrix Calculus,” p. 3, Interscience Publishers, Inc., New York, 1959. (2) Dobbins, Thomas O., Thermodynamics of Rocket Propulsion and Theoretical Evaluation of Some Prototype Propellant Combinations, WADC TR-59757, Wright Air Developmcnt Center, Air Research and Development Command, U. S. Air Force, Wright-Patterson Air Force Base, Ohio, 1959. (3) Naphtali, Leonard, M., J . Chem. P h y ~ .31, 263-4 (1959). (4) White, W. B., Johnson, S. M., Dantzig, G. B., Zbid., 28, 751-5 (1958). RECEIVED for review March 18, 1960 ACCEPTEDJanuary 19, 1961 Meeting-in-Miniature, ACS, New York, March 11, 1960.
