E. E. Stone
Research and Development Center General Electric Company Schenectady, New York
I
I
Complex Chemical Equilibria Application of Newton-Raphson method to solve non-linear equations
Chemical equilibria involving more than a few chemical species are not easily solved. Recently there have appeared a number of schemes which depend upon the use of high speed automatic computers to obtain answers within a reasonable period of time. However, for small problems or for purposes of illustrating these techniques, programming an automatic computer may take up a large portion of available time. The procedure given below will solve equilibria iuvolving up to six different chemical species with the aid of only a desk calculator in a period of a few hours. I t follows the approach used by White, et al.' (Of course, the procedure can also he programmed for solution by automatic computer.) The method is an iterative procedure involving preliminary estimates of the equilibrium concentrations followed by calculation of correction terms which are applied to the initial estimates; and the process is repeated until the degree of improvement falls below some designated value. Both gaseous and solid species can be included in the calculations; however, for ease of explanation only an ideal gaseous system will be described. I n any chemical reaction the total mass of the reacting elements remains unchanged, which in mathematical terms is
,,
where A = number of atoms of element jin a molecule of specie i, Xt = number of gram moles of specie i, and B, = total atomic weight of j in the system. The state of equilibrium between reacting chemical species is the extent or degree of reaction such that the sum of the free energies of the individual species is a minimum value. Mathematically this situation can be stated by use of the Gibbs free energy function, which for an ideal gaseous system of n different chemical species is
where R = gas constant, T = absolute temperature, P = total pressure, n = number of moles of specie i, n = number of moles of specie i, x, = mole fraction of specie i, and fi = f(T). Or, written in slightly different form more suitable for this procedure, 'DANTZIG, G., JOHNSON, S., AND WHITE, W., Management F.J., A N D WHITE,W. B.,J. Chem. Science 5,38(1958);KRIEGER, Ph~s.,16, 358 (1948); WHITE, W. B., JOHNSON, S. M., AND DANTEIG, G. B., 3. Chem. Phys. 28, 751 (1958).
where G , = molal Gibbs free energy, P I = partial pressure of i-th specie, and GO/RT is derived from tables of the function (Go- H#)/T. From Dalton's law for ideal gases, the partial pressure of the ith component is
Substituting (4) in equation (3) and for conditions of 1 atm the result is equation (5).
Then that set of non-negative x, values which minimizes equation (5) subject to the constraints of equation (1) will be the equilibrium values for the chosen conditions. From the mass balance constraints of the system there will he in. simultaneous equations in n unknowns. These are solved form of the 2:s by applying an elimination scheme such as Gauss's or Cholesky's; the resulting system of (n-m) unknowns is then substituted into equation (5). A minimum value of the function will occur where b+/bx, = 0. Now the problem becomes one of finding the roots of these (n - m) unknowns; the Newton-Raphson2 method is convenient t o use and lends itself readily to either hand or machine computation. After becoming familiar with the method it is possible to write the final terms by inspection of the original expressions. The method proceeds as follows: Let 6,' = b+/bz, = 0, xi0 = the approximate roots of the equations and hc = correction values such that x, = xO , ht or (in mathematical notation) these (n-m) equations are +'(XI.. . zn-,) = ~'(zP h,, . . . zon-, h,-,) (6) This equation is expanded according to Taylor's theorem and because the value of h , is small relative to that of st,the square and higher power terms can be neglected, resulting in the system
+
+
+
SCARBOROUGH, J. B., "Numerical Mathematical Analysis," The John Hopkins Press, Baltimore, Maryland, 1955, p. 203. Volume 43, Number 5, May I966
/
241
Initial values of x, are estimated and the h, corrections found by solving the system of (n-m) simultaneous equations using determinants or matrix algebra. Illustrative Example
While all texts in chemical and chemical engineering thermodynamics discuss chemical equilibria, they seldom go beyond the simple systems of two atoms or else assume the simultaneous side reactions are of negligible effect. For example, we may choose one of the important industrial processes in the manufacture of synthesis gas (amixtureof hydrogen and carbon monoxide) by the catalytic reforming of methane, the primary constituent of natural gas: CHI
+ HzO = 3Hz + CO
However, it is also possible to have other reactions occurring simultaneously, for example:
and then the products of these reactions can interact in several ways such as:
Each of these chemical reactions is reversible and what is observed is the net effect between the forward and reverse reaction. The extent or degree to which any reaction will proceed is a function of the thermodynamic properties of the reacting constituents. The figure illustrates the variation of the equilibrium
constant with temperature for seven possible reactions involving these same chemicals species. A more positive value of the equilibrium constant indicates a greater degree of conversion for the reaction. The amount of HI or CO that can be obtained at any temperature will be dependent not only on the reaction of primary interest, but also on the interaction between these chemical species via the "side reactions." For this particularsystem these reactionsarenot all independent and the overall system can be reduced to two equation^.^ In most other s y s t e m ~ a n d in this system under different conditions of temperature or pressureerroneous results can be obtained if all the important reactions are not included in the calculations. Such a study involving several simultaneous chemical reactions requires the application of unusual mathematical techniques. What is of interest here is that moderately complexsystems can easily be handled with no more than a desk calculator. When the Gibbs minimum free energy method is used the procedure does not require any special knowledge of the chemical system involved, but is it necessary to include all possible specie until computation indicates certain concentrations that are insignificant and can then be excluded from the mixture. A recent textbook example serves as an illustration. It involves the reaction of 1 mole of CHa with 5 moles of H 2 0 at 6OO0C, considering the species a t equilibrium to be CH4, HzO,COz, CO, and The A , table for these is the following: i
CH4
H1O
CO*
H1
CO
Bj
+
Values of B1 come from 1 CHI 5 HzOwhich are the atoms of reactants chosen for this case adjusted to a total of 100 grams i.e.: Substance
No. atoms
Atomic wt.
Grams
Total grams 106.11973
or per 100 grams the quantity of hydrogen is
and similarly for the C and 0 atoms. From the above mass balance table the constraint equations are XI
421
+ 2x1 xl
+ + 2.
z4
+ 2x8 +
z4
=
0.94233
(8)
+ 22%= 13.19265
(9)
4.71165
(10)
=
which are solved in terms of (n - m) = 5 - 3 = 2 unknowns. Any 2 of the xi's can be selected; for this
Figure 1.
242
/
nmp.mhul. 'K
Varimtion of equilibrium cmrtant with ternperdure.
Journal o f Chemical Education
3 See DODGE, B. F.,"Chemical Engineering Thermodynamics," MeGrrtw-Hill Book Co., Inc., New York, 1944, p. 526; or HouGEN, 0 . A., AND WATSON, K. M., "Chemical Process Principles," John Wiley & Sons, Inc., New York, 1947, p. 733. 4 The choice of using 5 males of H 2 0 per mole of CH4 was purely arbitrary and any value above about 3.5:l could have been used. Below an H>O:CH, ratio of 3.5: 1, solid carbon exists in the system and this introduces a modified calculating procedure (mathematically it is possible to calculate values for any H30:CH4ratio).
illustration xz and xa were chosen with the result shown below: 2,
st
= 0.94233 =
4.71165
- 2xr sa - z, = zn
from e q .
+ xa + xr
f r o m eq.
(8)
- 3x2 - 228
f r o m eq.
(9)
-3.76932 zs = 6.59632
- 22, -
(10)
2%=
14.13496
When the equivalent expressions involving only xz and x8 are substituted in the above equation and differentiated with respect to each of these unknowns the result is
Then the following matrix can be written G"/RT (600 C )
xi
The free energy term Go/RT is obtained from interpolation of free energy functions as illustrated in Table 1 (600°C = 873°K). From the above matrix the minimizing functions given in eqns. (12) and (13) can be written down by inspection, keeping in mind the xt with which the partial derivative is being taken. The constant is obtained as the sum of the (G8/RT)I multiplied by the appropriate x, coefficient while the in term comes from the individual x, expressions raised to that Dower determined bv the coefficient of the ~ a r t i c ular bx; term. However, h order to undemtani these manipulations it will help to follow the intermediate steps.
Table 2.
Specie
AH10 (oal/moie)
Free Energy Functions"
8W°K
T 900°K (oal/dep mole)
873'K
C/RTb 873'K
'From National Bureau of Standards Selected Values of Chemioal Thermodynamia Properties Series 111. b AG/T = C" - HO/T + ~ F I J O / T R. 1.98719eal/degmale.
-
The Taylor series expansions for these equations can also be written down by inspection. (The reader will
Detail of Calculations
Sequence
number
Table 1.
Procedure*
First iteratmn
Second iteration
Third iteration
Est. value of z* Est. value of zs ( 1 ) ( 2 ) - 3.76932 = ( E s t . a f 2,) 4.71165 - ( 1 ) - 2 ( 2 ) = ( E s t . of
+
31
32
(30)
+ ( 2 ) = improved xr
* N u m b e r s i n prentheses refer to sequence number. Volume 43, Number 5, May 1966
/
243
understand better the ease with which the final forms can be obtained after once performing the manipulations.) '*a
a=*
=
'
+
+ + +
1 - 3 . 7 6 ~ 3 2 zs ZI 1 9 4.71165 - 2 2 - 22. + 14.13496 - 322 - 22. 4 15.07729 - 2x3 - 228 ( 1 4 )
x2
*az3= w = -3.76932 1+ zp + z3 + axx
The solution of these equations is facilitated by a systematic procedure of bookkeeping as shown in Table 2, where solution by determinants was used. Initial values of x, and x3 (3.5 and 0.6) were selected (see sequence numbers 1 and 2) so that the remaining x, values (x, x,, and xs) will have positive values. Should negative numbers arise then the initial guesses must be adjusted to such values as will correct the situation. With each iteration the values of &' and &'converge toward zero (see sequence numbers 13 and 16). The first corrections for xzand x3were determined to be -0.41072 and 0.18562; the adjusted x, values obtained after three iterations were S~ecie
Moles st Eauilibrium
The author wishes to thank the General Electric Company for its permission to publish this work.
244 / lournol of Chernicol Education
