At the end of Section 15.2, the bubble-pressure method was briefly introduced to show how the fugacity coefficients are incorporated into a VLE calculation, without concentrating on the details of the iterations. Section 15.3 offered derivations of formulas for the fugacity coefficients that were presented without proof at the beginning of the chapter. Now, it is time to turn to the applied engineering objective: calculation of phase equilibria. Refer again to Table 10.1 on page 373, that lists the types of routines that are needed and the convergence criteria. Note that Table 10.1 is independent of the model used for calculating VLE. As an example of the iteration procedure for cubic equations of state, the bubble-pressure flow sheet is presented in Fig. 15.2. The flow sheet puts detail to the procedure discussed superficially in Example 15.2 and immediately preceding the example. Flow sheets for bubble temperature, dew, and flash routines are available in Appendix C. As with ideal solutions, the bubble-pressure routine is the easiest to apply, so we cover it in detail in the following examples. Iterative phase equilibrium calculations can be tedious and difficult to automate. We can facilitate the calculations to some extent by combining two copies of the PrFug spreadsheet into a single workbook, which we call Prmix.xlsx. The four possible K-value representations are included for convenient selection, as described in Example 15.2. This workbook forms only a starting basis with an emphasis on clearly showing the fundamental steps.
Figure 15.2. Bubble-pressure flow sheet for the equation of state method of representing VLE. Other routines are given in Appendix C.
The engineering objective is to use equations of state for bubble, dew, and flash calculations.
Flow sheets for bubble temperature, dew, and flash routines are in Appendix C.
Bubble-pressure calculations are enabled with the spreadsheet Prmix.xlsx.
Example 15.6. Bubble-point pressure from the Peng-Robinson equation
Use the Peng-Robinson equation (kij = 0) to determine the bubble-point pressure of an equimolar solution of nitrogen (1) + methane (2) at 100 K.
Solution
The calculations proceed by first calculating the short-cut K-ratio as in Example 15.2 on page 587. The ideal-solution (is) bubble pressure was 
; yN2is = 0.958. In fact, the K-values for the first iteration have already been determined in that example, in great detail. The values from that example are K1 = 1.955, K2 = 0.1106. The new estimates of vapor mole fractions are obtained by multiplying xi·Ki. In Example 15.2, the sum was found to be 
. These calculations are summarized in the first column of Table 10.1.
Noting that these sum to a number greater than unity, we must choose a greater value of pressure for the next iteration. Before we can start the next iteration, however, we must develop new estimates of the vapor mole fractions; the ones we have do not make sense because they sum to more than unity. These new estimates can be obtained simply by dividing the given vapor mole fractions by the number to which they sum. This process is known as normalization of the mole fractions. For example, to start the second iteration, y1 = 0.978/1.033 = 0.947. After repeating the process for the other component, the mole fractions will sum to unity. Since the result for the second iteration is less than one, the pressure guess is too high.
Normalization of mole fractions.
The third iteration consists of applying the interpolation rule to obtain the estimate of pressure and use of the normalization procedure to obtain the estimates of vapor mole fractions. P = 0.4119 + (1 – 1.033)/(0.956 – 1.033) · (0.45 – 0.4119) = 0.428 MPa. Since the estimated vapor mole fractions after the third iteration sum very nearly to unity, we may conclude the calculations here. This is the bubble pressure. Note how quickly the estimate for y1 converges to the final estimate of 0.945.
Example 15.7. Isothermal flash using the Peng-Robinson equation
A distillation column is to produce overhead products having the following compositions:
Suppose a partial condenser is operating at 320 K and 8 bars. What fraction of liquid would be condensed according to the Peng-Robinson equation, assuming all binary interaction parameters are zero (kij = 0)?
Solution
This is an isothermal flash calculation. Refer back to the same problem (Example 10.1 on page 382) for an initial guess based on the shortcut K-ratio equation. V/F = 0.25 ⇒ {xi} = {0.1829, 0.7053, 0.1117} and {yi} = {0.3713, 0.5642, 0.0648}. Substituting these composition estimates for the vapor and liquid compositions into the routine for estimating K-values (cf. Example 15.2 on page 587), we can obtain the estimates for K-values given below:
The computations for the flash calculation are basically analogous to those in Example 10.1, except that Ki values are calculated from Eqn. 15.20. A detailed flow sheet is presented in Appendix C. For this example, the K-values are not modified until the iteration on V/F converges. After convergence on V/F, the vapor and liquid mole fractions are recomputed using Eqns. 10.15 and 10.16, followed by recomputed estimates for the K-values. If the new estimates for K-values are equal to the old estimates for K-values, then the overall iteration has converged. If not, then the new estimates for K-values are substituted for the old values, and the next iteration proceeds just like the last. This method of iteratively solving for the vector of K-values is known in numerical analysis as the “successive substitution” method.
Using these x’s and y’s for guesses we find K = 1.7276, 0.8318, and 0.6407, respectively. These K-values are similar to those estimated at the compositions derived from the ideal-solution approximation, and will yield a similar V/F. Therefore, we conclude that this iteration has converged (a general criterion is that the average % change in the K-values from one iteration to the next is less than 10-4). Comparison to the shortcut K-ratio approximation shows small but significant deviations—V/F = 0.13 for Peng-Robinson versus 0.25 for the shortcut K-ratio method.
Based on this example, we may conclude that the shortcut K-Ratio approximation provides a reasonable first approximation at these conditions. Note, however, that none of the components is supercritical and all the components are saturated hydrocarbons.
It is tempting to expand further on Prmix.xlsx to facilitate greater automation and simple-minded application. An online supplement provides a very preliminary step in this direction through the use of macro’s. Ultimately, however, this literature comprises specialized research that is beyond our introductory scope. In general, the analysis requires detailed consideration of phase stability and criticality. References cited in Chapter 16 describe works by Michelsen and Mollerup, Eubank, and Tang that can help to create more reliable algorithms. It is a useful exercise to customize your workbooks to increase your confidence in achieving reliable solutions, but do not spend excessive time trying to program every possibility. Chapter 16 and the references cited there are recommended for advanced programming.
Example 15.8. Phase diagram for azeotropic methanol + benzene
Methanol and benzene form an azeotrope. For methanol + benzene the azeotrope occurs at 61.4 mole% methanol and 58°C at atmospheric pressure (1.01325 bars). Additional data for this system are available in the Chemical Engineers’ Handbook. Use the Peng-Robinson equation with kij = 0 (see Eqn. 15.9) to estimate the phase diagram for this system and compare it to the experimental data on a T-x-y diagram. Determine a better estimate for kij by iterating on the value until the bubble point pressure matches the experimental value (1.013 bar) at the azeotropic composition and temperature. Plot these results on the T-x-y diagram as well. Note that it is impossible to match both the azeotropic composition and pressure with the Peng-Robinson equation because of the limitations of the single parameter, kij.
Prmix.xlxx can be used to fit the kij.
The experimental data for this system are as follows:
Solution
Solving this problem is computationally intensive, but still approachable with Prmix.xlsx. The strategy is to manually set a guessed kij and then perform a bubble pressure calculation at the azeotrope temperature (331.15 K) and composition, xm = 0.614. The program will give a calculated pressure and vapor phase composition. The vapor-phase composition may not match the liquid-phase composition because the azeotrope is not perfectly predicted; however, we continue to manually change kij, and repeat the bubble pressure calculation until we match the experimental pressure of 1.013 bar. The following values are obtained for the bubble pressure at the experimental azeotropic composition and temperature with various values of kij.
The resultant kij is used to perform bubble-temperature calculations across the composition range, resulting in Fig. 15.3. Note that we might find a way to fit the data more accurately than the method given here, but any improvements would be small relative to estimating kij = 0. We see that the fit is not as good as we would like for process design calculations. This solution is so nonideal that a more flexible model of the thermodynamics is necessary. Note that the binary interaction parameter alters the magnitude of the bubble-pressure curve very effectively but hardly affects the skewness at all. Since this mixture is far from the critical region, a two-parameter activity model like van Laar or UNIQUAC would be recommended as shown in Fig. 12.1. The Peng-Robinson model with van der Waals’ mixing rule comes closest to the Scatchard-Hildebrand activity model. We observed that the Scatchard-Hildebrand model performed poorly for hydrogen bonding mixtures. Two approaches are common when precision is needed in the critical region for hydrogen bonding mixtures. Either a multiparameter activity model can be adapted as a basis for an advanced mixing rule, or hydrogen bonding can be treated explicitly. References to these are approaches are presented in the chapter summary.
Figure 15.3. T-x-y diagram for the azeotropic system methanol + benzene. Curves show the predictions of the Peng-Robinson equation (kij = 0) and correlation (kij = 0.084) based on fitting a single data point at the azeotrope; x’s and triangles represent liquid and vapor phases, respectively.
Example 15.9. Phase diagram for nitrogen + methane
Use the Peng-Robinson equation (kij = 0) to determine the phase diagram of nitrogen + methane at 150 K. Plot P versus x, y and compare the results to the results from the shortcut K-ratio equations.
Solution
First, the shortcut K-ratio method gives the dotted phase diagram in Fig. 15.4. Applying the bubble-pressure procedure with the program Prmix.xlsx, we calculate the solid line in Fig. 15.4. For the Peng-Robinson method we assume K-values from the previous solution as the initial guess to get the solutions near xN2 = 0.685. The program Prmix.xlsx assumes this automatically, but we must also be careful to make small changes in the liquid composition as we approach the critical region.
Figure 15.4. High pressure P-x-y diagram for the nitrogen + methane system comparing the shortcut K-ratio approximation and the Peng-Robinson equation at 150 K. The data points represent experimental results. Both theories are entirely predictive since the Peng-Robinson equation assumes that kij = 0.
Fig. 15.4 was generated by entering liquid nitrogen compositions of 0.10, 0.20, 0.40, 0.60, 0.61, 0.62…, 0.68, and 0.685. This procedure of starting in a region where a simple approximation is reliable and systematically moving to more difficult regions using previous results is often necessary and should become a familiar trick in your accumulated expertise on phase equilibria in mixtures. We apply a similar approach in estimating liquid-liquid equilibria.
Comparing the two approximations numerically and graphically, it is clear that the shortcut approximation is significantly less accurate than the Peng-Robinson equation at high concentrations of the supercritical component. This happens because the mixture possesses a critical point, above which separate liquid and vapor roots are impossible, analogous to the situation for pure fluids. Since the mixing rules are in terms of a and b instead of Tc and Pc, the equation of state is generating effective values for Ac and Bc of the mixture.
Instead of depending simply on T and P as they did for pure fluids, however, Ac and Bc also depend on composition. The mixture critical point varies from one component to the other as the composition changes. Since the shortcut (and also SCVP+) approximation extrapolates the vapor pressure curve to obtain an effective vapor pressure of the supercritical component, that approximation does not reflect the presence of the mixture critical point and this leads to significant errors as the mixture becomes rich in the supercritical component.
The mixture critical point also leads to computational difficulties. If the composition is excessively rich in the supercritical component, the equation of state calculations may obtain the same solution for the vapor root as for the liquid root and, since the fugacities are equal, the program will terminate. The program may indicate accurate convergence in this case due to some slight inaccuracies that are unavoidable in the critical region. Or the program may diverge. It is often up to the competent engineer to recognize the difference between accurate convergence and a spurious answer. Plotting the phase envelope is an excellent way to stay out of trouble.
The shortcut K-ratio method provides an initial estimate when a supercritical component is at low liquid-phase compositions, but incorrectly predicts VLE at high liquid-phase concentrations of the supercritical component.
Example 15.10. Ethane + heptane phase envelopes
Use the Peng-Robinson equation (kij = 0) to determine the phase envelope of ethane + n-heptane at compositions of xC7 = [0, 0.1, 0.2, 0.3, 0.5, 0.7, 0.9, 1.0]. Plot P versus T for each composition by performing bubble-pressure calculations to their terminal point and dew-temperature calculations until the temperature begins to decrease significantly and the pressure approaches its maximum. If necessary, close the phase envelope by starting at the last dew-temperature state and performing dew-pressure calculations until the temperature and pressure approach the terminus of the bubble-point curve. For each composition, mark the points where the bubble and dew curves meet with X’s. These X’s designate the “mixture critical points.” Connect the X’s with a dashed curve. The dashed curve is known as the critical locus of the mixture.
Prmix.xlsx facilitates bubble, dew, and flash calculations.
Solution
Note that these phase envelopes are similar to the one from the previous problem, except that we are changing the temperature instead of the composition along each curve. They are more tedious in that both dew and bubble calculations must be performed to generate each curve. The lines of constant composition are sometimes called isopleths. The results of the calculations are illustrated in Fig. 15.5. The results at mole fractions of 0 and 1.0 are indicated by dash-dot curves to distinguish them as the vapor pressure curves. Phase equilibria on the P-T plot occurs at the conditions where a bubble line of one composition intersects a dew line of a different composition.
Figure 15.5. High pressure phase envelopes for the ethane + heptane system comparing the effects of composition according to the Peng-Robinson equation. The theory is entirely predictive because the Peng-Robinson equation has been applied with kij = 0. X’s mark the mixture critical points and the dashed line indicates the critical locus. The first and last curves represent the vapor pressures for the pure components.
Some practical considerations for high pressure processing can be inferred from the diagram. Consider what happens when starting at 90 bars and ~445 K and dropping the pressure on a 30 mole% C7 mixture at constant temperature. Similar situations could arise with flow of natural gas through a small pipe during natural gas recovery. As the pressure drops, the dew-point curve is crossed and liquid begins to condense. Based on intuition developed from experiences at lower pressure, one might expect that dropping the pressure should result in more vapor-like behavior, not condensation. On the other hand, dropping the pressure reduces the density and solvent power of the ethane-rich mixture. This phenomenon is known as retrograde condensation. It occurs near the critical locus when the operating temperature is less than the maximum temperature of the phase envelope. Since this maximum temperature is different from the mixture critical temperature, it needs a distinctive name. The name applied is the “critical condensation temperature” or cricondentherm. Similarly, the maximum pressure on the phase envelope is known as the cricondenbar. Note that an analogous type of phase transition can occur near the critical locus when the pressure is just above the critical locus and the temperature is changed.
To extend the analysis, imagine what happens in a natural gas stream composed primarily of methane but also containing small amounts of components as heavy as C80. A retrograde condensation region exists where the heavy components begin to precipitate, as discussed in Example 14.12 on page 567. But a different possibility also exists because the melting temperature of the heavy components may often exceed the operating temperature, and the precipitate that forms might be a solid that could stick to the walls of the pipe. This in turn generates a larger constriction which generates a larger pressure drop during flow, right in the vicinity of the deposit. In other words, this deposition process may tend to promote itself until the flow is substantially inhibited. Wax deposition is a significant problem in the oil and natural gas industry and requires considerable engineering expertise because it often occurs away from critical points, as well as in the near-critical regions of this discussion. A wide variety of solubility behavior can occur, as we show in Chapter 16.
Leave a Reply