The only thermodynamic specification that is required for determining the saturation temperature or pressure is that the Gibbs energies (or fugacities) of the vapor and liquid be equal. This involves finding the pressure or temperature where the vapor and liquid fugacities are equal. The interesting part of the problem comes in computing the saturation condition by iterating on the temperature or pressure.
Phase equilibria involves finding the state where fL = fV.
Example 9.6. Vapor pressure from the Peng-Robinson equation
Use the Peng-Robinson equation to calculate the normal boiling point of methane.
Solution
Vapor pressure calculations are available in Preos.xlsx and PreosPropsmenu.m. Let us discuss Preos.xlsx first. The spreadsheet is more illustrative in showing the steps to the calculation. Computing the saturation temperature or pressure in Excel is rapid using the Solver tool in Excel.
On the spreadsheet shown in Fig. 9.7, cell H12 is included with the fugacity ratio of the two phases; the cell can be used to locate a saturation condition. Initialize Excel by entering the desired P in cell B8, in this case 0.1 MPa. Then, adjust the temperature to provide a guess in the two-phase (three-root) region. Then, instruct Solver to set the cell for the fugacity ratio (H12) to a value of one by adjusting temperature (B7), subject to the constraint that the temperature (B7) is less than the critical temperature.
Figure 9.7. Example of Preos.xlsx used to calculate vapor pressure.
In MATLAB, the initial guess is entered in the upper left. The “Run Type” is set as a saturation calculation. The “Root to use” and “Value to match” are not used for a saturation calculation. The drop-down box “For Matching…” is set to adjust the temperature. The results are shown in Fig. 9.8.
Figure 9.8. Example of PreosPropsMenu.m used to calculate vapor pressure.
For methane the solution is found to be 111.4 K which is very close to the experimental value used in Example 8.9 on page 320. Saturation pressures can also be found by adjusting pressure at fixed temperature.
Fugacity and P-V isotherms for CO2 as calculated by the Peng-Robinson equation are shown in Fig. 9.9 and Fig. 7.5 on page 264. Fig. 9.9 shows more clearly how the shape of the isotherm is related to the fugacity calculation. Note that the fugacity of the liquid root at pressures between B and B′ of Fig. 9.6 is lower than the fugacity of the vapor root in the same range, and thus is more stable because the Gibbs energy is lower. Analogous comparisons of vapor and liquid roots at pressures between C and C′ show that vapor is more stable. In Chapter 7, we empirically instructed readers to use the lower fugacity. Now, in light of Fig. 9.9, readers can understand the reasons for the use of fugacity.
Figure 9.9. Predictions of the Peng-Robinson equation of state for CO2: (a) prediction of the P-V isotherm and fugacity at 280 K; (b) plot of data from part (a) as fugacity versus pressure, showing the crossover of fugacity at the vapor pressure. Several isotherms for CO2 are shown in Fig. 7.5 on page 264.
The term “fugacity” was defined by G. N. Lewis based on the Latin for “fleetness,” which is often interpreted as “the tendency to flee or escape,” or simply “escaping tendency.” It is sometimes hard to understand the reasons for this term when calculating the property for a single root. However, when multiple roots exist as shown in Fig. 9.9, you may be able to understand how the system tries to “escape: from the higher fugacity values to the lower values. This perspective is especially helpful in mixtures, indicating the direction of driving forces to lower fugacities.
A stable system minimizes its Gibbs energy and its fugacity.
Just as the vapor pressure estimated by the shortcut vapor pressure equation is less than 100% accurate, the vapor pressure estimated by an equation of state is less than 100% accurate. For example, the Peng-Robinson equation tends to yield about 5% average error over the range 0.4 < Tr < 1.0. This represents a significant improvement over the shortcut equation. The van der Waals equation, on the other hand, yields much larger errors in vapor pressure. One problem is that the van der Waals equation offers no means of incorporating the acentric factor to fine-tune the characterization of vapor pressure. But the problems with the van der Waals equation go much deeper, as illustrated in the example below.
Adapting Preos.xlsx to a different equation of state.
Example 9.7. Acentric factor for the van der Waals equation
To clarify the problem with the van der Waals equation in regard to phase-equilibrium calculations, it is enlightening to compute the reduced vapor pressure at a reduced temperature of 0.7. Then we can apply the definition of the acentric factor to characterize the vapor pressure behavior of the van der Waals equation. If the acentric factor computed by the van der Waals equation deviates significantly from the acentric factor of typical fluids of interest, then we can quickly assess the magnitude of the error by applying the shortcut vapor-pressure equation. Perform this calculation and compare the resulting acentric factor to those on the inside covers of the book.
Solution
The computations for the van der Waals equation are very similar to those for the Peng-Robinson equation. We simply need to derive the appropriate expressions for a0, a1, and a2, that go into the analytical solution of the cubic equation: Z3+ a2 Z2 + a1 Z + a0 = 0.
Adapting the procedure for the Peng-Robinson equation given in Section 7.6 on page 263, we can make Eqn. 7.13 dimensionless:
where the dimensionless parameters are given by Eqns. 7.21–7.24; A = (27/64) Pr/Tr2; B = 0.125 Pr/Tr.
After writing the cubic in Z, the coefficients can be identified: a0 = –AB; a1 = A; and a2 = –(1 + B). For the calculation of vapor pressure, the fugacity coefficient for the van der Waals equation is quickly derived as the following:
Substituting these relations in place of their equivalents in Preos.xlsx, the problem is ready to be solved. Since we are not interested in any specific compound, we can set Tc = 1 and Pc = 1, Tr = 0.7. Setting an initial guess of Pr = 0.1, Solver gives the result that Pr = 0.20046.
Modification of PreosPropsMenu.m is accomplished by editing the routine PreosProps.m. Search for the text “global constants.” Change the statements to match the a and b for the van der Waals equation. Search for “PRsolveZ” Two cases will be calls and you may wish to change the function name to “vdwsolveZ.” The third case of “PRsolveZ” will be the function that solves the cubic. Change the function name. Edit the formulas used for the cubic coefficients. Finally, specify a fluid and find the vapor pressure at the temperature corresponding to Tr = 0.7.
The definition of the acentric factor gives
ω = –log(Pr) – 1 = –log(0.20046) – 1 = –0.302
Comparing this value to the acentric factors listed in the table on the back flap, the only compound that even comes close is hydrogen, for which we rarely calculate fugacities at Tr < 1. This is the most significant shortcoming of the van der Waals equation. This shortcoming becomes most apparent when attempting to correlate phase-equilibria data for mixtures. Then it becomes very clear that accurate correlation of the mixture phase equilibria is impossible without accurate characterization of the pure component phase equilibria, and thus the van der Waals equation by itself is not useful for quantitative calculations. Correcting the repulsive contribution of the van der Waals equation using the Carnahan-Starling or ESD form gives significant improvement in the acentric factor. Another approach is to correct the attractive contribution in a way that cancels the error of the repulsive contribution. Cancellation is the approach that historically prevailed in the Redlich-Kwong, Soave, and Peng-Robinson equations.
The Equal Area Rule
As noted above, the swings in the P-V curve give rise to a cancellation in the area under the curve that becomes the free energy/fugacity. A brief discussion is helpful to develop an understanding of how the saturation pressure and liquid and vapor volumes are determined from such an isotherm.
To make this analysis quantitative, it is helpful to recall the formulas for the Gibbs departure functions, noting that the Gibbs departure is equal for the vapor and liquid phases (Eqn. 9.26).
In the final equation, the second term in the right-hand side braces represents the area under the isotherm, and the first term on the right-hand side represents the rectangular area described by drawing a horizontal line at the saturation pressure from the liquid volume to the vapor volume in Fig. 9.6(a). Since this area is subtracted from the total inside the braces, the shaded area above a vapor pressure is equal to the shaded area below the vapor pressure for each isotherm. This method of computing the saturation condition is very sensitive to the shape of the P-V curve in the vicinity of the critical point and can be quite useful in estimating saturation properties at near-critical conditions.
Although Eqn. 9.49 illustrates the concept of the equal area rule most clearly, it is not in the form that is most useful for practical application. Noting that GL = GV at equilibrium and rearranging Eqn. 9.47 gives
You should recognize the first term on the right-hand side as (AL – Aig)T,V – (AV – Aig)T,V. You probably have an expression for (A – Aig)T,V already derived. Solving for P is iterative in nature because we must first guess a value for P to solve for VV and VL. Five or six iterations normally suffice to converge to reasonable precision.5 The method is guaranteed to converge as long as a maximum and minimum exist in the P-V isotherm. Therefore, initiation begins with finding the extrema, a form of “stability check” (see below) in the sense that the absence of extrema indicates a single phase. The search for extrema is facilitated by noting that the vapor maximum must appear at ρ < ρc while the liquid minimum must appear at ρ < ρc. If Pmin > 0, then initialize to P0 = (Pmin+Pmax)/2. Otherwise, initialize by finding VV and VL at P = 10–12. The procedure is applied in Example 9.8 below.
Example 9.8. Vapor pressure using equal area rule
Convergence can be tricky near the critical point or at very low temperatures when using the equality of fugacity, as in Example 9.6. The equal area rule can be helpful in those situations. As an example, try calling the solver for CO2 at 30°C. Even though the initial guess from the shortcut equation is very good, the solver diverges. Alternatively, apply the equal area rule to solve as described above. Conditions in this range may be “critical” to designing CO2 refrigeration systems, so reliable convergence is important.
Solution
The first step is to construct an isotherm and find the spinodal densities and pressures. Fig. 9.10 shows that Pmin = 7.1917, Pmax = 7.2291, Vmax = 117.98, and Vmin = 94.509. Following the procedure above, P0 = 7.2104. Solving for the vapor and liquid roots at P0 in the usual way gives Vvap = 129.842, and Vliq = 88.160. Similarly, (AL – Aig)T,V = –1.0652 and (AV – Aig)T,V. = –0.7973, referring to the formula given in Example 8.6 on page 317:
Figure 9.10. Illustration of use of the equal area rule for a small van der Waals loop.
This leads to the next estimate of Psat as,
Psat = [–1.0652 + 0.7973 – ln(88.160/129.842)][8.314(303.15)/(129.842 – 88.160)] = 7.2129
Solving for the vapor and liquid roots and repeating twice more gives: Psat = 7.21288, shown below. Note the narrow range of pressures.
Leave a Reply