Ind. Eng. Chem. Fundam. 1983,22, 216-219
216
Identification of Principal Equilibria in Complex Systems with Linear Programming Robert A. Alberty Department of Chemistry, Massachusetts Institute of Technology, CambrMge, Massachusetts 02 139
Linear programming provMes the means for rapidly surveying the effects of temperature, pressure, and elemental composition on the chemical equilibrium of a system involving many chemical species. This article provides an interpretation of the results of linear programming calculations for chemical systems and describes a procedure for investigating complex systems. Linear programming picks out the products of the reactions that will give the lowest Gibbs energy within the constraints on the system. The results of linear programming calculations are compared with the equilibrium compositions for systems containing C(graphite), CO(g), CO,(g), CH,(g), H2(g),and H,O(g).
Introduction In 1947, Dantzig developed the simplex method for solving the general linear programming problem (Dantzig, 1963). This method is widely used in management science for finding the optimal solution of linear problems where a quantity, such as cost of production, is to be minimized subject to linear constraints on material and equipment (Bradley et ai., 1977). Computer programs have been developed for the routine solution of these problems with very large numbers of variables and constraints. Smith and Missen (1969) recommended linear programming as a means for finding a suitable starting composition for the more dificult calculation of the equilibrium composition. By use of linear programming, the Gibbs energy of the system is minimized subject to the elemental abundance equations N i=l
akini = bk
(k = 1, 2, 3, .. ., M)
(1)
and ni > 0, where ni is the number of moles of species i, aki is the subscript of the kth element in the molecular formula of species i, and b is the number of moles of element k in the system. The system contains M elements and N species. Each species in a different phase counts as an additional species. Linear programming provides the means for minimizing the Gibbs energy, neglecting niIn (ni/n)terms in the equation for the chemical potentials of ideal gases and ideal solid solutions. N
G = C ni[AGofi+ R T In ( P / P o )+ R T In ( n i / n ) ] (2) i
Here AGOfiis the standard Gibbs energy of formation of species i, P is the total pressure, P is the standard-state pressure, and n is the total number of moles of perfect gas in the system. The R T In ( P / P o )term is added only for gaseous species. In order to make a linear programming calculation, it is necessary to specify the chemical species present in the system at equilibrium, the Gibbs energies of formation, that is, AGfi + R T In ( P / P o ) ,for these species at the temperature and pressure of the system and the composition of the system in terms of elements. The linear programming calculation yields the number niof moles of the species which will give the lowest value of G (defined by eq 2 without the niRT In (ni/n)terms) within the element-balance constraints. Any possible chemical change in a system may be represented by a set of R linearly independent chemical equations. The number R is given
by R = N - rank (A) where N is the number of species and A is the system formula matrix (Smith and Missen, 1969). The rank of the system formula matrix is ordinarily equal to the number of elements. However, the rank of the system formula matrix may be smaller; for example, if some elements are in a constant ratio in all the species in which they occur, it will be smaller. The rank of this matrix is larger if additional linear constraints on the numbers of moles of various species arise from the mechanism (Bjornbom, 1975). For a given set of conditions, the linear programming solution gives the numbers of moles of products of a single chemical reaction which will give the lowest Gibbs energy if the reaction goes to completion to produce unmixed products. The linear programming calculation is not concerned with the reactant species, except through their overall elemental composition. However, it is frequently desirable to think in terms of specific numbers of moles of reactant species, and in this case the calculation yields a balanced chemical equation and its Gibbs energy change. When the calculation is made at the standard-state pressure, the equilibrium constant obtained is K the constant written in terms of P/P", where Po is tKe standard-state pressure. When R T In ( P / P o )is added to AGOfi,the resulting equilibrium constant is Ky, which is expressed in terms of mole fractions y i in the gas phase at the total pressure P. In summary, the linear programming solution yields the number of moles of specified products that will give the lowest Gibbs energy for the system assuming that the reaction is complete and that entropies of mixing of products may be neglected. This solution may be interpreted as a balanced chemical equation if a set of reactants with the elemental composition of the system is placed on the left. The number of possible reactions goes up rapidly with the number of species, but linear programming provides an extremely powerful means for identifying the products with the lowest Gibbs energy that satisfy the elementbalance constraints. However, the utility of a linear programming calculation is not over at this point. If the reactant species are withdrawn one by one from the calculation, the calculation can be repeated to identify the products with the next smallest Gibbs energy. This process can be continued to identify further reactions until all possibilities have been exhausted. System Considered The system C(graphite), CO(g), H k ) , COdg), H@(g), and CH,(g) was chosen because equilibrium compositions
0196-4313/83/1022-0216$01.50/00 1983 American Chemical Society
Ind. Eng. Chem. Fundam., Vol. 22, No. 2, 1983
Table I. Reactions Obtained by Linear Programming at 1 Atm 588.7 K, H / O = 1
Table 11. Reactions Obtained by Linear Programming at 1 0 0 atm 588.7 K, H i 0 = 1 K, = 7 . 7 1 x 108 K, = 6 . 1 2 x lo8 K, = 1 . 7 9 x 108
1. 2. 3.
217
19. '/qC 20. C
+
+
l/zH,O
'/,H,O
+
+ '/402 '/do2
21. '/ZC + 1/2Hz0 +
=
'/402
= '/zH,O
+ '/,CO,
+ '/zCO, = '/,HZ + '/zCOz I/zCH,
H/O = 2 22. H , O = H , O 23. C + H,O = '/zCO, + 'IzCH, 24. '/,C + H,O = H, + '/zCO,
K, = 7.69 x 108 Ky =
x
6.12
loa
Ky =
5.64
x 10' K, = 1 K, = 0 . 6 3 K, = 5.38 x 10-3
4. 5. 6.
K, = 1 K, = 0 . 6 3 3 K, = 5.37 x
7. 8. 9.
K, = 3.43 K, = 2.17 K, = 1
H/O = 3 25. '/4C + H,O + '/zHz = H,O + '/,CH4 26. l'/.+C + H,O + '/,HZ = '/zCO, + 3/4CH4 27. H,O + '/zH2 = H,O + '/,HZ
K,
28. C
K, = 10.85 K , = 6.85 K, = 1
1588.7 K, H/O = 1
4.37 x 10s K, = 8.82 x 104 K, = 7 . 8 9 x 103
10.
11. 12.
13. C
+
+ H, + H,O = CO +
H/O = 2
H,O = CO
'IzCH, 14. l'/,C 1 5 . '/zC + HzO = '/zCO, + H, 16. C
+ H,O +
1 7 . 1/,C 18. 13/&
=
Hi0 = 3 1/,H, = CO + 3 / ~ H z
+ H,O + I/zHz = ' / S O 2 + VZH, + H,O + I/,H, = CO + 3/4CH4
K, = 1 . 0 5 x 103 K, = 4 2 . 8 K, = 1 8 . 9 5 K, = 1.05 x 103 K, = 1 9 . 0 K, = 8 . 6 3
have been calculated over wide ranges of temperature, pressure, and initial composition by Baron et al. (1976), assuming ideal gas behavior. They obtained the equilibrium compositions by selecting three independent equilibrium equations and solving them by an iterative method. Since there are 6 species and 3 element-balance equations, there are three linearly independent reactions. Possible chemical changes in this system may be represented by quite a number of sets of three independent chemical reactions. However, in order to represent the equilibrium composition of the system, the question arises as to whether some particular set might be more convenient and whether a single reaction or two might actually represent chemical change in the system within experimental error. The choice of such a reaction, or reactions, obviously depends on the temperature, pressure, and elemental composition of the system. It is obvious that the choice of a reaction, or reactions, to represent the equilibrium state of the system will depend on the temperature and pressure, but it is less clear how this choice depends on the elemental composition of the system, in this case the H/O ratio. For the system we are interested in here, only H and 0 have to be considered since there is always an excess of graphite. If we assume that the system can be represented by a single reaction that goes to completion, this must be a reaction in which the left and right sides have the correct elemental composition in terms of hydrogen and oxygen. Results of Calculations Linear programming calculations were carried out for the system C(graphite1, CO(g), H&), COdg), H2O(g), and CH4(g)at 588.7 and 1588.7 K (these temperatures correspond with the 600 and 2400 OF of Baron et al., 1976), 1
29. 30.
+ '/zH,O t '/4OZ = CO + '/,HZ 1'/4C + '/zH,O + = CO + l/qCH4 */,C + '/zH,O + = '/zCO + '/zH20 '/402
'/402
K, = 1 . 3 8 x 104 K , = 8.83 x 103 K y = 4.27
x 103
H/O= 2 31. C + H,O = CO + H, 32. l'/zC + H,O = CO + '/,CH4 33. '/zC + H,O = '/zCO, + H,
K, = 1 0 . 5 K, = 1 . 4 5 K, = 1 . 9 0
Hi0 = 3 34. C + H,O + '/,HZ = CO + 3/,H, 35. lY4C + H,O + 1/,H, = CO + 3/.,CH4 36. I/,C + H,O + '/%HZ= l/,CO, + 3/,H2
K , = 10.5 K, = 2.73 K, = 1 . 8 9
atm and 100 atm, and H/O = 1 , 2 , and 3. The standard Gibbs energies of formation that were used for CO(g), C02(g),H,O(g), and CH4(g)were -163.60, -400.62, -214.61, and -24.13 kJ mol-' at 588.7 K and -251.24, -396.42, -159.35, and 84.56 kJ mol-' at 1588.7 K. Tables I and I1 give the first three reactions obtained under each set of conditions, along with their K y values. Linear programming yielded 5-6 reactions for each set of conditions, but only the three with the largest Ky values are given in Tables I and II. Some of the 5-6 equations are not linearly independent because they are obtained for different systems. The initial system has 6 species and 3 linearly independent reactions, but only the one with the largest equilibrium constant is taken. Then the system is changed by leaving out one of the species on the right-hand side of the first equation. Now there are two linearly independent reactions, but only the one with the largest equilibrium constant is taken. This process is continued until all the species are used up. Since the calculations are made in this way, reactions such as numbers 4 and 9 in Table I are obtained; these reactions indicate that according to linear programming no chemical change can take place in the system which will lower the Gibbs energy, of course ignoring the Gibbs energy of mixing. The program used here (see Acknowledgment)yields the Gibbs energies of the products, but K y is readily calculated by subtracting the Gibbs energy of the assumed reactants. At the three H/O ratios the reactants are H/O = 1: &(graphite)
H/O = 2:
+ '/&O(g)
+
nC(graphite1 + H20(g)
H/O = 3: nC(graphite1 + H20(g)+ '/2H2(g) The number n of moles of carbon varies because an excess
218
Ind. Eng. Chem. Fundam., Vol. 22, No. 2, 1983
Table 111. Equilibrium Mole Fractions Calculated by Baron, Porter, and Hammond (BPH) and by Linear Programming
co
HZ
COZ
H2O
CH,
0.02
0.47 0.66
0.07
0.04
0.42 0.33 0.66 0.21
0.05
0.5 0.12
1 atm
588.7 K
Hi0
= 1
2 3 1588.7 K
H/O = 1 2 3
BPH rxl rx2 BPH rx4 rx5 BPH rx7 rx8 BPH rxlO BPH rxl3 BPH rx16
0.56 1.00 0.52 0.80
0.40
0.33 0.19 0.5 0.29 0.20 0.60
0.33 0.33 0.50 0.50 0.60 0.60
0.66 0.66 0.50 0.50 0.40 0.40
100 atm 588.7 K
Hi0
= 1
2
3 1588.7 K
Hi0 = 1 2
3
BPH rxl9 rx20 BPH rx22 rx23 BPH rx25 rx26 BPH rx28 BPH rx31 BPH rx34
0.01 0.01 0.66 0.66 0.49 0.50 0.40 0.40
of C(graphite) was always present, so that various amounts of this reactant appear in the products. Since the results are compared with the equilibrium calculations of Baron et al. (1976), the standard Gibbs energies of formation at 577.6 and 1577.6 K were calculated using their empirical expressions for these quantities as a function of temperature. They had determined these empirical expressions from the JANAF tables (Stull and Prophet, 1971). The equilibrium mole fractions of the reacting gases, as calculated by Baron et al. (1976), are given in Table 111. These equilibrium values are compared with what is expected from the reactions in Tables I and 11, results which are given in the rows labelled “rx”. At 588.7 K, 1 atm, and H/O = 1, the first reaction has a very large equilibrium constant and would yield ycoZ= 0.33 and Y H ~ O= 0.66, in comparison with yco = 0.42 and yHp = 0.47 calculated by Baron, Porter, and fiammond. However, the equilibrium constant for the second reaction in Table I is smaller by only a factor of 1.26, and it would yield ycoz = 0.66 and ycb = 0.33. The equilibrium values are of the magnitude expected from the first reaction alone, but the second reaction does have an effect and helps to explain the equilibrium values. The third reaction in Table I produces hydrogen, but very little hydrogen is formed because the equilibrium constant for this reaction is smaller than the first by a factor of 4.3. At 588.7 K, 1 atm, and H/O = 2, the reaction with the largest equilibrium constant is that forming H20, which corresponds to no reaction at all. The equilibrium mole fraction of H 2 0 is 0.56, however, rather than unity. This is because the equilibrium constant for the second reaction is 0.633. A t 588.7 K, 1atm, and H/O = 3, the reaction with the largest equilibrium constant yields yH2O= 0.8 and yCH4= 0.20, in comparison with the equilibrium values of 0.52 and 0.29, respectively. Again the equilibrium constant for the
0.30 0.33 0.44 0.50 0.53 0.60
0.42 0.33 0.66 0.22 0.50 0.12 0.33 0.50 0.01
0.49 0.66 0.56 1.00 0.54 0.66
0.08 0.33 0.22 0.50 0.32
0.02
0.50 0.02
0.01
0.02
0.03
0.01
0.02
0.05
second reaction is close, and helps to explain the equilibrium composition. A t 1588.7 K and 1 atm, the first reaction obtained by linear programming and given in Table I represents the equilibrium composition at each H/O. The equilibrium constant for the second reaction in each case is enough smaller so that competition with the first reaction does not have to be considered. The equilibrium constants expressed in terms of mole fractions at 100 atm are summarized in Table 11. These values were calculated by adding R T In 100 to the stand’md Gibbs energies of formation of the gaseous species. The same sequence is found for the first three reactions at both pressures, except for the third reaction at 1588.7 K at H/O = 1, and for the second and third reactions at H/O = 3. The comparison of the linear programming results with the equilibrium mole fractions is shown in Table 111.
Discussion On the basis of these linear programming computations the equilibrium composition may be obtained more accurately by calculating the extent of reaction from the principal equilibrium expression and then introducing the second equilibrium expression to obtain the mole fractions of further species. However, it is not the objective of this article to pursue the complete calculation of the equilibrium mole fractions. Linear programming does not provide the equivalent of the various programs for solving this problem, but rather provides a rapid and economical way of surveying wide rang= of conditions for complex systems. The above calculations have been carried out for integer ratios of H/O, but calculations can be made for any ratio, even though the results may not be as easily interpreted in terms of equilibrium constants. The program used provides not only optimal solutions but also upper and lower limits on the objective coefficients (Gibbs energies of formation in this case) and constraints
219
Ind. Eng. Chem. Fundam. 1983, 22, 219-223
agement science. The program used here was implemented on a minicomputer. Linear programming is not a substitute for the algorithms (reviewed recently by Smith, 1980) that have been developed for the computation of equilibria in complex systems when precise compositions are required and there are adequate resources and staff. However, linear programming provides a rapid and efficient way of surveying wide ranges of independent variables and getting a better intuitive understanding of the principal reactions. Acknowledgment The author is indebted to Professor Thomas L. Magnanti of the Sloan School of Management at MIT for making available a linear programming routine on their Prime computer. Financial support was received from the Dreyfus Foundation. Literature Cited
(moles of various elements present) and yields "shadow prices" for changes in the constraints. The "shadow prices" are the changes in Gibbs energy when the number of moles of an element in the system is changed by one unit, everything else being held constant. Thus a great deal more information is obtained from the linear programming calculation than the optimal solution. When it is of interest, Gibbs energies and numbers of moles of various elements can be changed to values beyond these limits to see what new optimal solution is obtained. Conclusions These calculations show that linear programming provides a useful tool for surveying the effects of temperature, pressure, and reactant composition on the equilibrium of a system for which the number of species is greater than one larger than the number of linear constraints. In the case examined there are 6 species and 3 linear constraints so that three independent chemical reactions are required to represent all possible chemical changes. However, under some conditions the equilibrium composition may be represented by a single chemical reaction, and in others by two chemical reactions. These reactions are conveniently obtained from linear programming. Linear programming may be used for systems with a very large number of species requiring a very large number of independent chemical reactions. The computer program used will handle 400 species and 100 linear restrictions. Linear programming is frequently available in connection with large computer systems because of its use in man-
Baron, R. E.; Porter, J. H.; Hammond, 0. H. "Chemical Equilibria in CarbonHydrogen-Oxygen Systems"; The MIT Press: Cambridge, MA 1976. Bjornbom, P. H. Ind. Eng. Chem. Fundam. 1975, 74, 102. Bradley, S. P.; Hax, A. C.; Magnanti, T. L. "Applied Mathematical Programming"; Addlson-Wesiey Publishing Co.: Reading, MA, 1977. Dantzig, G. B. "Linear Programming and Extensions"; Princeton University Press: Princeton, NJ, 1963. Smith, W. R.; Missen, R. W. Can. J . Cbem. Eng. 198% 46, 269. Smith, W. R. I n "Theoretical Chemistry, Advances and Perspectives"; Academlc Press: New York, 1980 Voi. 5, p 185. Stuil, D. R.; Prophet, H. "JANAF Thermochemical Tables", 2nd ed.;NSRDSNSB 37, U S . Government Printing Office, Washington, DC, 1971.
Received for review March 11, 1982 Accepted January 3, 1983
Equilibrium Constant for the Methyl tert -Buty I Ether Liquid- Phase Synthesis by Use of UNIFAC F. Colombo, L. Corl, L. Dalloro, and P. Delogu" Centro Ricerche MONTEDIPE, Via San Pietro 50, Bollate, Milano, Italy
A predictive method is used for evaluating the equilibrium constant of a reaction carried out in the liquid phase. This method uses the values of the standard free energy of formation available in the literature for the chemicals involved and the UNIFAC estimates of activity coefficients to describe the liquid-phase nonideality. Some experimental values of the equilibrium concentrations were obtained in a micro-pilot plant. The predicted equilibrium constants agree quite well with our experimental values.
Introduction Methyl tert-butyl ether is a chemical whose production will greatly increase in the near future because of its good antiknock properties. It is industrially manufactured from the isobutene contained in olefinic C4 cuts and from methanol through the reversible reaction i-C4H8 CHSOH * t-C,HgOCH3
+
This reaction is exothermic; therefore the equilibrium is shifted to the right at low temperatures and to the left at high temperatures. In all industrial processes the reaction is carried out in the liquid phase, but no data on the equilibrium constant are available. Classical methods of thermodynamics allow the computation of the gas-phase equilibrium constant; the derivation of the liquid-phase constant, in terms of mole fractions, requires knowledge of the liquid-phase activities 0196-4313/83/ 1022-021 9$01.50/0
for all components of the reaction mixture. The reaction is usually performed in the presence of inert linear C4 olefins; therefore the system is a multicomponent one, and its behavior is not simple to predict. We have tried to evaluate the equilibrium constant only by means of predictive methods, using available literature values of AGf" for the chemicals involved and UNIFAC predictions of activity coefficients (Fredenslund et al., 1977) to describe the liquid phase nonideality. The results are compared with the experimental values measured in our laboratory. Thermodynamic Framework Tabulated values of the standard free energy of formation (AGof,298)iallow the computation of the standard free energy of reaction for the MTBE synthesis in the gas phase n
(AG
0 1983 American
298)reaction
=
Chemical Society
2 v i (AG i=l
f,298)i
(1)
