3530
Ind. Eng. Chem. Res. 1999, 38, 3530-3534
GENERAL RESEARCH A Method For Avoiding Trivial Roots in Isothermal Flash Calculations Using Cubic Equations of State G. V. Pasad and G. Venkatarathnam* Refrigeration and Airconditioning Laboratory, Department of Mechanical Engineering, Indian Institute of Technology, Chennai 600 036 India
The variation of isothermal compressibility with temperature is unique for gases and liquids. Based on this relationship, a simple method is presented for identifying trivial roots which occur during isothermal flash calculations. The usefulness of the method are demonstrated for some mixtures. The results obtained with the proposed method are compared with those available in literature. 1. Introduction The simplicity of the cubic equations of state such as the Peng-Robinson, Soave-Redlich-Kwong, and Harmen-Knapp have made them popular for determining the properties of both the liquid and vapor phases of mixtures even at high pressures.1 It is also well-known that a single root is obtained with the cubic equations of state even within the vapor dome of some mixtures, particularly at high pressures.2 Isothermal flash calculations performed using the cubic equations of state can sometimes predict the wrong phase because of the poor initial approximation of the phase compositions (K values). Reid et al.2 suggest that the best way to avoid the problem is to make sufficiently accurate initial guesses in iterative calculations. Sometimes the calculations may converge to a physically nonexistent condition when there is a large difference between the K values used at the start of iterations and the actual K values.3 The problem of converging to physically nonexistent conditions is generally known as the trivial root problem. The problem of physically nonexistent convergence and trivial roots has been studied by many authors. Gundersen4 showed that the Soave-Redlich-Kwong equation of state becomes monotonic at moderate and high pressures. He showed that this may happen for one or both the phases well inside the two-phase region. Gundersen suggested that the maximum and minimum points be used instead of the actual volume (compressibility factor) roots of the equation of state. He proposed an algorithm for solving the cubic equations to overcome the trivial root problem. Gundersen compared the results obtained with his method with that obtained using MINSIM program developed at the Oklahoma State University for a eight component mixture of nitrogen and different hydrocarbons. Gundersen states that the algorithm fails outside the two-phase region when the pressure is greater than 60% of the true * To whom all correspondence should be addressed. Email:
[email protected].
critical pressure. The critical region could be reached only by sequential calculations starting at low pressures (1 atm). Poling et al.1 studied the problem of trivial roots and spurious derivatives and proposed a simple method for finding trivial roots. They presented the variation of molar enthalpies and dK/dT with temperature for an equimolar mixture of methane, ethylene, and ethane at 30 atm. They showed that negative values are predicted for dK/dT and for the heat of vaporization between 244 and 250 K, because of the occurrence of trivial roots. They suggested that the compressibility factor (Z) be considered satisfactory for determining liquid and vapor phases below the pseudocritical pressure if the isothermal compressibility of the liquid and vapor phases satisfy the following criteria:
βL < 0.005 atm-1
(1)
3.0 0.9 < βV < P P
(2)
Poling et al.1 state that there is no theoretical significance for eqs 1 and 2; however, real liquid and vapor phases satisfy the conditions for a wide range of pressure and temperature. They recommend that the compositions be adjusted to move out of the trivial root region when liquid-phase compositions are desired and to reduce the pressure to generate the correct vaporphase properties. The adjustment procedures are used to obtain the initial estimates for use in early iterative calculations. The adjustments are made until the assumed phase compositions are close to true values and conditions 1 and 2 are satisfied. Poling et al.1 found that the success of the method is not strongly dependent on the specific limits on isothermal compressibility for the examples presented by them. When 0.005 in eq 1 was changed to 0.03, the method was still successful for the distillation calculations presented by them. The fundamental equations that are used in the flash calculations are (a) the fugacity equality criteria and (b) the material balance relations. Different methods can
10.1021/ie980661t CCC: $18.00 © 1999 American Chemical Society Published on Web 09/07/1999
Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 3531
be used for updating the fugacity and the composition of the two phases in the iterative calculations. Veeranna et al.3 studied two different methods for solving the flash calculations. They showed that the scheme where the material balance equation is in the inner loop and the fugacity equality is checked in the outer loop is sensitive to the initial K-values. On the other hand, the scheme where the material balance equation is in the outer loop and the fugacity equality is in the inner loop is sensitive to the initial vapor fraction assumed. They studied the computational difficulties arising out of (a) trivial solutions, (b) deviation function becoming very small leading to false convergences, (c) physically nonexistent convergences, and (d) stationary points on the objective function (material balance relation). According to Veeranna et al.,3 the above problems arise because of the nonlinear nature of the resulting mathematical equations which can also be satisfied by some physically nonexistent solutions. They proposed an algorithm for avoiding trivial roots and physically nonexistent converges. In their method, they predict the correctness of the solution by comparing the density of the phases with the critical density as follows:
Fsat,V < Fc
(3)
Fsat,L > Fc
(4)
They proposed a flash algorithm based on the above method. They calculated the vapor fraction of the eightcomponent nitrogen-hydrocarbon mixture studied by Gundersen4 using the Soave-Redlich-Kwong equation of state using their proposed flash algorithm. Their results match closely with that of Gundersen,4 except at 172 atm, where there is a small deviation. By use of the critical density criteria and the flash algorithm proposed by Veeranna et al.,3 it is possible to reach correct solutions, irrespective of the starting K-values. The method, however, requires correct estimation of critical density. The correct critical density of a mixture can be estimated by using either rigorous or empirical methods.2 The computation time required to determine the true critical point using the rigorous method is not trivial. The empirical methods, on the other hand, have been tested largely with hydrocarbon mixtures, and the accuracy of the predicted values is not known when the methods are applied to other fluids.2 Another approach that can be used to avoid trivial roots is to use pseudo roots for an unfeasible phase to evaluate K-values, as proposed by Mathias et al.5 The method is based on a heuristic approach and aims at maintaining continuity in the pseudo property and its derivatives to aid the convergence of iterations in process calculations. Very recently, Zhao and Saha6 extended the method of Mathias et al.5 They presented analytical expressions for implementing their method, which is based on the variation of complex roots with pressure for a cubic equation of state at constant temperature and composition. In this paper, a simple method is proposed to identify the trivial vapor and liquid roots. The method has been used to predict the phase split of some fluid mixtures. The results obtained are compared with those available in literature. 2. Analysis The isothermal compressibility of liquids is, in general, about an order of magnitude lower than that of
Figure 1. Variation of isothermal compressibility with reduced temperature and pressure for a van der Waals fluid.
Figure 2. Variation of isothermal compressibility with reduced temperature and pressure for nitrogen.
the vapor at pressures well below the critical pressure. The isothermal compressibility criteria of Poling et al.1 is thus an excellent criteria for identifying the correct vapor and liquid roots. Poling et al.’s conditions (eqs 1 and 2), however, have no theoretical basis. The isothermal compressibility of real gases is quite different from that of the ideal value (βideal ) 1/P). The deviation of isothermal compressibility of a real fluid from that of an ideal fluid is a complex function of the pressure and temperature, as well as the nature of fluid. However, the variation of isothermal compressibility with temperature is quite unique for all gases and gas mixtures, and can be used to identify the correct phase at all pressures and temperatures below the critical pressure. Figure 1 shows the variation of Pβ with reduced temperature for a van der Waals fluid at different reduced pressures. It can be observed from Figure 1 that the nature of variation of isothermal compressibility with temperature is different for liquids and gases. Therefore the term (∂β/∂T)P of the vapor will be negative and that of the liquid phase positive for pure fluids at all temperatures and pressures less than the critical pressure. Figure 2 shows the variation of Pβ with reduced temperature for nitrogen at different reduced temperatures and pressures, evaluated using the PengRobinson equation of state. It is evident from Figures 1 and 2 that the nature of variation of the isothermal compressibility with temperature is quite unique and different for liquids and gases and is independent of the fluid and the type of equation of state used. The
3532 Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999
Figure 3. Variation of (∂β/∂T)P with temperature for nitrogen above the critical pressure.
following generalizations can now be made from the above study:
(∂β/∂T)P,L > 0
(5)
(∂β/∂T)P,V < 0
(6)
The above expressions are valid at all temperatures and pressure below the critical pressure only since the isothermal compressibility-temperature curves show points of inflection above the critical pressure, as shown in Figure 3. It must be remembered that the isothermal compressibility becomes infinite at a pressure below the true critical pressure for some mixtures.7 This point is known as the mechanical critical point. The expressions are therefore applicable to mixtures only below the mechanical critical pressure. Figure 4 shows the variation of (∂β/∂T)P with temperature for a mixture of methane (36.2 mol %), ethane (9%), ethylene (37.4%), propane (1.7%), and propylene (15.7%) at different reduced pressures below the pseudocritical pressure. The different parameters were obtained using the Peng-Robinson equation of state, and the binary interaction parameters presented by Knapp et al.8 It can be observed from Figure 4 that the general nature of the variation of β with temperature for the mixture is the same as that for pure fluids shown in Figures 1 and 2. The points of inflection seen in Figure 4 corresponds to the bubble and dew point temperatures, which are due to the difference in the composition of the phases in the single and the two-phase (equilibrium) conditions. The overall nature, however, is the same for mixtures and pure fluids in both the single and two-phase conditions. It can now be concluded that the roots of a cubic equation are nontrivial only when they satisfy the conditions listed in eqs 5 and 6. The criteria for detecting the trivial vapor and liquid roots for pressure below the mechanical critical pressure can now be expressed as follows. Trivial liquid root:
(∂β/∂T)P,L < 0
(7)
Trivial vapor root:
(∂β/∂T)P,V > 0
(8)
The isothermal compressibility can be evaluated in
Figure 4. Variation of (∂β/∂T)P with temperature and pressure for mixture M1. Table 1. Steps in Algorithm for Isothermal Flash Calculations 1. Estimate initial K-values (ideal). 2. Assume initial vapor fraction of 0.5. 3. Calculate the phase composition, component fugacities, and new K-values. 4. Repeat step 3 until (knew - kold) < allowable error. 5. If (∂β/∂T)V > 0.0, decrease a vapor fraction; reinitalize K-values; go to step 3. 6. If (∂β/∂T)L < 0.0, increase a vapor fraction; reinitialize K-values; go to step 3. 7. If (∑yi - ∑xi) > 0.0, increase vapor fraction; go to step 3. 8. If (∑yi - ∑xi) < 0.0, decrease vapor fraction; go to step 3.
terms of pressure and temperature as follows:
β)-
1 ∂v 1 P ∂Z ) 1v ∂P T P Z ∂P T
( )
[
( )]
(9)
The term (∂β/∂T)P can be expressed directly in terms of pressure and temperature as follows:
(∂T∂β)
P
∂ 1 ∂v ) βv ∂v v ∂T P T
[( )]
(10)
The determination of the partial derivative (∂v/∂T)P is quite straightforward for all cubic equations of state. Analytical expressions can be obtained easily for (∂β/ ∂T)P since (∂v/∂T)P is a function of specific volume and temperature only. Table 1 shows the algorithm for the isothermal flash calculations that incorporates the criteria for detecting trivial roots. In our algorithm, the vapor fraction is assumed to be 0.50 at the start of the iterations. Ideal K-values (Pvp/P) are used to start the iteration. The K-values are updated continuously until the fugacity equality condition is satisfied. At this point, the possible occurrence of trivial volume roots of the two phases is checked using conditions specified by eqs 7 and 8. If a trivial vapor root is detected, the vapor fraction is decreased or the calculation is moved toward the liquid region. If a trivial liquid root is detected, the vapor fraction is increased or the calculation is moved toward the vapor region. When the calculation moves out of the trivial root region, the algorithm checks for component balance and adjusts the vapor fraction to achieve the final convergence. In the above method, the conventional method of assigning the maximum of the three volume roots to the vapor phase and the minimum to the liquid phase has been used. The main strength of the proposed algorithm is the negligible computational
Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 3533
Figure 5. Variation of objective function F(VF) with vapor fraction for mixture M1 at 30 bar and 215 K.
Figure 6. Variation of (∂β/∂T)P with vapor fraction for mixture M1 at 30 bar and 215 K.
Table 2. Compostion of Mixtures Studied (Mole Fraction) fluid nitrogen methane ethane ethylene propane propylene n-butane n-pentane n-hexane n-heptane
mixture M1 0.362 0.09 0.374 0.017 0.157
mixture M2 0.0054 0.728 0.0546 0.0306 0.0307 0.0688 0.0438 0.0375
overhead in checking for trivial roots and the ease with which the same can be incorporated into traditional flash algorithms. Unlike Poling et al.’s criteria, the method is based on a fundamental property of liquids and gases and can be applied to all mixtures. The approach presented here has been used to study the occurrence of trivial roots in some mixtures (see Table 2); the results of which are described in the next section.
Figure 7. Variation of objective function F(VF) with vapor fraction for mixture M2 at 32 bar and 340 K.
3. Results and Discussion Figure 5 shows the variation of objective function F (F ) ∑yi - ∑xi) at different vapor fractions for mixture M1 at 30 bar and 215 K. The fugacity equality criteria is satisfied at all of the vapor fractions in Figure 5. In this case, trivial roots occur between a vapor fraction of 0.73 and 1.0. Normal flash algorithms can predict the solution correctly when the vapor fraction assumed is between 0 and 0.73, but fail in the trivial root region. Figure 6 shows the variation of (∂β/∂T)P with vapor fraction. The term (∂β/∂T)P of the vapor phase is positive between a vapor fraction of 0.73 and 1.0, indicating trivial vapor roots in that region (eq 8). The algorithm described in Table 1 moves the calculation into the nontrivial region to arrive at the correct vapor fraction of 0.01. Figures 7 and 8 show the variation of objective function F and (∂β/∂T)P, respectively, with the vapor fraction for mixture M2 at a pressure of 32 bar and 340 K. It can be seen that the occurrence of trivial liquid root between vapor fractions of 0.0 and 0.44 can be correctly predicted using the condition specified by eq 7. The proposed algorithm increases the vapor fraction to cross the trivial root region, leading to the correct solution. The solution obtained using our algorithm (VF ) 0.866) is same as that obtained by other workers3,4 using their algorithms. The results obtained with the
Figure 8. Variation of (∂β/∂T)P with vapor fraction for mixture M2 at 32 bar and 340 K.
algorithm presented in Table 1 coincides with that presented by other workers3,4 for mixture M2 at all pressures below the pseudocritical pressure. It must, however, be pointed out that both Veeranna et al.3 and Gundersen4 did not use binary interaction parameters in their computations. The results obtained by us with the Peng-Robinson equation of state match with that obtained by the authors3,4 using the Soave-RedelichKwong equation of state only when a similar assumption is made, particularly at high pressures. The proposed algorithm, however, fails to predict the correct roots of mixture M2 for pressures above the mechanical critical pressure and more sophisticated methods need to be employed.
3534 Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999
4. Conclusions
Subscripts
The occurrence of trivial vapor and liquid roots can be correctly identified using the sign of (∂β/∂T)P using eqs 7 and 8, below the mechanical critical pressure. The flash algorithm presented leads to correct solutions for mixtures at all pressures below the mechanical critical pressure, with very little computational overhead for checking the occurrence of trivial roots. To simplify the computations, pseudocritical pressure was used in our calculations in place of mechanical critical pressure.
c ) critical i ) component identifier L ) liquid phase r ) reduced sat ) saturated V ) vapor phase vp ) vapor pressure
Acknowledgment The authors are grateful to Profs. S. Srinivasa Murthy and M. S. Ananth, Indian Institute of Technology, Chennai, and Prof. Dr.-Ing. L. R. Oellrich, Institut fu¨r Technische Thermodynamik und Ka¨ltetechnic, Universita¨t Karlsruhe, Germany for reviewing the work and suggesting some useful changes. Nomenclature F ) objective function in flash calculations (F ) ∑yi - ∑xi) K ) equilibrium constant for component i P ) pressure T ) temperature x ) mole fraction of the component in the liquid phase y ) mole fraction of the component in the vapor phase Z ) compressibility factor
Literature Cited (1) Poling, B. E.; Grens, E. A.; Prausnitz, J. M. Ind. Eng. Chem. Process Des. Dev. 1981, 20, 12. (2) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases and Liquids, 4th ed.; McGraw-Hill: New York, 1989. (3) Veeranna, D.; Husain, A.; Subrahmanyam, S.; Sarkar, M. K. Comput. Chem. Eng. 1987, 11, 489. (4) Gundersen, T. Comput. Chem. Eng. 1982, 6, 245. (5) Mathias, P. M.; Boston, J. F.; Watanasiri, S. AIChE J. 1984, 30, 182. (6) Zhao, E.; Saha, S. Ind. Eng. Chem. Res. 1998, 37, 1625. (7) Rowlinson, J. S. Liquids and Liquid Mixtures, 2nd ed.; Butterworth: London, 1971. (8) Knapp, H.; Doring, R.; Oellrich, L.; Plocker, U.; Prausnitz, J. M. Vapor-Liquid Equilibria for Mixtures of Low Boiling Substances; DECHEMA Chemistry Data Series, Vol. VI; DECHEMA: Frankfurt, 1981.
Received for review October 19, 1998 Revised manuscript received May 17, 1999 Accepted May 28, 1999
Greek Symbol β ) isothermal compressibility
IE980661T
