Simultaneous Chemical and Phase Equilibrium Calculation Robert V. Sanderson' and Henry H. Y. Chien Monsanto Co., St. Louz's., M o . 65166
A new, generalized method is presented for the calculation of phase and chemical equilibria. The method requires specification of reactions and chemical equilibrium K-values. Solution of equilibrium, material balance, and phase relationships i s b y the nonlinear equation-solving method of Marquardt. The solution i s iterative, but i s very fast. The method can b e used for reactor simulation because of the solution speed. A sample problem involving an esterification reaction with phase equilibrium i s presented.
C h e m i c a l equilibrium is often a n important considerat,ion in the design of chemical reactors. Xumerous examples could be cited to show its usefulness. The equilibrium conversion of a reactant' for a given feed a t specified reaction condit'ions represent's the maximum conversion attainable regardless of reaction kinetics. If this maximum conversion is less than that economically necessary for a process, reaction conditions must be changed or further process development terminated. The design of reactors involving fast or high temperature reactions sometimes uses the concept of approaches to equilibrium. For reactions at, a high temperature, the temperature a t which the equilibrium is calculated is assumed to approach the exit temperature of the reactor but differs from it by a fixed quantity; for fast reactions, reaction extents are assumed to be a t a percentage of equilibrium reaction extents. I n some simulations, equilibrium reactions may be part, of the assumed process model. They are used to extrapolate experimental data for process design or other uses. I n many chemical processes, reactioiis occur in a two-phase system. For such a system, chemical and phase equilibria should always be computed simultaiieously to include the effect of their strong interaction. The condensation or vaporization of products from the reaction phase often changes drastically the equilibrium composition distribution. Or, in some cases, the reactants may be in different phases and the phase equilibrium relationship is a n integral part of the model. d s a good esample of simultaneous chemical and phase equilibria, the design of air pollution control equipment invariably uses favorable chemical equilibrium in a vapor liquid contacting device. The scrubbing of SOn with a n aqueous basic solution is one such example. Despite the apparent usefulness, generalized methods of solution for a rigorous model have not been reported in the literature. A recent article by Zeleznik and Gordon (1968) reviewed numerous methods for calculat'ion of chemical equilibrium. I n general, published methods of calculating the equilibrium composition of a chemical system have fallen in one of two classes. One solves algebraically the equilibrium and material balance equations; the other developed by Clasen (1965) aiid White et al. (1958) minimizes the free eiiergy of the system with geometric programming. The first is more direct for data correlation and for process simulation To whom correspondence should be addressed.
and the second is more suitable for large systems and for estimation. The algebraic solution of the mass and equilibrium relations produces the same results as the minimization of free energies under conditions as follows: First, the equilibrium constants must be defined in terms of free energies. Also, the reactions must be chosen to represent equilibrium among all chemical species present'. Kandiiier and Uriiikley (1950) outline a method for defining the maximum number of linearly independent reactions. The third condition is that rigorous phase equilibrium relationships are included in the calculations. The approach developed in this work belongs to the first of the above classes. It calculat'es the equilibrium composition of a mixture of a given temperature, pressure, and initial composition. A prior specification of possible reactions and their equilibrium constants is required. The method is designed for single-phase tems as well as multiphase systems. The method can be used for systems with one or more reactions; however, the reactions must, be linearly iiidepeiident. This new method for chemical aiid phase equilibria calculations has been programmed for the 11331 370/165 computer. It is now a block in the lfonsanto FLOWTRAK process simulation system. Previous investigators have published chemical equilibrium calculat'ion methods of varying generality. Goldwasser (1959) presented a method which solved the equations for a onephase system by direct iteration. Smith and Xissen (1968) modified a met'liod originally presented by Cruise (1964) and Villars (1959) in which the extents of reactions were searched to minimize the system free energy. Their methods were applied successfully to single-phase systems. Methods for single-phase systems do iiot work well for two-phase systems. For example when t,he highly nonlinear phase equilibrium relationships are added t'o the system, the method of searching the extents fails to converge to a solution. The problem is usually formulated as a 11011linear equation-solving problem with linear constraints of the material balances. Algorithms for this kind of problem have not been very effective. I n this paper, we reformulated the problem iiito a n unconstrained noiilinear problem which easily can be sdved numerically. T h e concept of converting constrained problems to uncoiistrained ones is not new. However, its application using penalty functions has been rather disappointing. This paper Ind. Eng. Chem. Process Des. Develop., Vol. 1 2 , No. 1 , 1973
81
demonstrates a successful application of the concept to a general class of problems.
The constraints for El’s are that
2 0
Fout i
Problem Definition and Algorithm Development
i = 1, . . . , n
(8)
m
or
For a reaction system involving the following reactions:
pij€j
2
-Fin
i = 1, . . . , n
i
(9)
j=1 n
j=1,2
...,m
1=1
where i and j are the component and reaction indices, respectively, p the stoichiometric coefficient and A i the chemical species, the reaction equilibrium constants are defined by:
The equilibrium constants, K j , may be determined experimentally or calculated by
Kj
=
exp (-AGio/RT) n
where
AGj” =
pijGio i=l
Symbols are defined in the nomenclature section. The material balance equations for the reactions are Fout i = F i n i
+
m
i
pijtj
=
1, 2, . .
. )
n
j=1
(3)
The phase equilibrium equations and the flash equations are 2. = a
pout i
i = 1,2, . . . , n
~
2 Pout
(4)
i
i=l
zi
=
CYy,
+ (1 - a ) &
i
=
1, 2, , . ., n
(5)
in the extents space. Generalized methods for solving multivariable nonlinear equations with linear constraints have been few and inconsistent’. We have attempted to find a suitable method for the above procedure. Results have been disappointing because of the nonlinearity of the problem. The linear constraints (Equations 9) are normally introduced when Equations 3 are enforced in each iteration. If, however, Equations 3 were included in the set of dependent equations and Fout were included in the set of independent variables, we would have a (m n) dimensional unconstrained nonlinear equation problem. One basic difference between this approach and a penalty function approach in removing the constraints is that this method increases the complexity of the set of independent variables whereas the other increases the complexity of the dependent variables. We believe our approach is more effective for difficult problems. The unconstrained nonlinear equations are solved using the Narquardt’s method (1963). An algorithm for solving this problem is as follows:
+
1. Estimate tj’s and Fout 2. Perform flash calculation on FOuti’s to get xi’s and yi’s 3. Repeat above steps if Equations 2 and 3 are not satisfied
T o remove the linear constraints, we have introduced to the problem two undesirable qualities :
The liquid activity coefficients, y i , may be correlated by any one of several well-known equations-e.g., Van Laar, Wilson, or nonrandom in two-liquid ( S R T L ) equations. An accurate description of the liquid activity coefficients is important for a valid solution to the problem. However, the numerical method described in this paper is independent of the correlation used. The vapor fugacity coefficients, qi, may be calculated using any one of several equations of state, e.g., Redlich-Kwong or BWR (Benedict-Webb-Rubin) equations. Equations 2-7 are solved for 4 n m 1 variables. The variables are
+ +
Fout i 2 i r Y i , zi €j
and
i
=
1, . . ., 71
i = 1, . . . ,n
T o improve the solution speed, we have developed and tested successfully the use of approximate analytical derivatives. The derivation is presented below. Numerical Considerations
j = 1, . . . , m
The material balance Equations 3 are redefined to be:
CY
This problem may be solved as a n m dimensional nonlinear equation problem. The extents ( E j ) are the independent variables and Equations 2 are the dependent equations. The algorithm may read as follows: 1. Estimate tj’s 2. Calculate FoUti’swith Equations 3 3. Perform flash calculat’ion on Fouti’s with Equat’ions 4-7 to obtain xi’s and yi’s 4. Substitute xi’s (or yi’s) into Equations 2 5 . Repeat steps until Equations 2 are satisfied 82 Ind.
1. The dimensionality of the problem has been increased significantly. This increase usually means longer computer time and larger storage size. 2. The interaction between the variables has been increased significantly. An error in the j t h Equation 2 does not affect the corresponding ti directly, but’ only i’s and then a n increase of errors causes a change in FOut in Equations 3. These errors force a change in ti. The interaction results usually in longer computer time.
Eng. Chem. Process Des. Develop., Vol. 12, No. 1 ,
1973
i
=
1,
,n
(10)
)=1
T h e values of St’s a t solution are zero. The equilibrium relationships (Equations 2) are redefined to be:
r
7,
1
for K j
>
1.0 and
r n
for K ,
1
< 1.0.
Reasons for rewriting the equilibrium equations in a kinetic form are to avoid division by zero and to avoid nonlinearity of division. Since K , values may range between very large and very small numbers, two alternate definitions as in Equations 11 and 12 are used to help scaling. Scaling of the functions is important in the convergence of the l l a r q u a r d t algorithm. Each equation in Equations 10-12 is scaled by using a constant scale factor in the denominator. The scale factor found most effective is the sum of all absolute terms on the right side of the equation and
success. The constants in the derivatives are updated with each new base point. If a vapor phase is present, the above derivative calculation is done using vapor composition. Since the vapor phase nonideality is usually less significant’,the approximate derivatives will be more accurate. The convergence of the problem will also be faster. The reactions must be linearly independent for this approach to produce a solution. As explained by Kandirier and Brinkley (1950), this can be attained by ensuring that the rank of the stoichiometric coefficient matrix must equal the number of reactions. Algorithm
The algorithm proceeds as follows: 1. Initialize Fo,t
2.
1
K For example the scale factor for Equations 10 is
3. 4. 5. 6.
The values of the variables in the scale factor are those a t the base point. They are constants but are updated frequently. The partial derivative matrix, D, used by the blarquardt (1963) algorithm has elements as follows:
i
= Fin i
=
l,,..,n
=
+
+
The flash calculation is done using a nonlinear one-variable equation solving algorithm for a to satisfy Equation 7 . f
S = 1, . , . , m
i
j = l,...,m 0 If better estimates are availa,ble, they should be used. Calculate Si;i = 1 , . , . , m n, using the vapor-liquid flash calculation to get xi’s and yi’s. Calculate scale factors Calculate D Obtain new values for t j and Fout using Narquardt’s method Obtain new Si;i = 1 , . . . , m n using the vapor-liquid flash calculation to get xi’s and yi’s. If Si’s are not close to zero, return to step 3 and repeat. i j
+ l , m +2, ..., m + n
= m
azo--b F oass ut,t-m
as,
ayk
%zazEz
- k=l
a=O
S=
m
+ 1 . .. m + n
where
ass
~-
bFout, 1-rn
(13) -
as,
Direct iteration method is used for yz and q l for nonideal solutions. Example
or
To illustrate the new method, a problem was solved involving phase equilibrium and esterification reaction. The system contained ethanol, acetic acid, ethyl acetate, and water. The esterification reaction is CHaCHZOH Ethanol 61.i
Equations 14 arid 15 are valid if and only if the solution is ideal-i.e., the phase equilibrium constants are not functions of composition. These approximate analytical derivatives have been used, however, for nonideal systems with reasonable
+ CHaCOOH & CHsCOOCH2CHs + HzO -icetic acid
Ethyl acetate
I n any phase equilibrium calculation involving these chemical species, this reaction must, also be considered. Data required for the phase and chemical equilibrium include vapor pressures, liquid activities, and a n equilibrium constant for the reaction. For a n adiabatic calculation, enthalpy data would also be required. The vapor pressure and liquid activity data shown in Table I were used for t,his esample problen i. The reaction equilibrium constant was coniputed from Gibbs free energy data in Table 11. I n dealing with acetic acid, we know that dimers and trimers of acetic acid form in the vapor phase. The effect of this phenomenon was accounted for in the Wilso~i liquid activity coefficient model (Wilson, 1964) as used here. Ind. Eng. Chem. Process Des. Develop., Vol. 1 2 , No. 1 , 1973
83
Table 1. Vapor Pressure and liquid Activity Data for System EtOH-HAC-EtOAc-HzO (Suzuki et al., 1970)
Table 111. Esterification Equilibrium Stream names
Feed
Liquid
Vapor
0.0 50.0000 0.0 50.0000 -~ 100.00 0.0 358.15 101.325
10,2965 6,7876 5.0415 14,8625 ~. 36,99
16,9205 15.9954 22.1756 7.9204 63.01
0.0 2303.5 0.0 3002,6 5306,l
185.5 312.7 444.2 _ _892.5 . 1834.9
304.8 736.9 1963.7 475.6 3471.1
Mol y* 0.0 27.8373 Water Ethanol 50,0000 18,3508 0.0 13.6300 Ethyl acetate 50.0000 40.1819 Acetic acid Totals 100.00 100,OO
26.8529 25,3846 35.1927 12.5697 -~
Antoine constants Component
b
C
1440 52 1644 05 1238.71 1668 21 b
-60 44 -39 63 -56.15 -45 15
a
Ethanol Acetic acid Ethyl acetate Water
9 9 9 10 loglo P o
95614 6845 22298 09171 = a
-
___
T+c
Wilson’s parameters System
Azi
AI2
Acetic acid-ethanol Acetic acid-water Acetic acid-ethyl acetate Ethanol-water Ethyl acetate-ethanol Ethyl acetate-water
0 0 0 0 0 0
27558 26838 61790 15347 55046 12353
2 1 0 0 0 0
28180 22642 89277 92038 76670 14907
j= 1
Kg mol Water Ethanol Ethyl acetate Acetic acid Totals Fraction vapor, mol Temp, K Press, kX/m2 Kg Water Ethanol Ethyl acetate Acetic acid Totals
0.0 358.15 101.325
~
Table II. Reaction Equilibrium Constants for EtOH HAC EtOAc HzO (Stull et al., 1969) r AG,’ K
e
+
1.101 x 107
300 400 500 600
9.879 X IO6 8.960 X lo6 8.206 X lo6
+
Tvt %
Water Ethanol Ethyl acetate Acetic acid
82,62 19.51 8.63 5.18
The specific problem was to flash a stream containing initially
~
Totals
0.0 43.4 0.0 56.6 100.0
~
10.1 17.0 24.2 48.6 100.0
~
1.0000 358.15 101.325
100.00
~
8.8 21.2 56.3 13.7 100,o
residence time to attain equilibrium. The solution speed of this calculation method allows use of the method in this way. Nomenclature (SI Units)
Mol % Ethanol Ethyl acetate Acetic acid Water
50
0 50
0
The system temperature was 358’K. The pressure was 1.013 X lo5 N / m 2 (1 a t m ) . Table I11 gives results of the equilibrium calculations as produced by the Monsanto FLOWTRAX block. The calculated extent of reaction was 27.2 kg mol, and the output stream fraction vapor mas 63%. This solution was obtained in 10 iterations which required 0.3 see on the I B M 370/165 computer. Conclusions
A general calculation method has been developed to aid the process designer in translating thermodynamic data into usable process information. The method can be used with chemical equilibrium K-values defined from free energies or from laboratory measured equilibrium compositions. This flexibility plus the rigor of the method have proved useful in solving problems in which data have been provided in either or both of the above ways. This calculation method can be used by itself to calculate the equilibrium composition of a given mixture of chemical species. It can also be used as part of a continuous material and energy balance calculation. I n this instance, the equilibrium calculation would describe reactions with sufficient 84 Ind. Eng. Chem. Process Des. Develop., Vol. 12, No. 1, 1973
D
= matrix of partial derivatives (Jacobian) required by Marquardt’s algorithm Si = fugacity of species i, X / m 2 Fin i = initial amount of species i, kg mol Fout = final amount of species i (or estimate), kg mol Gio = standard Gibbs free energy of formation of species i, J AGjo = standard Gibbs free energy of reaction for reactionj, J ki = phase equilibrium ratio for species i ( = yi/’zi) K , = chemical equilibrium constant’ for reaction j (see Equation 2) m = number of reactions n = number of components p i = partial pressure of species i P = system pressure, Y / m 2 R = universal gas constant Si = error in equation i, described in Equation 10-12 T = temperature, K xi = mole fraction of species i in liquid phase y i = mole fraction of species i in vapor phase z i = mole fraction of species i in total system including liquid and vapor
GREEKLETTERS CY
yi pij
& qi
= = = = =
mole fraction vapor in the system activity coefficient of species i stoichiometric coefficient of species i in reaction j extent of reactionj, kg mol fugacity coefficient of species i ( = S / P )
SUBSCRIPTS
i ,j
= =
refers to species number refers to reaction number
t
=
refers t’o total moles
literature Cited
Clasen, R . J., “Th:,Sumerical Solut,ion of the Chemical Equilibrium Problem, Rand Corp. AIemorandum RL1-4345-PRj January 1965. Cruise, D. R., J . Phys. Chem., 68, 379i (1964). Goldwasser, S. R., Znd. Eng. Chem., 51, 595 (1959). Kandirier, H. J., Brinkley, S. R., ibid., 42, 850 (1950).
RIarquardt, D. W., J . SOC.Indus. A p p l . Math., 1 1 , 431 (1963). Smith, W. R., Jlissen, R. W., Can. J . Chene. Eng.,,?6, 269 (1968). Stull, D. R., Westrum, E. F., Jr., Sinke, G. C., The Chemical Thermodynamics of Organic Compounds,” pp 229, 423, 449, 431,’Wiley, New York, Y.Y., 1969. Suzuki, Isao, Hiromasa, Komatsu, JIitsuho, Hirata, J . Chem. Eng. J a p . , 3, 152 (1970). Villars, D. S., J. Phys. Chem., 63, 521 !1959). White, W. B., Johnson, S. AI., Dantzig, G. B., J . Chem. Phys., 28, i 5 l (19511). Wilson, G. AI., J . Amer. Chem. SOC.,86, 127 (1964). Zeleznik, F. J., Gordon, S.,Ind. Eng. Chern., 6 0 (6), 27 (1968). RECEIVED for review April 26, 1972 ACCEPTED J ~ l 26, y 1972
Rate of Diffusion-Controlled Reaction Between a Gas and a Porous Solid Sphere Reaction of SO, with CaCO, Robert 1. Pigfordl Cniversity of California, Berkeleg, Calij. 94720
Glenn Sliger The M. W . Kellogg Co., Houston, Tex. 770g6
The reaction rate i s assumed to b e governed b y two diffusional processes: the diffusion of SO2 through the pores, and the diffusion of SO2 through a layer of solid reaction product which forms on the active solid surface progressively as the reaction proceeds. One or the other predominates, depending on particle size, gas composition, and temperature. The theory accounts qualitatively for the observations of Borgwardt ( 1 970)and Coutant et a/. (1 970).
O n e of the major sources of atmospheric pollution in the United States is the coal-fired power plant, lyhich frequent,ly consumes sulfur-laden fuel and produces flue gas containing a few thousand parts per million SO,.I n recent’ years there has been increasing interest in chemical processes for the removal of sulfur compounds from the flue gas sbream. Many of t’hese involve treatment of the effluent gas stream from the plant with alkaline liquid reagents in large scrubbers before the gas is discharged. Such processes are capable of nearly complete recovery of gaseous sulfur but the equipment may be very costly and may be hard to add to many existing power plants. -1process which can be applied to esisting plants involves the injection of finely po1vdered limestone or Dolomite into the fire box of the boiler, where the solid particles decompose, forming calcium oxide. This react,s n-ith sulfur dioside arid with escess osggen, forming calcium sulfite and calcium sulfate. The reaction product must be collected with the fly ash in dust collection equipment. I n the evaluation of the limestone injection processes, i t is necessary to know t’he rates of the chemical reactions between the gas and the solid in order to determine the quantity of To whom correspondence should be addressed.
solid required, its particle size, the temperatures a t the injection point and downstream, and the time of esposure of the solid to the gas. Previous studies of the rate phenomena have led to somewhat confusing results. Borgwardt’ (1970) measured the rates of reaction using beds of solid particles of CaO through which the gas passed. His times of esposure were on the order of one to two minutes, much longer than those which are of pract’ical interest. He assumed that’ the heterogeneous reaction rate was governed b y the speed of a chemical reaction and eraluated a n apparent first-order rate constant, which he reported varied with the particle size and r i t h the elapsed time of t,he reaction. Such ambiguous results are cause for caution in the application of his data. Coutant et al. (1970) carried out a more estensive and realistic series of esperiments a t the Battelle Memorial Institute, in which particles of limestone (type 2061 , screened to average size of about 0.00895-cm diameter) \yere injected into a current of hot gas which flowed upward t’hrough a n insulated tube. By varying the height of the particle injection point, various times of esposure could be obtained, from a fraction of a second to about t’hree seconds. The particles were collected and analyzed to determine the estent to which Ind. Eng. Chern. Process Des. Develop., Vol. 12, No. 1 , 1973
85
