5474
Ind. Eng. Chem. Res. 1997, 36, 5474-5482
Phase Equilibrium Calculations for Chemically Reacting Systems Roumiana P. Stateva† and William A. Wakeham* Department of Chemical Engineering and Chemical Technology, Imperial College of Science, Technology and Medicine, London SW7 2BY, U.K.
An overall strategy for phase equilibrium calculations within chemically reacting multicomponent nonideal systems is presented. The strategy involves three stages and employs an indirect method to determine the phase configuration of the reacting system at equilibrium. The latter is based on an existing technique for the solution of the phase equilibrium problem extended to include chemically reacting systems. Independently, a new algorithm for chemical equilibrium calculations in any prescribed number of phases is advocated. A detailed analysis of its computational efficiency on the basis of multiplication operations required is provided. Results for three chemically reacting systems, two of which exhibit complex phase equilibria, are given and compared with the results of other authors. Introduction Process simulation is an important component of modern chemical engineering. Missing or inadequate physical properties of the fluid systems can undermine the accuracy of a model within the simulation, or even inhibit the performance of the simulation completely. Because many chemical process simulations include distillation, stripping, evaporation, and extraction, the separation of fluid mixtures is a fundamental aspect of chemical process design. To accomplish the development of such a process effectively, reliable thermodynamic models are required to describe phase, and possibly simultaneous chemical, equilibrium between these fluids. In particular, process simulations need to be able to reliably and efficiently predict whether or not a given phase will split into multiple phases, then to identify the phases that exist at equilibrium, and finally to calculate the distribution of components within these phases. This task is a key step in any phase equilibrium calculation and is especially important for high-pressure phase equilibria where even a simple binary solution can exhibit complex phase behavior. This problem becomes even more acute and difficult when some or all of the components may undergo chemical reactions simultaneously with phase equilibrium. For such circumstances the global minimum of the Gibbs free energy for the system always yields the true equilibrium state. However, the sheer complexity of minimizing the Gibbs energy function for an unknown number of phases makes the process fraught with numerical difficulties. As a result there has been an incentive to seek alternative ways of solving the same general task. In this paper we present a new approach to the solution of problems of this kind that involve fluid phases. Specifically we further extend the application of an existing technique for the solution of the phase equilibrium to include chemically reacting systems. Independently, we provide a new algorithm for the solution of a chemical equilibrium problem in any prescribed number of phases. The judicious application of these two procedures within an overall strategy will be shown to provide an * Author to whom correspondence should be addressed. Telephone: +44 171 594 5005. Fax: +44 171 594 8802. E-mail:
[email protected]. † Permanent address: Institute of Chemical Engineering, Bulgarian Academy of Sciences, acad. G. Bontchev str., bl. 103, Sofia 1113, Bulgaria. S0888-5885(97)00264-9 CCC: $14.00
effective and robust means of solving the chemical and phase equilibrium (CPE) problem. Formulation of the Problem and Prior Work Many recent papers have addressed the problem of the evaluation of the simultaneous CPE calculations, among which should be mentioned the recent series of publications by McDonald and Floudas (1995b, 1997) and the thorough review article by Seider and Widagdo (1996). At constant temperature and pressure the condition of equilibrium in any thermodynamic system is that the Gibbs energy function for the system attains its global minimum:
min G(n) subject to An - b ) 0 0 e n e nT
}
(1)
where G(n) is the total Gibbs free energy of a system: NP N
G(n) ≡
∑ ∑ j)1 i)1
NP N
nijG h ij )
(
nij G0ij + RT ln ∑ ∑ j)1 i)1
fij f 0ij
)
(2)
and nij and G h ij are the mole numbers and partial molar Gibbs free energy (chemical potential) of species i in phase j, N and NP are the numbers of species and phases, fij is the fugacity of species i in phase j, Gij0 is the Gibbs free energy of formation of species i in the standard state for phase j, and f ij0 is the associated standard-state fugacity. There are essentially two approaches that can be taken in order to demonstrate whether an equilibrium solution corresponds with a global minimum of the Gibbs free energy (McDonald and Floudas, 1997). The separate optimization problems that are used for these two approaches are (1) the minimization of the Gibbs free energy for the system, and (2) the minimization of the tangent plane distance function to verify that an equilibrium solution corresponds to the global minimum. A number of computational algorithms for the calculation of phase, and CPE, which use either the first or the second approach have been proposed before. Their merits and disadvantages have been the subject of a © 1997 American Chemical Society
Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5475
thorough discussion (McDonald and Floudas, 1995a,b, 1997; Seider and Widagdo, 1996), and the analysis demonstrates that methods employing thermodynamic stability analysis can be very effectively used to improve the chances of finding the true solution. These methods have been found to be more reliable and robust than minimizing the Gibbs function directly, since theoretically to guarantee a global minimum for G, using global optimization, the formulation can contain a large number of variables, making a global optimization algorithm computationally expensive, particularly for larger systems (McDonald and Floudas, 1997). The present work makes a fresh attempt at the solution of the CPE problem and advocates an overall strategy that overcomes some of the problems of the existing algorithms. In its concept it is similar to the method of Gautam and Seider (1979) and the approach of Castier et al. (1989) in that it uses the tangent plane criterion to test for stability. The key point is that the stability analysis is used only once and its judicious application provides, at an early stage of the computations, information for the phase configuration of the chemically reacting system at equilibrium. The paper also introduces a new algorithm for the solution of a chemical equilibrium problem in any prescribed number of phases. The overall strategy involves the following three stages: (1) In the first stage a candidate equilibrium solution is obtained, applying an efficient and easy to implement algorithm for simultaneous CPE calculations in a homogeneous phase. (2) In the second stage, the stability of the candidate solution is checked, applying a rigorous analysis, on the basis of the minimization of a modified tangent-plane distance function (Stateva and Tsvetkov, 1994). The method attempts to locate all minima, but since global optimization methods are not applied, there is no theoretically-based guarantee that all solutions are obtained. However, the approach has a built-in element of self-recovery so that it is, to some extent, robust against the unlikely event of missing a minimum in the stability test. An automated procedure for determining a set of initial estimates is employed (Stateva and Tsvetkov, 1991). If the candidate solution is unstable, the phase configuration of the reacting system at equilibrium is determined by performing a simple sequence of two-phase flash calculation upon the candidate solution. (3) At stage three an algorithm for CPE calculations is applied upon the candidate solution. The identity and number of the phases, and the initial estimates for the component distribution within the phases, are determined at stage two. The Chemical and Phase Equilibrium Calculations Smith (1980) classifies chemical equilibrium formulations to be either stoichiometric or nonstoichiometric. Stoichiometric formulations require a knowledge of the stoichiometric coefficients of a linearly independent set of reactions and utilize the extents of the reactions as independent variables. On the other hand, in the nonstoichiometric formulations, elemental abundance constraints are used in the calculations. In the classical approach, for a chemical equilibrium in a single phase, M independent chemical reactions are defined and the mole numbers are replaced by
M
Fout ) Fin i i +
νijξj ∑ j)1
i ) 1, 2, ..., N j ) 1, 2, ..., M (3)
It should be noted that the stoichiometric coefficient νij is positive if the component i is a product of reaction j, negative if the component i is a reactant of reaction j, and zero if the component i is inert. With ξj as the independent variables the mass-action equations result in N
R
K j)
∏ i)1
() fi
νij
j ) 1, 2, ..., M
f 0i
(4)
In the case of, for example, a two-phase vapor-liquid system involving N components and M reactions, the equilibrium constant KR for reaction j is defined by N
∏ i)1
N
(xiPφLi )νij )
(yiPφVi )ν ∏ i)1
ij
) KRj j ) 1, 2, ..., M
(5)
The constant KRj may be determined experimentally or calculated from thermochemical data. The equilibrium compositions and the number of phases of the system can be obtained by solving a set of mass balance, phase, and chemical equilibrium equations for specified components and reactions. This is called the K-value method, because distribution coefficients KRj (by analogy with the phase equilibrium) are introduced. Numerous alternatives to the classic “Kvalue” approaches have been introduced. Seider et al. (1991) and Seider and Widagdo (1996) provide a review of the literature on the existing solution methods. One generalized algorithm for solving CPE in twophase systems providing considerable improvements of the alternative K-value approach has been proposed by Xiao and co-workers (1989). Using effective K-values based on single-phase compositions, and a rearrangement of the chemical reaction and phase equilibrium calculation loops proposed originally by Sanderson and Chien (1973), they reported accelerated convergence of these methods when applied to an esterification reaction and several electrolytic solution equilibria problems. In this work an efficient modification of the method of Xiao and co-workers is advocated. To set the stage, the mathematical background of the algorithm of Xiao et. al. (1989) will be presented briefly. Then the new modification will be introduced and an analysis of its computational efficiency will be carried out. The KZ Algorithm of Xiao et al. (1989). The phase equilibrium and the flash equations for a two-phase vapor-liquid (VL) system, are
zi )
Fout i
i ) 1, 2, ..., N
N
(6)
Fout ∑ i i)1 zi ) Ryi + (1 - R)xi yi ) Kixi N
S(R) ≡
i ) 1, 2, ..., N i ) 1, 2, ..., N
(Ki - 1)zi
∑ i)11 + R(K
i
- 1)
)0
(7) (8) (9)
5476 Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997
Xiao et al. argued that it is efficient for an algorithm to reduce the number of K-value evaluations as much as possible and suggested moving the entire flash loop to the outside and keeping the inner loop for the chemical equilibrium computations. They rearranged eqs 3 and 5 and 7-9 to yield
Introducing eq 3 into eq 15, the following system of M nonlinear equations is obtained: N
Φj(ξ1, ξ2, ..., ξM) ) N
ln(
N
(zi) ∏ i)1
νij
) Kzj
where N
Kzj ) KRj
j ) 1, 2, ..., M
(10)
(
)
Pφi
∏ i)1 1 + R(K
i
N
-νij
(11)
- 1)
N
(zi)-ν )ν 0 ∏ i)1 i)1 ij
ij
ij
ij
j ) 1, 2, ..., M (12a) for Kzj g 1 and N
∏ i)1
Sj ≡ (
(zi)-νij)νij0 ∏ i)1
(
ij
ij
j ) 1, 2, ..., M (12b) for Kzj < 1. The material balance equation, eq 3, was redefined to be M
SM+i ≡ Fout ) Fin i i +
Fout ∑ i (ξ1, ξ2, ..., ξM))) - ln Kz ) 0 i)1 j
j ) 1, 2, ..., M (16)
In this formulation, effects of multiple phases and their nonideality on the reaction equilibria have been merged into a set of new parameters Kzj, j ) 1, 2, ..., M, which represent “effective” chemical equilibrium constants. Xiao et al. (1989) called their new algorithm the KZ algorithm and, in order to avoid division by zero, they rewrote eq 10 into a kinetic form:
Sj ≡ Kzj(
νij(ln Fout ∑ i (ξ1, ξ2, ..., ξM) i)1
νijξj ∑ j)1
i ) 1, 2, ..., N
(13)
The system of eqs 12 and 13 was solved for the extents and Fout (total number of the iterative variables was M i + N) by a modified Marquardt method (Xiao et al., 1989). The New Modification of the Chemical Equilibrium Computations. The present paper advocates a modification of the KZ algorithm specifically for the part concerning the chemical equilibrium calculations. As a result of this modification the initial system, eqs 12 and 13 with M + N nonlinear equations, is partitioned into two independent parts, namely a system of M nonlinear equations and N linear expressions. They are solved subsequently. The partitioning is realized by introducing the following transformations: After taking the logarithm, eq 10 can be rewritten in a different form:
The unknown variables of the system are ξj, j ) 1, 2, ..., M, while the parameters Kzj can be calculated directly from eq 11. The linear functions, eq 3, constitute the second independent part. Thus Fout i , i ) 1, 2, ..., N can be calculated directly when ξj, j ) 1, 2, ..., M are known. The Algorithm. The new modification of the chemical equilibrium computations is introduced into an algorithm for simultaneous CPE calculations. It follows the “inside-out” organization originally advocated by Boston and Britt (1978) for flash calculations, and further improved by Xiao et al. (1989). That is, the phase equilibrium calculations with constant K-values are performed in an outer loop and an inner loop is used for chemical equilibrium computations. To ensure objectivity of the comparison of the effectiveness and efficiency of the two methods, all further developments and discussions in this section will be kept within the framework of the paper of Xiao et al. (1989). If applied to CPE calculations only, the present algorithm requires input information, and initial estimates. The input includes the following parameters: temperature, pressure, number of components, number and identity of the equilibrium phases, number of independent chemical reactions and their chemical equilibrium constants, and stoichiometric coefficients. Here we present, as an example, a two-phase CPE calculation in order to be consistent with the original discussion of the KZ algorithm, because it allows a meaningful evaluation of the improvements addressed. Its main steps can be briefly summarized as follows: 1. Input T, P, zi, N, NP, M, KRj , and νij (i ) 1, 2, ..., N; j ) 1, 2, ..., M). 2. Estimate Ki, i ) 1, 2, ..., N. (The next section of the paper explains how to obtain initial estimates for the Ki.) 3. Estimate R. The phase fraction is estimated automatically by a secant method within the following boundaries:
Rmin )
1 1 e R e Rmax ) 1 - Kmax 1 - Kmin
as suggested by Whitson and Michelsen (1989). 4. Calculate xi,yi and the fugacity coefficients φi for the components in the liquid and vapor phases, applying the corresponding thermodynamic models. 5. Find Kzj, j ) 1, 2, ..., M from eq 11. 6. The Jacobian of system (16) is known analytically:
( ) N
N
ln Kzj -
νij ln zi ) 0 ∑ i)1
j ) 1, 2, ..., M (14)
∂Φj
N
) ∂ξl
Substituting eq 6 into eq 9 yields
νij ∑ i)1
∑ νkl
νil
k)1
-
Fout i
N
∑
l ) 1, 2, ..., M
(17)
Fout k
k)1
N
ln Kzj -
∑ i)1
N
Fout ∑ i )) ) 0 i)1
νij(ln Fout - ln( i
j ) 1, 2, ..., M (15)
and we use a Newton’s method to solve eq 16 for ξi, j ) 1, 2, ..., M. Upon convergence new values for ξ1, ξ2, ..., ξM are obtained.
Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5477
7. Calculate Fout from eq 3. i 8. Calculate new values for znew from eq 6. i 9. Check whether eq 9 is satisfied. If satisfied, continue with step 10; otherwise return to step 4. 10. Calculate new values for Ki and check the following criterion: N
∑ i)1
(
)
Knew - Kold i i Knew i
2
e 10-12
i ) 1, 2, ..., N
If the criterion is satisfied, end of the algorithm; otherwise return to step 3. The algorithm has an analogous structure in the case of three-phase CPE calculations, but is not given here owing to a lack of space. In order to further simplify the presentation we denote the two-phase algorithm by Algorithm II and the three-phase by Algorithm III. If Algorithm II is applied for simultaneous CPE calculations in a homogeneous phase, then only steps 5-8 will be required. These chemical equilibrium calculations within a single phase we denote by Algorithm I. We will comment briefly on the efficiency of the present algorithm. The system Φ (eq 16) to be solved in the chemical equilibrium calculations consists of M nonlinear equations as compared to the system of eqs 12 and 13, solved by Xiao et al. (1989), which contains M + N nonlinear equations. The corresponding Jacobian in the former case is of dimension (M × M) while in the latter case (M + N) × (M + N). This constitutes already a considerable potential savings in computational effort. It is obvious that the set of nonlinear equations solved in our case, eq 16, differs in form from the set of nonlinear equations solved by Xiao et al. Equation 16 seems superficially to be more cumbersome than the equivalent system of Xiao et al. (1989). Also it might seem that they could exhibit more complex solution space than the original set (e.g., changing the number of its solutions) and lead to an increase in the computational efforts and therefore computational time to find the true solution. However, the solution space will not be changed since the transformation of eq 10 to eq 16, namely taking the logarithm and then introducing two linear functions, will not lead to a change in the number of its solutions, since the logarithm and the linear functions are monotonic. With regard to the numerical effort and the computational time required to find the solution of eqs 16 and eqs 12 and 13, respectively, we have made a comparison of the KZ algorithm and Algorithm II. This was done on the basis of the number of multiplication operations required and carried out by a computer on one iteration step of the chemical equilibrium computations and from just one set of initial estimates. The comparison demonstrated that the number of operations required by the present algorithm is usually much less than the number of operations required by the KZ algorithm, the latter being especially pronounced when N is big (the number of components is large). This result quantified our observation that one of the main advantages of our method over the method of Xiao et al. arises from the fact that the dimension of the system to be solved is considerably decreased. Finally, an additional advantage of the present modification arises from the fact that we use the Newton method to solve system Φ. When the conditions for the
convergence of the Newton method hold then it converges at a second order (Kantorovich and Akilov, 1984). The Overall Strategy First Stage. The first stage of the overall strategy requires the application of Algorithm I upon the initial system. The output is a homogeneous liquid or vapor phase with the distribution of the components denoted as the candidate solution z*, and the extent of the reaction ξ1. Second Stage. The second stage consists of two substages. Firstly a rigorous stability analysis is run upon the candidate solution and if it is proved to be unstable, a sequence of two-phase liquid-liquid and/or liquid-vapor flashes follow to determine the number and identity of the phases at equilibrium.. This “indirect” method is based on the assumption that the components in the candidate solution will not further react appreciably among themselves, because the system is near equilibrium. The correctness of this assumption and the certainty of obtaining the true phase configuration could be tested ex post facto. Furthermore, it is backed by the experimental results of Sherman et al. (1996), and others, who demonstrate that the phase behavior of reacting systems can be estimated from the binary information on its nonreacting pairs, and is presented numerically in the Test Examples section of the paper. The rigorous method for thermodynamic stability analysis, performed once and only once upon the candidate solution, is based on the tangent-plane criterion (Baker et al., 1982; Michelsen, 1982), but uses a different objective function (Stateva and Tsvetkov, 1994). For an N-component mixture, the essential point of the technique is to locate all the zeros (y*) of a functional Φ(y): N
Φ(y) )
[ki+1(y) - ki(y)]2 ∑ i)1
(18)
where
ki(y) ) ln φi(y) + ln yi - hi
i ) 1, 2, ..., N (19a)
hi ) ln zi + ln φi(z)
i ) 1, 2, ..., N (19b)
and kN+1(y) is equated to k1(y). From eq 18 it follows directly that when y ) y*, where k1(y*) ) k2(y*) ) ... ) kN(y*), then the functional Φ(y*) ) 0, which is its minimum value. The zeros of Φ(y), y*, correspond with points on the Gibbs energy surface, where the local tangent hyperplane at each y* is parallel to that at z. To each vector y*, a quantity k* corresponds, such that
k* ) k* for any i (y*) ) ln φi(y*) + ln y* i - hi i ) 1, 2, ..., N (20) Furthermore, k* geometrically represents the distance in the Gibbs energy hypersurface between the hyperplane at z and the hyperplane at any y*. This distance can be either positive or negative. When all calculated k* are positive, the chosen phase configuration is stable. The essential step of the technique is the determination of all the zeros of Φ(y). The functional Φ(y), which is such that its zeros are its minima, as well as the fact that it is easily differentiated, means that it is more effective to employ a nonlinear minimization technique to determine its stationary points. This is because
5478 Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997
multidimensional minimization techniques are more robust than root-finding procedures and because their convergence is superlinear (Fletcher, 1987). We use here the BFGS method. The implementation of the procedure requires the identification of multiple initial estimates of the stationary points. An automated procedure for the generation of these initial estimates has been proposed by Stateva and Tsvetkov (1991) which begins from an original set founded upon the almost pure components. If one or more of this original set of estimates should fail to lead to a new stationary point, another set is generated by means of a line search technique (Fletcher, 1987) from which the Broyden, Fletcher, Goldfarb, and Shauno method can be restarted. In this way a method has been created which finds almost all of the zeros of the functional Φ(y), but there is no theoretically-based guarantee that it will always find them all. As mentioned earlier, the method leads not only to the zeros of Φ(y) (each of which corresponds with a possible phase) but also to the corresponding values k* for each y*. The signs of the k* determine the stability of solution tested. If it is not stable, then the density of each potential phase is compared with a density intermediate between the liquid and vapor. These comparisons serve to characterize each potential phase as either liquid-like or vapor-like and normally result in only one vapor-like phase. The stationary points located on the first substage are used further on the second substage as starting points (initial estimates) for a search for an equilibrium solution which will have a lower Gibbs energy. To do that, all possible y* pairs are taken together since there is no way of telling which minima of the functional correspond to the global solution. Thus, the chance of using guesses that may in the end lead to local optima or even infeasible solutions is eliminated. It is possible that some guesses may lead to a phase configuration physically inconsistent with the overall composition specified for the mixture and they are rejected. For the remaining pairs, two-phase flash calculations (in any order) are performed by any available technique; for example the “negative-flash” method. These calculations lead to values of the Gibbs energy and that with the lowest value is selected for further processing, which determines whether the inclusion of a third phase leads to further stability. If not, then the two-phase result is the final, stable solution. The three-phase system can be subjected to a further test of the same kind if necessary. It should be mentioned that the various twophase flash calculations can be performed in parallel, if it is desired. The following should be highlighted as well. The commonly accepted methods for phase equilibrium calculations assume that if a postulated solution is found to be unstable, then a phase must be added or deleted. However, the true equilibrium solution may in fact be a solution with the same number of phases but either with a different identification (VL vs LL) or with a different distribution of the components within those phases (e.g., in a four-component system there might be several three-phase solutions found, of which only one corresponds to the global minimum). The present approach overcomes this by selecting on each step the solution with the minimum Gibbs energy for further processing. This simple, easy to implement, and self-recovering procedure provides excellent initial estimates for the subsequent two/three CPE calculations, thus improving, almost perfecting, their convergence characteristics.
In summary, stage two of the overall strategy performs firstly a rigorous thermodynamic stability analysis to test the stability of the candidate solution z*. If z* is stable then the chemical reaction takes place in a homogeneous either liquid or vapor phase. If the candidate solution is unstable then z* is tested whether it is stable toward a VLL split (for example), applying a sequence of two-phase flash calculations only. Once the phase configuration at equilibrium is identified, the calculations proceed with stage three. Third Stage. The distribution of the components among the equilibrium phases and the extent of the independent reactions are determined in this stage using Algorithm II and Algorithm III. The identity and number of the phases, and the initial estimates, required by either Algorithm II or Algorithm III applied in stage three are known from stage two. Summary of the Overall Strategy. The generalized scheme of the overall strategy to the CPE calculations is summarized as follows: 1. Given z, T, P, and information for the independent chemical reactions. 2. Run Algorithm I. The output is the candidate solution z* and ξh1. 3. Run a stability analysis upon the candidate solution z*. If the candidate solution is stable, the reaction takes place in a homogeneous phase. END. Otherwise, continue with step 4. 4. Run two-phase flash calculations upon z*, using as initial estimates the stationary points obtained in step 3. Choose among the VL and LL solutions the one with the minimum G. Let it be a LL type. 5. Test this particular solution toward a VLL split by performing a vapor-liquid flash on its lighter phase. If the phase split upon convergence is outside the physically acceptable bounds, the equilibrium phase configuration is liquid-liquid. Continue with step 6. Otherwise, the equilibrium phase configuration is VLL. Intial estimates for Algorithm III are found by performing an additional LL flash on the liquid phase of the above VL solution. Continue with step 7. 6. Run Algorithm II upon z*. Use as initial estimates the converged mole fractions of the LL solution determined in step 4. The application of Algorithm II yields a solution, characterized with ξh2 and distribution of the components within the phases which corresponds with the equilibrium one. The final extent of the reaction is ξhfinal ) ξh1 + ξh2. (See Appendix A.) END. 7. Run Algorithm III upon z*. Use as initial estimates the converged mole fractions of the solutions determined in step 5. The application of Algorithm III yields a solution which is characterized with ξj3 and distribution of the components within the phases which corresponds with the equilibrium one. The final extent of the reaction is ξhfinal ) ξh1 + ξh3. (See Appendix A.) END. Test Examples To test the robustness of the new strategy for the solution of the CPE problem, examples of chemical reactions with complex phase behavior must be examined. There are, however, only a few examples available in the literature, and two of them are presented here. We apply in all our calculations the MHV-2 model (Dahl and Michelsen, 1990), and the modified UNIFAC model (Larsen et al., 1987; Dahl et al., 1991), for the excess Gibbs energy in connection with the mixing rule. The robustness (in terms of locating the correct solution) of the thermodynamic stability analysis and the phase
Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5479 Table 1. Equilibrium Phase Compositions (Mole Fractions) for the Esterification Reaction (Eq 21) at T ) 355 K and P ) 0.1013 MPaa present study
Table 2. Equilibrium Liquid Phase Compositions (Mole Fractions) for the Methanol Synthesis (Eq 22, 22a) at T ) 473.15 K and P ) 30 MPa
Xiao et al. (1989)
component
feed
liquid
vapor
liquid
vapor
EtOH HAC EtAC H2O
0.5 0.5
0.0675 0.2243 0.1545 0.5537
0.1029 0.0554 0.4813 0.3604
0.0686 0.2376 0.1373 0.5565
0.0862 0.0624 0.4551 0.3963
a Calculated values of the phase split and extent of reaction: This study: R ) 0.767; ξfinal ) 0.4243. Xiao et al. (1989): R ) 0.877; ξfinal ) 0.8321.
identification procedures, and their reliable application to nonreacting highly nonideal systems, has been demonstrated extensively in a number of publications (e.g., Stateva and Tsvetkov, 1994; Stateva, 1995), so no details of stage two will be given here. Furthermore, it must be stressed that the exactitude of the thermodynamic model is itself of no significance here. Results for three test problems involving different reactive systems, with different number of components and reactions, are presented. All the systems have been studied by other researchers and can be regarded as a suitable test bed for comparisons. Esterification of Acetic Acid with Ethanol. An equimolar mixture of ethanol and acetic acid reacts reversibly according to the following esterification reaction:
C2H5OH + CH3COOH S CH3COOC2H5 + H2O (21) EtOH HAc EtAc The esterification reaction does not offer any significant challenges but has been extensively studied and used by many authors as a bench-mark example (see for example Sanderson and Chien, 1973; Castillo and Grossman, 1981; Castier et al., 1989; Xiao et al., 1989; McDonald and Floudas, 1997; Perez Cisneros et al., 1997). Almost all known methods for CPE calculations in this system predict two-phase vapor-liquid equilibrium. In the present study the reaction equilibrium constant is calculated using thermochemical data for all compounds from Stull et al. (1969). The effects of dimerization of acetic acid in the vapor phase has not been taken into account. The results obtained at T ) 355 K and P ) 0.1013 MPa are presented in Table 1 and compared with the results of Xiao et al. (1989). The latter use the UNIQUAC model for the liquid phase activity coefficients, with its parameters generated by the original UNIFAC model, and assume the equilibrium vapor phase to be ideal. Methanol Synthesis. The methanol synthesis reaction is studied at T ) 473.15 K and two values of pressure, P ) 10.13 and 30 MPa. The reactions involved are:
CO + 2H2 S CH3OH
(22)
CO2 + H2 S CO + H2O
(22a)
As suggested by Castier et al. (1989) and Gupta et al. (1991), the influence of the presence of an inert oil in the initial system on the phase equilibrium of the reaction is analyzed. The inert oil is modeled as n-octadecane C18H38. The values for the equilibrium constants of reaction 22 and reaction 22a are calculated from the expressions
component
feed
CO H2 CO2 CH3OH H2O CH4
0.15 0.74 0.08 0.03
liquid (this study)
liquid (Castier et al.)
liquid (Gupta et al.)
traces 0.0948 traces 0.6371 0.2488 0.0193
traces 0.0962 0.0002 0.6354 0.2436 0.0246
traces 0.0402 0.0001 0.6709 0.2732 0.0156
Table 3. Equilibrium Vapor Phase Compositions (Mole Fractions) for the Methanol Synthesis (Eq 22, 22a) at T ) 473.15 K and P ) 30 MPaa component
feed
CO H2 CO2 CH3OH H2O CH4
0.15 0.74 0.08 0.03
vapor (this study)
vapor (Castier et al.)
vapor (Gupta et al.)
1.33 × 10-5 0.6493 traces 0.2120 0.0464 0.0923
6 × 10-5 0.6589 0.0005 0.2053 0.0473 0.0875
0.0001 0.6693 0.0004 0.2040 0.0347 0.0915
a
Calculated values of the phase split and extents of reactions: final final This study: R ) 0.4968; ξ(22) ) 0.23; ξ(22a) ) 0.08. Castier et al. (1989): R ) 0.5112; extents of reactions not reported. Gupta et final final al. (1991): R ) 0.5258; ξ(22) ) 0.2298; ξ(22a) ) 0.0799. Table 4. Equilibrium Phase Compositions (Mole Fractions) for the Methanol Synthesis in the Presence of a Heavy Hydrocarbon C18H38 at T ) 473.15 K and P ) 10.13 MPaa component
feed
CO H2 CO2 CH3OH H2O CH4 C18H38
0.1071 0.5286 0.0571
a
0.2143 0.0214 0.0715
water-rich phase
hydrocarbonrich phase
vapor phase
1.27 × 10-10 0.0071 2.96 × 10-6 0.2870 0.7047 0.0011 2.7 × 10-6
4.8 × 10-9 0.06 3.18 × 10-12 0.1418 0.007 0.0210 0.7702
5.63 × 10-8 0.5328 7.27 × 10-12 0.2274 0.1635 0.0752 0.0010
final final ψ1 ) 0.4843; ψ2 ) 0.3780; ξ(22) ) 0.1642; ξ(22a) ) 0.057.
suggested by Bissett (1977) and Cherednichenko (1953), which reveal KRj to be a function of temperature. At T ) 473.15 K and P ) 30 MPa, and in the absence of an inert oil, a two-phase VL solution is found to be the equilibrium one. The results are compared with the phase prediction and calculations of Castier et al. (1989) and Gupta et al. (1991) and are presented in Tables 2 and 3. The former authors use the Soave-RedlichKwong cubic EOS to model the equilibrium phases, while the latter use the Trebble-Bishnoi EOS (1988) and calculate the chemical equilibrium constants for the reactions using the data from Reid et al. (1977). At T ) 473.15 K and P ) 10.13 MPa, and with the inert oil present, a three-phase vapor-liquid-liquid solution is obtained. The results of this study are in a good agreement with those reported by the other two groups of authors who also predict a three-phase VLLE (see corresponding references), but the latter are not shown owing to a lack of space. Table 4 gives the equilibrium three-phase results, obtained with Algorithm III, which was initialized with the estimates shown in Table 5. The computation demonstrates stable convergence characteristics, and the desired solution is found in four iterations. Propene Hydration. The hydration of propene to 2-propanol is presented at T ) 353.15 K and P ) 0.1013 MPa and in the presence of a scarcely miscible inert component n-nonane to study how the latter influences the phase equilibrium (Castier et al., 1989). The value
5480 Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 Table 5. Initial Estimates (Mole Fractions) for Algorithm III for the Methanol Synthesis Reaction component
water-rich phase
hydrocarbon rich phase
vapor phase
CO H2 CO2 CH3OH H2O CH4 C18H38
1.72 × 10-11 0.0038 7.45 × 10-12 0.2286 0.7671 0.0005 8.13 × 10-8
1.2 × 10-9 0.0480 1.02 × 10-10 0.1269 0.0131 0.0139 0.7981
1.18 × 10-8 0.3711 2.19 × 10-10 0.2286 0.3517 0.0451 0.0035
propene water 2-propanol n-nonane a
R ) 0.628;
liquid
vapor
0.313 0.626
9.3 × 0.9978 0.0021 5.5 × 10-6
0.3445 0.4035 0.1545 0.0975
0.061 ξfinal
) 0.0973.
Table 7. Equilibrium Compositions (Mole Fractions) for the Propene Hydration (Eq 23) at T ) 355 K and P ) 0.1013 MPa (Case b)a component
z
propene water 2-propanol n-nonane
0.25 0.5
a
0.25
aqueous phase
organic phase
vapor phase
1.29 × 10-4 0.9998 8.15 × 10-6 4.85 × 10-6
0.0200 0.0007 1.49 × 10-5 0.9793
0.4980 0.4042 6.30 × 10-4 0.0972
ψ1 ) 0.3004; ψ2 ) 0.2064; ξfinal ) 3.16 × 10-4.
for the equilibrium constant of reaction 23 is taken from Vidal (1974). CH3 CH propene
CH2 + H2O water
CH3 CH CH3 OH 2-propanol
aqueous phase
organic phase
vapor phase
propene water 2-propanol n-nonane
8.26 × 10-5 0.9970 0.0029 5.73 × 10-6
0.0121 0.0008 0.0052 0.9818
0.3063 0.4035 0.1551 0.1351
Table 9. Equilibrium Phase Compositions (Mole Fractions) for the Propene Hydration (Eq 23) at T ) 353.15 K and P ) 0.1013 MPa (Case c)a propene water 2-propanol n-nonane
z
10-5
component
component
Table 6. Equilibrium Phase Compositions (Mole Fractions) for the Propene Hydration (Eq 23) at T ) 353.15 K and P ) 0.1013 MPa (Case a)a component
Table 8. Initial Estimates (Mole Fractions) for Algorithm III
(23)
The results refer to a mixture whose initial mole numbers are 10 mol of propene and 20 mol of water. Applying the new strategy to simultaneous CPE we have determined that (a) amounts of n-nonane below 1.96 mol in the initial system yield a vapor-liquid equilibrium; (b) amounts of n-nonane between 1.96 and 381.2 mol in the initial system yield a vapor-liquidliquid system; (c) amounts of n-nonane above 381.2 mol in the initial system yield a liquid-liquid system. The above results are in a good agreement with the results of Castier et al. (1989), who predict VLE below 2.072 mol of n-nonane, LLE above 378 mol of n-nonane, and a VLLE for n-nonane between these two quantities. Our calculations slightly overpredict the VLL region in comparison with that of Castier et al. (1989) who model the vapor phase as an ideal vapor, employ the modified UNIFAC (Larsen et al., 1987) for the liquid phases, and check the stability of each new phase configuration after deleting or adding phases. The calculations for case (a) are demonstrated for 1.95 mol of n-nonane in the system; the mixture composition is at the boundary between the VL and VLL regions. Table 6 gives the equilibrium two-phase results, obtained with Algorithm II. The calculations for case (b) are performed for 10 mol of n-nonane. Table 7 shows the equilibrium three-phase results, obtained with Algorithm III, which was initialized with the estimates shown in Table 8. The path to the solution for case (c) proved to be more complex than for cases (a) and (b) for Castier et al. (1989), since it required the addition and subsequent removal of a vapor phase. In our case the
a
R ) 0.04;
z
aqueous phase
organic phase
0.0233 0.0466
1.3 × 0.9970 0.0030 5.7 × 10-6
0.0191 0.0008 0.0052 0.9749
0.9301 ξfinal
10-4
) 5.1 × 10-3.
true LL configuration at equilibrium was easily identified once the LL phase split, which was with the lowest Gibbs energy among the three two-phase splits performed upon the candidate solution z*, proved to be stable to further splits. Table 9 gives the equilibrium two-phase results, obtained with Algorithm II for 400 mol of n-nonane in the system. Discussion The three stages of the overall strategy presented here realize a judicious combination of a stability analysis, phase identification routines, and calculations of simultaneous CPE. The strategy provides an efficient and robust means for the solution of the CPE problem and gives the following advantages: (1) In stage one Algorithm I is used. It applies the Newton method and usually is extremely fast. Though there is currently no practical evidence that the Newton method is likely to fail because the conditions for its convergence are not satisfied, it is only fair to say that we have not been stimulated to look for counter examples where it will not work. (2) In stage two the hybrid stability routine leads to a prediction of the number and identity of the phases at equilibrium at an early stage of the calculations (see example Propene Hydration, case (c). (3) In stage three of the overall strategy either Algorithm II or Algorithm III is employed. Both algorithms demonstrate excellent convergence characteristics firstly because they exploit the set of initial estimates, provided by the hybrid stability routine of stage two, and secondly because they involve the Newton method in the chemical equilibrium calculations loop. Thus, even for the cases of a multicomponent system with more than one chemical reaction (e.g., the methanol synthesis; two reactions, and seven components, of which six are reacting), the three-phase CPE calculations converged in only four iterations. (4) Finally, the present strategy is not limited by the type of thermodynamic model applied. It works both for symmetric and asymmetric models. Among the shortcomings of the overall strategy the following should be mentioned. In stage two a rigorous thermodynamic stability analysis is applied (as the first step of the hybrid routine) to test the candidate solution z* for stability. Just as for any other technique, which does not use a global optimization method, its robustness (in the sense of finding all stationary points) cannot be proved mathematically. Although a theoretically-based guarantee of finding all stationary points cannot be given,
Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5481
the technique does provide an automatic means of finding almost all of them. Furthermore, if the initial stability test fails for some reason to find all zeros of the functional, then the process of the pairing of solutions may in very many cases provide a means of finding this same phase configuration by a different route. For example, among the pairs of phases flashed in the second part of the technique, one may fail to find a solution with a positive phase split and the true, minimum Gibbs energy simply because of the absence of one initial stationary point. In such a case, the subsequent process will use the pair of phases with the current lowest Gibbs energy values for subsequent processing, and the process will then recover the missing phase in one or more steps. It follows that the progression of the procedures in stage two has an element of self-recovery built into it. There is, however, a price to pay for the robustness of stage two, because the stability test is time consuming. The second step of the hybrid routine involves two-phase flash calculations which practically converge in less than three iterations and can be performed in parallel. We have not given a theoretical proof of the “indirect method” of stage two. Thus, the correctness of the conclusion about the number and identity of the phases at equilibrium is tested ex post facto. We tried the new strategy extensively on a large body of examples of chemically reacting systems. The results obtained have been compared with the results of the calculation of other authors and, where possible, with experimental data. In all cases studied the strategy predicted the correct phase configuration of the reacting system at equilibrium. The component distribution within the phases, however, was found for some of the examples to differ quantitatively, which was a consequence of the different thermodynamic models applied. The overall strategy can be applied successfully to chemically reacting systems with multiple liquid phases. Conclusions The overall strategy presented in this paper advocates a different route to the solution of the simultaneous chemical and phase equilibrium problem. Among its strong points should be mentioned in particular the following: The phase configuration of the reacting system is determined at an early stage by a hybrid technique, which includes a thermodynamic stability analysis and a sequence of two-phase flash calculations upon the candidate solution. This indirect method proved to be extremely reliable, and in addition offered a flexibility and automatic determination of a set of excellent initial estimates for the initialization of the chemical and phase calculations algorithms. The new method for the solution of the chemical equilibrium problem requires modest computational efforts because the number of its independent variables is equal to the number of independent chemical reactions, which makes it particularly appropriate for the cases of chemical reactions involving many components and not too many reactions.
fi ) fugacity of component i in the mixture, Pa Fin i ) initial amount of species i, mol ) final amount of species i, mol Fout i G ) total Gibbs free energy, J G0i ) standard Gibbs free energy of the formation of species i, J ∆G0i ) standard Gibbs free energy change of reaction j, J ki(y) ) function of the chemical potential difference, eq 18 k* ) a number, corresponding to y* KRj ) chemical equilibrium constant for reaction j Kzj ) constant for reaction j defined by eqs 10 and 11 Ki ) K-value for component i Kmax ) the maximum K-value (among all) at equilibrium Kmin ) the minimum K-value (among all) at equilibrium M ) number of independent chemical reactions n ) number of moles N ) number of components in the reacting system NP ) number of phases P ) pressure, Pa S(R) ) error function, eq 9 T ) temperature, K xi ) mole fraction of species i in liquid phase yi ) mole fraction of species i in vapor phase y* ) zero of the functional Φ(y) for the initial system zi ) mole fraction of species i in total system z* ) mole fraction vector of the candidate solution Greek Letters R ) phase split for a two-phase flash φi ) fugacity coefficient, component i Φ(y) ) functional, eq 18 Φ ) system of M equations, eq 16 νij ) stoichiometric coefficient of component i in reaction j ψ ) phase split for the three-phase flash ξj ) extent of reaction j Abbreviations EOS ) equation of state EtAc ) ethyl acetate EtOH ) ethanol HAC ) acetic acid LLE ) liquid-liquid equilibrium LLLE ) liquid-liquid-liquid equilibrium VLE ) vapor-liquid equilibrium VLLE ) vapor-liquid-liquid equilibrium Superscripts final ) final new ) value of the variable obtained at the current iteration old ) value of the variable obtained at the previous iteration L ) liquid V ) vapor * ) corresponding to a zero of the functional Φ(y) Subscripts i ) component index, i ) 1, 2, ... j ) reaction index, j ) 1, 2, ...
Acknowledgment
Appendix A
Dr. Roumiana P. Stateva would like to acknowledge with gratitude the financial support of the Royal Society.
In what follows we shall prove that the extent of the chemical reaction ξfinal is a summation of the extents of the reactions obtained by application of Algorithms I upon z and Algorithm II upon z* (Algorithms I upon z and Algorithm III upon z*). The material balance of the chemical equilibrium is given according to eq 3:
Nomenclature A ) atom matrix (elemental abundance matrix) b ) elemental abundance vector
5482 Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 M
Fout ) Fin i i +
∑ j)1
i ) 1, 2, ..., N
νijξj
j ) 1, 2, ..., M The iteration formula for the calculation of the chemical equilibrium is M
+ Fki ) Fk-1 i
νijξk-1 ∑ j j)1
i ) 1, 2, ..., N (A1)
F0i ) Fin i ξk-1 is the iteration step on the extent of the reaction ξj j and ξ0j , the initial estimate for ξj. The expression for Fki is
Fki ) Fk-1 + i
Fk-2 + ∑νijξk-2 + ∑j νijξk-1 j i j i ) ... ) ∑j νijξk-1 j
) F0i +
) ∑j νijξ0j + ... + ∑j νijξk-1 j )) F0i + ∑νij(ξ0j + ... + ξk-1 j j
) Fin i +
) ∑j νij(ξ0j + ... + ξk-1 j
(A2)
If the iteration process is converged on the kth iteration and in addition Fout ) Fki , then i
Fout ) Fin i i +
) ∑j νij(ξ0j + ... + ξk-1 j
or the value of the reaction extent is obtained according to
ξfinal ) ξ0j + ... + ξk-1 j j
(A3)
It is very important to stress that the manner in which the values of the extents of the reactions are obtained is not discussed since the above derivation does not depend on that. Literature Cited Baker, L. E.; Pierce, A. C.; Luks, K. D. Gibbs Energy Analysis of Phase Equilibria. Soc. Pet. Eng. 1982, 22, 731. Bisset, L. Equilibrium Constants for Shift Reactions. Chem. Eng. 1977, 84, 155. Boston, J. F.; Britt, H. I. A Radically Different Formulation and Solution of the Single-Stage Flash Problem. Comput. Chem. Eng. 1978, 2, 109. Castier, M.; Rasmusen P.; Fredenslund, Aa. Calculation of Simultaneous Chemical and Phase Equilibria in Non-Ideal Systems. Chem. Eng. Sci. 1989, 44, 237. Castillo, J.; Grossmann, I. E. Computation of Phase and Chemical Equilibria. Comput. Chem. Eng. 1981, 5, 99. Cherednichenko, V. M. Ph.D. Thesis, Karpov’s Institute of Physical Chemistry, Moscow, USSR, 1953. Dahl, S.; Michelsen, M. L. High-Pressure Vapour-Liquid Equilibrium with UNIFAC-based Equation of State. AIChE J. 1990, 36, 1829. Dahl, S.; Fredenslund, Aa.; Rasmusen, P. The MHV-2 Model; a UNIFAC-based Equation of State Model for Prediction of Gas
Solubility and Vapor-Liquid Equilibrium at Low and High Pressures. Ind. Eng. Chem. Res. 1991, 30, 1936. Fletcher, R. Practical Methods for Optimization; 2nd ed.; Wiley: New York, 1987. Gautam, R.; Seider, W. D. Computation of Phase and Chemical Equilibria. AIChE J. 1979, 25, 991. Gupta, A. K.; Bishnoi, P. R.; Kalogerakis, N. A Method for the Simultaneous Phase Equilibria and Stability Calculations for Multiphase Reacting and Non-Reacting Systems. Fluid Phase Equilib. 1991, 63, 65. Kantorovich, L. B.; Akilov, G. P. Functional Analysis; Science: Moscow, 1984. Larsen, B. L.; Rasmussen, P.; Fredenslund, Aa. A Modified UNIFAC Group Contribution Model for Prediction of Phase Equilibriua and Heat of Mixing. Ind. Eng. Chem. Res. 1987, 26, 2274. McDonald, C. M.; Floudas, C. A. Global Optimisation for the Phase Stability Problem. AIChE J. 1995a, 41, 1798. McDonald, C. M.; Floudas, C. A. Global Optimisation for the Phase and Chemical Equilibrium Problem. Application to the NRTL Equation. Comput. Chem. Eng. 1995b, 19, 111. McDonald, C. M.; Floudas, C. A. GLOPEQ: A New Computational Tool for the Phase and Chemical Equilibrium Problem. Comput. Chem. Eng. 1997, 21, 1. Michelsen, M. L. The Isothermal Flash Problem. Part I: Stability. Fluid Phase Equilib. 1982, 9, 1. Perez Cisneros, E. S.; Gani, R.; Michelsen, M. L. Reactive Separation SystemssI. Computation of Physical and Chemical Equilibrium. Chem. Eng. Sci. 1997, 52, 527. Reid, R. C.; Prausnitz, J. M.; Sherwood, T. K. The Properties of Gases and Liquids; McGraw-Hill Book Company: New York, 1977. Sanderson, R. V.; Chien, H. H. Y. Simultaneous Chemical and Phase Equilibrium Calculation. Ind. Eng. Chem. Process Des. Dev. 1973, 12, 80. Seider, W. D.; Widagdo, S. Multiphase Equilibria of Reactive Systems. Fluid Phase Equilib. 1996, 123, 283. Seider, W. D.; Brengel, D. D.; Widagdo, S. Nonlinear Analysis in Process Design. AIChE J. 1991, 37, 1. Sherman, S. R.; Eckert, Ch. A.; Scott, L. S. Modeling Multicomponent Equilibria from Binary Equilibrium Data for Reacting Systems. Chem. Eng. Process. 1996, 35, 363. Smith, W. R. The Computation of Chemical Equilibria in Complex Systems. Ind. Eng. Chem. Fundam. 1980, 19, 1. Stateva, R. P.; Tsvetkov, St. G. A Rigorous Approach to Thermodynamic Stability Analysis as a First Step When Solving the Isothermal Multiphase Flash Problem. Technol. Today 1991, 4, 233. Stateva, R. P. Predicting and Calculating Complex Phase Equilibrium in Supercritical Fluid Systems with a New Technique. Collect. Czech. Chem. Commun. 1995, 60, 188. Stateva, R. P.; Tsvetkov, St. G. A Diverse Approach for the Solution of the Isothermal Multiphase Flash Problem. Application to Vapor-Liquid-Liquid Systems. Can. J. Chem. Eng. 1994, 72, 722. Stull, D. R.; Westrum, E. F.; Sinke, G. C. The Chemical Thermodynamics of Organic Compounds; Wiley: New York, 1969. Trebble, M. A.; Bishnoi, P. R. Extension of the Trebble-Bishnoi Equation of State to Fluid Mixtures. Fluid Phase Equilib. 1988, 40, 1. Vidal, J. Thermodynamic: Methods Apliquees au Rafinage et au Genie Chimique; Editions Technip: Paris, 1974; Vol. 2. Whitson, C. H.; Michelsen, M. L. The Negative Flash. Fluid Phase Equilib. 1989, 53, 51. Xiao, W.; Zhu, K.; Yuan, W.; Chien H. H. Algorithm for Simultaneous Phase and Chemical Equilibrium Calculations. AIChE J. 1989, 35, 1813.
Received for review April 4, 1997 Revised manuscript received August 29, 1997 Accepted August 29, 1997X IE9702643 X Abstract published in Advance ACS Abstracts, November 1, 1997.
