DENSITY MARCHING METHOD FOR DRAWING ... - ACS Publications

Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX-XXX

pubs.acs.org/IECR

Density Marching Method for Drawing Phase Envelopes. 3. P‑xy Diagrams of Binary Mixtures G. Venkatarathnam* Refrigeration and Air Conditioning Laboratory, Department of Mechanical Engineering, Indian Institute of Technology Madras, Chennai 600036, India S Supporting Information *

ABSTRACT: Some of the most widely used commercial process simulators and property programs cannot calculate Pxy diagrams of even simple binary mixtures satisfactorily when the specified temperature is greater than the critical temperature of one of the two components because of the limitations of the methods used by them currently. In this paper, we present a new method for drawing P-xy diagrams of binary mixtures that overcomes the limitations of the present methods. In the proposed method, the density of the vapor/liquid phase is varied systematically and the saturation point calculated at each of the specified density (instead of the conventional specified concentration) to generate P-xy diagrams. Two different approaches have been presented for the calculation of saturation points at specified density. The methods presented are an extension of the density marching method presented earlier for drawing two-phase and three-phase isopleths. A number of zeotropic and azeotropic mixture examples are presented to demonstrate the methods.

1. INTRODUCTION The study of phase equilibria of mixtures is important in the design of a variety of industrial and scientific applications. The estimation of the phase boundaries between different types of phase equilibria (vapor−liquid, liquid−liquid, liquid−liquid− vapor, as well as those involving the solid phases) in binary mixtures is the starting point of many phase equilibria studies. The calculation of phase boundaries at constant temperature/ pressure and representing them on a pressure−concentration (or temperature−concentration) plane is therefore a standard component of any thermodynamic property program or a process simulator. Most programs are capable of calculating the P-xy/T-xy diagrams of binary mixtures in the vapor−liquid region fairly well, including those of azeotropic mixtures. However, there are some exceptions. Figures 1, 2, and 3 show the P-xy diagram of a binary mixture of methane and ethane at a temperature of 250 K calculated using Aspen Plus,1 NIST Refprop,2 and Aspen Hysys3 respectively. The temperature has been chosen such that it is higher than the critical temperature of methane (190.56 K). The P-xy plots generated by Aspen Plus version 9 and NIST Refprop2 have been reproduced as is in Figures 1 and 2. Aspen Hysys version 9 gives an error and does not produce any plot, but does provide a table of the calculated values. Figure 3 has been drawn by us using the values in the table. The correct P-xy diagram produced with the methods developed in this work is shown in Figure 4 for comparison. It is evident from Figures 1−4 that even some of the most widely used commercial © XXXX American Chemical Society

process simulators and property programs cannot calculate P-xy diagrams out of the box in cases where the specified temperature is greater than the critical temperature of one of the two components of even simple binary mixtures. It is possible to overcome this situation by limiting the concentration steps to a value close to that of the critical point of the mixture, iteratively. Most commercial programs also have trouble producing the correct P-xy and T-xy diagrams in mixtures that exhibit vapor− liquid−liquid equilibria when the specified temperature is higher than the critical temperature of the low boiler. Figures 5 and 6 show the P-xy diagram of nitrogen-ethane at a temperature of 135 K (critical temperature of nitrogen2 = 126.19 K) generated using the Aspen Plus process simulator and the method described in this work, respectively. It is evident that the Aspen Plus process simulator is not able to reproduce the correct P-xy diagram correctly in this case also. Use of smaller concentration step size can overcome some of the problems in Figure 5. However, most current commercial programs are incapable of detecting the small region shown in the inset of Figure 6. It is evident from the above discussion that there is a need for new approaches for calculating P-xy diagrams in commercial Received: Revised: Accepted: Published: A

August 1, 2017 October 19, 2017 October 23, 2017 October 23, 2017 DOI: 10.1021/acs.iecr.7b03188 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 4. P-xy diagram of methane−ethane mixture at a temperature of 250 K determined using the density marching method presented in this paper with the Peng−Robinson equation of state.

Figure 1. P-xy diagram of methane−ethane binary mixture at a temperature of 250 K produced by the Aspen Plus process simulation software version 9.01 with the Peng−Robinson equation of state.

Figure 2. P-xy diagram of of methane−ethane mixture at a temperature of 250 K produced with the NIST Refprop program2 with the Helmholtz equation of state.

Figure 5. P-xy diagram of nitrogen-ethane estimated using the Aspen Plus process simulator, version 91 with the Peng−Robinson equation of state at a temperature of 135 K.

phase4 and three-phase isopleths.5 The main philosophy of this approach to draw space curves that loop back is to use a different independent variable (say α) to express the variables of the space curve (x and y) as functions of α (say x = f(α) and y = g(α)). The independent variable α is varied monotonically to generate different values of x and y, which when plotted yield space curves in the x−y plane. Venkatarathnam4,5 showed that the density of the vapor phase is a good independent variable for drawing isopleths. The methods are not repeated here; please refer to our earlier papers. In this paper, the density marching method is extended to the calculation of P-xy diagrams including three-phase lines. Different examples are presented to demonstrate the method. The method has immediate application in process simulators and thermodynamic property programs.

Figure 3. P-xy diagram of of methane−ethane mixture at a temperature of 250 K drawn using the data estimated by the Aspen Hysys program version 93 with the Soave−Redlich−Kwong equation of state.

process simulators and property programs. A new method that does not require any heuristics, known as the density marching method was presented by the author earlier to estimate twoB

DOI: 10.1021/acs.iecr.7b03188 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


energy analysis of Baker et al.8 to identify the locus of unstable points and the vapor−liquid−liquid coexistence temperature/ pressure. Michelsen and Møllerup9 present a slightly different method for drawing the phase diagrams of binary mixtures. They suggest the solution of the following equations: ln K i + ln ϕi v − ln ϕi l = 0

yi − K ixi = 0

i = 1,2

(4)

i = 1,2

(5)

x1 + x 2 − 1 = 0

(6)

y1 + y2 − 1 = 0

(7)

ln Ki, P, T, yi (or xi) are the variables. Usually one of ln Ki, T, or P is specified. Michelsen and Møllerup9 showed that the maxima and minima seen in the bubble line of a P-xy diagram similar to that shown in Figure 7 is due to one of the phases being unstable

Figure 6. P-xy diagram of nitrogen−ethane calculated using the density marching method of this work using the Helmholtz energy equation of state of the NIST Refprop program2 at a temperature of 135 K.

2. LITERATURE A number of authors have studied the methods for drawing Pxy and T-xy diagrams. Li and Ngheim6 proposed a general purpose method for phase envelopes including the construction of P-xy and T-xy diagrams. The following equations are solved simultaneously in the case of P-xy and T-xy diagrams at specified temperature (P-xy) or specified pressure (T-xy), and a specified vapor mole fraction v/f using the Newton−Raphson method. ln K i + ln ϕi v − ln ϕi l = 0 N

∑ i=1

zi(K i − 1) =0 v 1 + f (K i − 1)

⎛v⎞ v −⎜ ⎟ =0 f ⎝ f ⎠specified

i = 1 ... N

(1) Figure 7. P-xy diagram of refrigerants R14 (tetrafluoromethane) and R23 (trifluoromethane) estimated with the NIST Refprop program2 using the Helmholtz energy equation of state at a temperature of 146 K.

(2)

indicating the occurrence of three-phase (vapor−liquid−liquid) equilibria. They presented two different methods to determine the three-phase line from the unstable P-xy diagram shown in Figure 7. Cismondi and Michelsen10 use a slightly different approach for calculating the phase boundaries of binary fluids. They solve the following equations simultaneously:

(3)

In the above equations, v/f is the vapor mole fraction, zi is the composition, and Ki the distribution coefficient of the ith component of the feed. The design variables are Ki, either T or P. In their method, Li and Ngheim6 use the pressure/ temperature (P-xy/T-xy) as the initial specification variable. The specification variable is changed at each step as the phase boundaries are generated, to any one of the calculation variables. The step size of the specification variable is so chosen that the calculations converge in 3 iterations. If the number of iterations required for the Newton−Raphson method to solve eqs 1−3 is greater than 4, the calculation is repeated with a smaller step size of the specified variable, or using the specified variable from a previous step. The specification variable for the next step is chosen to be the most sensitive variable of the previous calculation using the Jacobian of eqs 1−3 used in the previous iteration. The above method is essentially a modification of the method for drawing saturation envelopes by Michelsen7 to P-xy and T-xy diagrams. Li and Ngheim6 provided examples of zeotrope and azeotrope binary mixtures. The solution of eqs 1−3 proposed by Li and Ngheim6 cannot predict the three phase pressure (temperature) when two immiscible liquid phases are present. Li and Ngheim6 therefore subject the phase boundaries to the Gibbs

ln P l(X , T , v l) − ln P(Y , T , v v) = 0

(8)

ln f1l (X , T , v l) − ln f1v (Y , T , v v) = 0

(9)

ln f 2l (X , T , v l) − ln f 2v (Y , T , v v) = 0

(10)

g (X ) − specified value = 0

(11) v

v

The variables of the calculation: lnx1, lny2, lnv , lnv . g(X) can either be pressure (T-xy diagram) or temperature (P-xy diagram), or any of the variables of the calculation, or their combinations (viz., yi − xi, vv/vl etc.). The main advantage in using molar volume (vl, vv) is that the equation of state need not be solved for density as in the case of eq 1 particularly advantageous in complex multiparameter equations of state such as the Helmholtz energy equation of state. The use of natural log of volume or concentration ensures that volume/ concentration never becomes negative. With a bad starting value, the concentration of the first component of the x-phase C


Article


3. CALCULATION OF DEW AND BUBBLE POINTS In commercial programs such as the Aspen Plus and the Aspen Hysys, the P-xy diagram is generated by calculating the bubble point pressure at specified liquid concentration. The concentration of the dew point is obtained from the K-value. An alternate approach used by the NIST Refprop program involves calculation of both the bubble and dew point pressures separately at specified concentration and temperature. The concentration of the low boiler is normally varied in equal steps from 0 to 1 to generate the P-xy diagram in this case. Conventional general purpose multicomponent dew/bubble point routines can therefore be used straight away in the above two approaches. It is also possible to use a different approach for calculating bubble/dew points and also draw the P-xy diagram. In the present work, the dew/bubble point pressure is calculated at a specified vapor phase density instead of specified vapor or liquid phase composition. Both successive substitution and Newton−Raphson based methods can be employed. Because successive substitution methods converge linearly, Newton− Raphson methods are used after a few successive substitution iterations (typically 4 in our program) to speed up convergence. If good guess values are available after a few density steps (three or four steps in our program), a second- or third-order polynomial is fit based on the last few points, and the same is used to extrapolate the values of variables w.r.t. density of the vapor/liquid phase to generate the guess values of the next iteration using Taylor series expansion. The density of the vapor phase is increased in steps to draw the entire P-xy diagram. More details of the procedure for drawing the entire phase envelope are presented in the next section. The first step is to estimate the dew/bubble point temperature at the specified density using successive substitution method or a Newton−Raphson method. 3.1. Successive Substitution Method for Determining Saturation Pressure/Temperature of a Binary Mixture at Specified Density of Vapor Phase. Most process simulators and thermodynamic programs use the traditional dew and bubble point temperature/pressure algorithms written for multicomponent mixtures to determine the saturation pressure/temperature. These programs require the concentration of the feed, and either pressure or temperature as the input, depending on T-xy or P-xy calculation. None of them can use temperature/pressure and density of vapor phase as input to determine the pressure/temperature (depending on P-xy or T-xy calculation) and the concentration of the vapor and liquid phases as the output. The following successive substitution method can be used to determine the dew point pressure at a specified density of the vapor phase, temperature, and concentration: 1. Specify the density of the vapor phase at which the dew/ bubble point pressure is required. 2. Assume a K-value for all the components (say ideal Ki) 3. Assume the concentration of the vapor phase (y1, y2). 4. Calculate the dew point pressure from the guess concentration of vapor phase (Y), the density of the vapor phase and the specified temperature from the equation of state:

or the second component of the y-phase can, however, become greater than one. As in the case of Li and Ngheim6 and Michelsen,7 the extrapolation is carried out w.r.t the most sensitive variable. The authors state that the step size is varied such that iterations converge in 3 to 4 iterations. If larger number of iterations is required, the step size of the sensitive variable is reduced, or if the convergence occurs in 2 or 3 iterations, the step size is increased by a known fraction. Though being good methods, the main deficiency of the above methods is in the change of the specification variable along the length of the phase boundaries, and the choice of the step size based on the number of iterations required for converging any iteration. Consider the variation of the concentration of methane in the vapor and liquid phases shown in Figure 4 for a simple mixture of methane and ethane. Extrapolation of the liquid concentration w.r.t vapor concentration can therefore be problematic in the retrograde condensation region, particularly close to saddle points (∂P/ ∂y = 0). This problem is normally overcome by (a) redoing the calculations with very small steps, (b) using the same specification variable as in the previous step, (c) choosing the most appropriate variable w.r.t to which different variables are extrapolated (to serve as the guess values for the next step), and a host of other heuristics not fully described in the open literature. Consequently, not all commercial programs are able to replicate to the same degree the success of the authors of the original methods in drawing phase envelopes at all conditions, sometimes even after many decades of the publication of the original articles. Cismondi and Michelsen10 presented a novel alternate approach to that presented by Møllerup and Michelsen.9 In their approach, they first calculate all the possible critical lines (VL, LL, high pressure FF etc.) and also identify the type of mixture according to the classification of Konynenburg and Scott.11 They use the critical point data as the starting point to construct the different phase boundaries. This is a very novel and probably the most appropriate way of calculating phase boundaries of binary mixtures. This approach, however, requires routines for calculating critical points, and critical end points using different equations of state. Aspen Plus, a very widely used process simulator, as of version 9, does not have the ability to calculate even a simple vapor−liquid critical point. The need to calculate not only vapor−liquid critical lines but also critical end points, and high pressure (fluid−fluid) critical lines before a simple P-xy diagram, can be calculated makes this method not very suitable for use in commercial property programs and process simulators. Further, the algorithms in commercial process simulators are expected to work with activity coefficient methods for the liquid phase and an equation of state for the vapor phase. Because a critical point cannot be predicted in this case, the method of Cismondi and Michelsen10 that relies on the critical point to initiate the calculations is less suitable for use in process simulators. Very recently, Bell and Deiters12 presented a new approach for drawing P-xy and T-xy diagrams. In their approach, a system of differential equations are solved simultaneously to generate the P-xy and T-xy diagrams. They use an isochoric approach originally used by Quiñ ones-Cisneros and Deiters13 for isopleths in their work. The P-xy and T-xy diagrams were calculated by varying either concentration or pressure (temperature). They indicate that more exotic algorithms are required in the case of phase envelopes that exhibit multiple minima or maxima in the tracing variable.

P = F1(Y , ρ v , T ) D

(12) DOI: 10.1021/acs.iecr.7b03188 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research 5. Calculate the fugacity coefficient of the vapor phase (ln ϕvi ) from the composition of the vapor phase, density, and temperature using the equation of state v

v

ln ϕi = F2(Y , ρ , T )

(13)

6. Calculate the composition of the liquid phase from that of the given vapor phase and the assumed K-value Ki y1 K1

+

y2 K2

(14)

(19)

P(X , T , ρ l ) = P(X , T , ρ v )

(20)

v ln ρ v = ln ρspecified

(21a)

l ln ρ l = ln ρspecified

(21b)

And another specification equation, chosen from the following to accommodate P-xy/T-xy diagram specification: Saturation pressure calculation

(15)

ln T = ln Tspecified

8. Calculate the new K-values from the fugacity coefficients ln K inew = ln ϕi l − ln ϕi v

ln f 2l = ln f 2v

Or

7. The fugacity coefficient of the liquid phase (ln ϕli) is then determined using the specified temperature, pressure (calculated in step 5) and the normalized liquid phase composition calculated in step 4 using an equation of state or an activity coefficient method. ln ϕil = F2(X , P , T )

(18)

Along with one specification equation, chosen from the following:

yi

xi =

ln f1l = ln f1v

(22a)

Or (16)

Saturation temperature calculation

9. If ||ln Ki − ln Knew i ||2 < ϵ, the calculations are assumed to have converged, else the Knew value is used as the Ki value i in the subsequent iteration. A value of 10−6 has been used for ϵ in our calculation. 10. Because x1 + x2 = 1, calculate the new vapor phase composition of the low boiler (y1) using the following expression: y1 1 − y1 =1 new + K1 K 2new (17)

P(Y , T , ρ v ) = Pspecified

(22b)

The calculation variables are ln y1, ln x2, ln ρv, ln ρl, and ln T. The same routine can also be used to determine the saturation pressure/temperature at specified concentration by using other specification equations such as ln y1 = ln y1,specified or ln x1 = ln x1,specified. The Jacobian of eqs 18−22a has been presented in Supporting Information (Appendix-A). The above approach is useful only when an equation of state is used for both vapor and liquid phases. A slightly different approach is useful when either an activity coefficient or an equation of state method is used for the liquid phase. In this approach, the following equations are solved simultaneously:

Calculate the composition of the second component using y2 = 1 − y1 . 11. If the total iterations is less than a few iterations (say 4− 6), go to step 4, else exit to the Newton−Raphson method. The above method is used in our program to calculate the saturation temperature or pressure for the first two or three points before switching to the full Newton−Raphson based calculation. The method can also be easily modified to accommodate a specified liquid phase density in place of specified vapor phase density. We have found the method to be particularly useful in drawing the phase boundaries in the high pressure liquid−liquid region because the change in pressure between steps is very large. The first few successive substitution calculations normally results in very stable convergence of the subsequent Newton−Raphson iterations in the liquid−liquid phase boundary calculations. The density of the lighter liquid phase in the previous iteration is always used as the specification variable in the case of a liquid−liquid boundary. Only a minor modification of step 4 needs to be made to determine the temperature at specified pressure to determine the saturation temperature. 3.2. Newton−Raphson Method for Determining the Saturation Pressure/Temperature of a Binary Mixture at Specified Density of Vapor/Liquid Phase. The successive substitution method described in the previous section converges linearly, and can sometimes be slow. The convergence can be speeded up by using the Newton− Raphson method to solve the governing equations. It is convenient to use the following equations:

ln K1 + ln ϕ1v − ln ϕ1l = 0

(23)

ln K 2 + ln ϕ2v − ln ϕ2l = 0

(24)

2

∑ i=1 2

∑ i=1

K iξi 1 + (K i − 1)

() v f

ξi 1 + (K i − 1)

() v f

−1=0 (25)

−1=0 (26)

Along with one specification equation, chosen from the following: v ln ρ v = ln ρspecified

(27a)

l ln ρ l = ln ρspecified

(27b)

Or

And another specification equation, chosen from the following to accommodate P-xy/T-xy diagram specification: Saturation pressure calculation ln T = ln Tspecified

(28a)

Or E


Article


even with the guess concentration is not close to the final result, and makes the method most suitable for use in process simulators.

Saturation temperature calculation v

P(Y , T , ρ ) = Pspecified

(28b)

The following are the variables of the above equations: ln K1, ln K2, ln ξ1, ln ξ2, and ln T. The concentration of the vapor and liquid phases (yi, xi) are related to the variable ξ as follows:

xi′ =

ξi v 1 + (K i − 1) f

(29)

yi′ = K ixi′ yi =

xi =

4. CALCULATION OF THE P-XY DIAGRAM BY THE DENSITY MARCHING METHOD The density marching method consists essentially of increasing the specified density of the vapor phase by a constant value and solving for the saturation pressure or saturation temperature. Calculations are always started at zero concentration of the low boiler, or the fluid with the lowest saturation pressure at the specified temperature. The method can be used with both zeotropes and azeotropes. The complete method is described below: 1. Check for a Bancroft point. Almost all mixtures with a Bancroft point exhibit azeotropic behavior.14 However, all azeotropes do not have a Bancroft point. If a Bancroft point is not found, check for the occurrence of azeotropic composition by calculating the dew point of the mixture at y1 = 0.01 using conventional methods at the normal boiling point temperature and at a pressure of 1 atm. If the saturation pressure at y1 = 0.01 is less than that at y1 = 0 then a negative azeotrope is predicted. Similarly if the saturation temperature at y1 = 0.01 is less than that at y1 = 0, then a positive azeotrope is predicted. The liquid phase density is chosen as the independent (specification) variable when an azeotrope is predicted and the vapor phase density is used when a zeotrope is predicted. 2. Increase the density of the vapor/liquid phase in small steps, say min (0.05 mol/L,(ρpure − ρpure 1 2 )/nsteps), with nsteps being the minimum number of steps desired. When ρv ≪ 1, choose ρv = ρv(previous step) × 1.1. 3. The direction of change of density in steps 1−2 is decided by the difference in the density of the vapor/ liquid phase of the high boiler and that at y1 = 0.01. 4. Estimate the dew/bubble point pressure/temperature using the successive substitution in the first three steps followed by one of the two Newton−Raphson methods described in the previous section. Extrapolate the different variables ρv, ρl, ln xi, ln yi, ln ξi etc. with respect to the density of the phase that is changed using the values converged in the last three or four steps using Taylor series expansion. Use direct Newton−Raphson calculations after the first few steps. 5. Saturation pressure or temperature calculations normally do not fail when reasonable step size (say 0.1 mol/L) is chosen, and the guess values of the next step are estimated from those of previous 3 or 4 steps. Any failure is normally due to either very large steps or the occurrence of a different phase equilibrium (say V-L1L2, V-L2 equilibrium etc.). Whenever saturation point calculations do not converge, the following are checked for: (a) Critical point (|ln K1| < ϵ, or ln K1 (this step) × ln K1 (previous step) < 0, and |ρl − ρv| < ϵ). It is important to check the absolute value of the density difference because the liquid and vapor phases switch when the critical point is crossed. (b) Azeotropic point |ln K1| < ϵ, or ln K1 (this step) × ln K1 (previous step) < 0, and |ρl − ρv| > ϵ)

(30)

yi′ y1′ + y2′

(31)

xi′ x1′ + x 2′

(32)

In the above expression, v/f is the specified vapor fraction. At convergence, ξi = yi for a dew point calculation and ξi = xi for a bubble point calculation. The use of ln ξi as the variable ensures that ξi is always positive. There is no upper limit on the value of ξi. This makes the Newton−Raphson iterations quite stable, and in most cases failure occurs only when Ki → 1 (critical point, azeotropic point). The Jacobian is slightly different from the usual. For example, ⎛ ∂ln ϕ v (P , T , n ̅ v ) ⎞ i ⎟ ⎜ ∂ln ξ1 ⎠ln K , ξ ⎝ i

2

⎡ ⎛ ⎛ ∂ln ϕi v ⎞ ⎤ ∂ln ϕi v ⎞ = y1y2 ⎢n v ⎜ ⎟ − nv⎜ ⎟ ⎥ ⎢ ⎝ ∂n1 ⎠ ⎥ n ∂ ⎝ 2 ⎠P , T ⎦ ⎣ P ,T ⎛ ∂ln ϕi v ⎞ ⎛ ∂P ⎞ +⎜ ⎟ ⎜ ⎟ ⎝ ∂P ⎠T , n v ⎝ ∂ln ξ1 ⎠ln K , ξ i

2

(33)

The last term on the right-hand side is evaluated from the variables of the phase whose density is specified. For example, when the density of the vapor phase is specified, ⎛ ∂P ⎞ ⎛ ∂P ⎞ = −⎜ ⎟ ⎟ ⎜ ⎝ ∂ln ξ2 ⎠ln K , ξ ⎝ ∂ln ξ1 ⎠ln K , ξ i

2

i 1

⎡ ⎛ ∂P ⎞ ⎛ ∂P ⎞⎤ = y1y2 ⎢n v ⎜ v ⎟ − n v ⎜ v ⎟⎥ ⎢⎣ ⎝ ∂n1 ⎠ ⎝ ∂n2 ⎠⎥⎦

(34)

Similarly, when the density of the liquid phase is specified, ⎛ ∂P ⎞ ⎛ ∂P ⎞ = −⎜ ⎟ ⎟ ⎜ ⎝ ∂ln ξ2 ⎠ln K , ξ ⎝ ∂ln ξ1 ⎠ln K , ξ i

2

i 1

⎡ ⎛ ∂P ⎞ ⎛ ∂P ⎞⎤ = x1x 2⎢n l ⎜⎜ l ⎟⎟ − n l ⎜⎜ l ⎟⎟⎥ ⎢⎣ ⎝ ∂n1 ⎠ ⎝ ∂n2 ⎠⎥⎦

(35)

The expressions for the complete Jacobian matrix are presented in Supporting Information (Appendix-B). The main advantage in using lnξi instead of the composition of the vapor (yi) or the liquid phase (xi) is that the sum (ξ1 + ξ2) can be greater than one during the iterations. Convergence is therefore guaranteed F


Article


variables. The density marching method presented in Section 4 is equally applicable to the isochoric approach of Deiters and his co-workers12,13 for determining the bubble/dew point.

(c) Three-phase line The saturation pressure/temperature calculation methods presented in Sections 3.1 and 3.2 do not check for the occurrence of the coexistence of two liquid phases. Hence the calculations can sometimes fail close to the three phase bubble point pressure/temperature, particularly in mixtures that exhibit lower critical solution temperature (LCST) behavior (L1 = L2 − V). In mixtures that exhibit a K-point (L1 − L2 = V), the calculations can cross the three-phase pressure/temperature, but the predicted phases can turn out to be unstable. The calculations can proceed to trace the liquid− liquid boundary, as shown in Figure 8. The

5. RESULTS AND DISCUSSION Figure 9 shows the variation of the density of the liquid phase (ρl), ln K value of methane and the concentration of methane

Figure 9. Variation of the density of the liquid phase (ρl), ln K value of methane and the concentration of methane in the vapor phase (yi) for a mixture of methane and ethane at a temperature of 250 K estimated with the Peng−Robinson equation of state. The P-xy diagram of this mixture is shown in Figure 4.

in the vapor phase (yi) with the density of the vapor phase (ρv) for a mixture of methane and ethane at a temperature of 250 K estimated with the Peng−Robinson equation of state. It can be seen that all the variables except concentration y vary monotonically with the density of the vapor phase even though the pressure−concentration diagram shown in Figure 4 loops back. The variation of concentration with density of the vapor phase is also unimodal. Hence the density of the vapor phase is the right marching (scan) variable for drawing P-xy diagram of the methane−ethane mixture. The monotonic variation of the other variables also allows easy extrapolation using Taylor series expansion. Even if conventional methods are used for drawing P-xy diagram, it is beneficial to guess the variables of the next saturation point calculation by extrapolating them w.r.t the density of the vapor phase, as in the case of isopleths presented in our earlier works on drawing two-phase4 and three-phase isopleths.5 In Figure 9, the density of the vapor phase was varied in steps of 0.1 mol/L. The calculations easily pass the critical point and the entire curve is traced if the calculations are allowed to be carried further. Unlike isopleths, there is no need to trace the curve beyond the critical point. Conditions 5(a) of Section 4 are satisfied when the critical point is crossed. The dew point calculations (v/f = 1) are therefore terminated, and a critical point calculation performed. Figure 10 shows the P-xy diagram of a binary mixture of refrigerants R14 and R23 calculated with the Helmholtz energy equation of state of the NIST Refprop program2 using the density marching method presented in this work at a temperature of 146 K. The calculations are done in the following stages: (a) Calculate the three-phase bubble point pressure first. (b) Increase the density of the vapor phase in steps until the pressure is lower than the three-phase bubble point

Figure 8. P-xy diagram of nitrogen-ethane mixture estimated with the Helmholtz energy equation of state of the NIST Refprop program2 at a temperature of 125 K using the density marching approach of this work when the occurrence of VLLE is not considered.

pressure at the intersection point is the threephase pressure.Two different approaches can be used to prevent the calculations converging to unstable condition shown in Figure 8: (i) Perform a three-phase bubble point calculation right at the beginning or (ii) when the density of the specified phase is greater than that of pure fluids perform a three-phase calculation. Our preference is to do a three-phase bubble point calculation right at the beginning since in most cases because the total time taken to estimate the three-phase bubble point, critical/azeotropic point, as well as calculation of a full P-xy diagrams is of the order of a 0.01−0.03 s only for most mixtures on an Intel i5 based personal computer used in this work. Occasionally, failure may occur due to a problem with the equation of state model (including wrong binary interaction parameters, wrong mixing rule, hydrogen bonds, etc.), discontinuous P-xy/T-xy diagram observed in high pressure fluid−fluid equilibria conditions, etc. One can switch to conventional methods by changing the specification equation in eq 28a to ln y1 = ln y1,specified in all such cases. In our experience, these cases are quite rare, and not of much importance in process simulators, and industrial practice. Quiñones-Cisneros and Deiters13 use a different approach for determining the phase equilibria. This approach, known as the Isochoric approach, uses the product of density and concentration of each component of the two phases as the G


Article


Figure 11. P-xy diagram of refrigerants R22 and R115 at a temperature of 250 K estimated with the Helmholtz energy equation of state of state of the NIST Refprop program.2

Figure 10. P-xy diagram of a binary mixture of refrigerants R14 and R23 calculated with the Helmholtz energy equation of state of the NIST Refprop program2 using the density marching method presented in this work at a temperature of 146 K.

pressure, or the density of the vapor phase is less than that at the three-phase pressure to essentially calculate the boundaries of the vapor and liquid-1 (V-L1) phases. (c) Restart the density marching calculations at the second liquid phase concentration and the density of vapor phase at the three-phase pressure to calculate the boundaries of the vapor and second liquid (V-L2) phases. Calculations are stopped when the density of the vapor phase equals the density of pure low boiler at the specified temperature. (d) Restart the density marching calculations at the first liquid phase concentration and the density of the first liquid phase at the three-phase pressure to calculate the boundaries of the two liquid (L1-L2) phases. If the density of the lighter liquid becomes more than that of the heavier liquid at any point, switch to marching of the lighter liquid phase. This step is important to prevent large pressure steps close to the L1-L2 critical point. This switching, should, however be avoided while crossing the critical point in vapor−liquid boundaries. In our program, the two cases are distinguished using the phase identification method of Venkatarathnam and Oellrich.15 In cases where the specified temperature is greater than that of the critical temperature of the low boiler (Figure 6), the calculation of the vapor and second liquid phase (V-L2) boundary is stopped at the critical temperature of the mixture. When a vapor−liquid critical point is encountered while drawing a P-xy diagram at temperatures less than the critical temperatures of both fluids, then the calculations need to be restarted from both ends as shown in ref.16 as the critical pressure of the mixture is lower than that of either pure fluids. Figure 11 shows the P-xy diagram of an azeotropic mixture of refrigerants R22 and R115 at a temperature of 250 K estimated with the Helmholtz energy equation of state of state of the NIST Refprop program.2 Figure 12 shows the variation of the density of the vapor phase and pressure w.r.t. the density of the liquid phase. It can be seen that the density of the vapor phase varies unimodally with the density of the liquid phase. The density of the liquid phase should therefore be varied systematically in the case of positive azeotropes. Figure 13

Figure 12. Variation of density of vapor phase and pressure with the density of the liquid phase for the azeotropic mixture of refrigerants R22 and R115 shown in Figure.x.

Figure 13. Variation of the vapor and liquid phases with the density of the liquid phase for the azeotropic mixture of refrigerants R22 and R115 shown in Figure 11.

shows the variation of the mole fraction of R22 with the density of the liquid phase. It can be seen that the mole fraction of R22 H


Article


Figure 6 but also includes the liquid−liquid equilibria phase boundaries. The inset shown in Figure 6 is not shown in Figure 15. Figure 16 shows the variation of concentration of the two

varies monotonically with that of the density of the liquid phase. The calculations pass the azeotropic composition where ln Ki = 0 without any difficulty when different values of liquid phase density are specified to generate the complete phase envelope. Similarly, all the variables of calculation are extrapolated w.r.t. to the density of the liquid phase in the case of azeotropes. This approach can be used in traditional Pxy diagram calculations as well. Figure 14 shows the P-xy diagram of a negative azeotropic mixture of refrigerants R134a and R152a at a temperature of

Figure 16. Variation of mole fraction of nitrogen of the two liquid phases on the liquid−liquid equilibria boundaries for a mixture of nitrogen and ethane (see Figure 15).

liquid phases with the density of the two phases. It can be seen that the concentration of the L1 phase loops backs when plotted against density of that phase. Thus, L2 phase is the correct phase for density marching. However, this information is not known apriori. Calculations are normally performed with the systematic increase of the density of the lighter liquid (L1) phase from the three-phase line. Figure 17 shows the variation

Figure 14. P-xy diagram of an azeotropic mixture of refrigerants R134a and R152a at a temperature of 250 K estimated with the Helmholtz energy equation of state of the NIST Refprop program.2 The complete diagram is generated by varying the density of the liquid phase systematically.

250 K estimated with the Helmholtz energy equation of state of the NIST Refprop program.2 The complete diagram is generated by varying the density of the liquid phase systematically from that of pure R152a to that of pure R134a. Even though the bubble and dew lines are very close, the calculations are completed using the methods presented in this work very easily. Figure 15 shows the P-xy diagram of a mixture of nitrogen and ethane at a temperature of 135 K. The figure is the same as

Figure 17. Variation of molar density of the two liquid phases with mole fraction of nitrogen in the L1 liquid phase for a mixture of nitrogen and ethane shown in Figure 15.

of density of the two phases with the concentration of nitrogen in the L1 liquid phase. Calculations start at a concentration of 95.6 mol % and proceed toward the critical point. The concentration of nitrogen in the L1 liquid phase decreases with the calculation steps. Thus, the two curves are calculated starting from the right-hand side. It can be observed that at the density of the L1 phase becomes more than that of the L2 phase when nitrogen concentration in the L1 phase is less than 84 mol %. The marching variable is now changed from density of L1 liquid phase to that of the L2 liquid phase, and the rest of

Figure 15. P-xy diagram of nitrogen−ethane calculated using the density marching method of this work using the Helmholtz energy equation of state of the NIST Refprop program2 at a temperature of 135 K showing both vapor−liquid and liquid−liquid equilibria regions. I


Article

Industrial & Engineering Chemistry Research the calculations performed until the critical point is reached. The looping back of the concentration-density line in Figure 16 is thus not a problem even if calculations are started at specified values of L1 phase density. Figures 14−16 also demonstrate the possibility of drawing the total phase diagram of nitrogen− ethane mixture with the methods proposed in this work. 5.1. Special Cases. The methods presented in this work assume that the saturation point of at least one pure fluid lies on the phase envelope. Thus, all calculations begin at the saturation point of a pure fluid. In some cases, the P-xy diagram is not connected to either pure fluid saturation point. This occurs in the case of mixtures that exhibit gas−gas equilibria. Michelsen and Møllerup9 have presented approaches for these cases, as well as similar cases that occur during the calculation of T-xy diagrams. Another assumption that is made is that only one critical point occurs at the specified temperature. Two critical points can occur at the same temperature as in the case of ethane− carbon dioxide mixture at a temperature of 290 K as shown by Michelsen and Møllerup9 as well as Bell and Deiters.16 In this particular case, the critical pressure of mixtures is lower than that of the two pure fluids. In such cases, the calculations are restarted at the opposite end to complete the diagram. Because this particular case has been well discussed in the literature, the same is not repeated again here. The method presented in this work also fails in some azeotropic mixtures in which the density of the liquid phase does not vary monotonically with concentration (say carbon dioxide−ethane mixture at temperature of 285 K). Calculations terminate at the azeotropic composition in such cases. Calculations can be restarted from the opposite end with either the density of the vapor or the liquid phase as the marching variable in all such cases to complete the diagram, as in the case of carbon dioxide−ethane mixtures with two critical points discussed above The present method has not been tested for all the special cases that can occur. No phase equilibria calculation method can be guaranteed to work universally (see for example, the difficulty of determining phase stability in the case of a mixture of methane and hydrogen sulfide in Michelsen17), including the present work. Improvements to phase equilibria calculation methods for special cases, and the development of new methods that overcome existing limitations form an important part of current research on phase equilibria. Finally, it should be noted that the methods presented in this work need to be modified to generate T-xy diagrams because molar density is inversely proportional to temperature. There are other problems as well such as the nonlinear variation of density (specific volume) of the mixture with concentration at specified pressure in some cases. The development of density marching methods suitable for drawing general purpose T-xy diagrams is a topic of continuing research.

■

1. The vapor phase density is the most appropriate independent marching (scan) variable for drawing P-xy diagrams of zeotropic binary mixtures and the liquid phase density that of azeotropic mixtures. 2. The use of an independent variable ln ξi in saturation temperature/pressure calculation ensures good convergence even when ideal K values are used in the first few steps of the P-xy/T-xy calculations. ξi converges to yi or xi (corresponding to v/f = 1 or 0, respectively).

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.7b03188. The expressions for calculating the Jacobian of eqs 18−22a, and eqs 23−28a (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*G. Venkatarathnam. E-mail: [email protected]. ORCID

G. Venkatarathnam: 0000-0003-3470-886X Notes

The author declares no competing financial interest.

■

■ ■

6. CONCLUSIONS A new approach for calculating both simple and complex P-xy diagrams has been presented in this work. A new successive substitution method and a new Newton−Raphson based method have also been presented for estimating saturation points of binary mixtures. Many examples have been presented to demonstrate the methods developed. The methods presented have immediate application in process simulators and thermodynamic programs. The following conclusions can be drawn from this work:

■

NOMENCLATURE f i = Fugacity of the ith component (MPa) Ki = K-value of the ith component n = Number of moles ni = Number of moles of the ith component P = Pressure (MPa) T = Temperature (K) v/f = Vapor fraction (mol/mol) v = Specific molar volume (l/mol) xi = Mole fraction of the ith component of the liquid phase X = Array containing the mole fraction of the two components in the liquid phase yi = Mole fraction of the ith component of the vapor phase Y = Array containing the mole fraction of the two components in the vapor phase zi = Mole fraction of the ith component of the feed GREEK LETTERS ϕi = Fugacity coefficient of the ith component ρ = Molar density (mol/L) ξi = Independent variable in eqs 23 and 26 SUPERSCRIPTS l = Liquid phase new = New calculation pure = Pure fluid v = Vapor phase REFERENCES

(1) Aspen Technologies. Aspen Plus Process Simulator, Version 9.0; Aspen Technologies: Bedford, MA, 2016. (2) Lemmon, E. W.; Huber, M. L.; McLinden, M. O. Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 9.1; NIST Standard Reference Database 23, National Institute of Standards and Technology: Gaithersburg, MD, 2013. (3) Aspen Technologies. Aspen Hysys Simulation Program, Version 9.0; Aspen Technologies: Bedford, MA, 2016. J


Article

Industrial & Engineering Chemistry Research (4) Venkatarathnam, G. Density Marching Method for Calculating Phase Envelopes. Ind. Eng. Chem. Res. 2014, 53 (9), 3723. (5) Venkatarathnam, G. Density Marching Method for Calculating Phase Envelopes. 2. Three-Phase Envelopes. Ind. Eng. Chem. Res. 2014, 53 (30), 12122. (6) Li, Y.-K.; Nghiem, L. X. The Development of a General Phase Envelope Construction Algorithm for Reservoir Fluid Studies. In Proc. of the 57th Annual Fall Technical Conference and Exhibition of the Society of Petroleum Engineers of AIME, New Orleans, LA, Sept. 26−29, 1982; Paper No. SPE 11198, pp 1−16. (7) Michelsen, M. Saturation Point Calculations. Fluid Phase Equilib. 1985, 23, 181. (8) Baker, L. E.; Luks, K. D.; Pierce, A. C. Gibbs Energy Analysis of Phase Equilibria. SPEJ, Soc. Pet. Eng. J. 1982, 22 (5), 731. (9) Michelsen, M. L.; Møllerup, J. M. Thermodynamic Models: Fundamentals & Computational Aspects; Tie-Line Publications: Holte, Denmark, 2004. (10) Cismondi, M.; Michelsen, M. Automated Calculation of Complete Pxy and Txy Diagrams for Binary Systems. Fluid Phase Equilib. 2007, 259 (2), 228. (11) Konynenburg, P. H. V.; Scott, R. L. Critical Lines and Phase Equilibria in Binary Van Der Waals Mixtures. Philos. Trans. R. Soc., A 1980, 298 (1442), 495. (12) Bell, I.; Deiters, U. Isochoric Thermodynamics and Binary Mixture P-X and T-X Diagrams. In 14th Joint European Thermodynamics Conference, Budapest, Hungary, May 21−25, 2017; pp 4−6. (13) Quiñones-Cisneros, S. E.; Deiters, U. K. An Efficient Algorithm for the Calculation of Phase Envelopes of Fluid Mixtures. Fluid Phase Equilib. 2012, 329, 22. (14) Morrison, G.; McLinden, M. O. Azeotropy in Refrigerant Mixtures. Int. J. Refrig. 1993, 16 (2), 129. (15) Venkatarathnam, G.; Oellrich, L. R. Identification of the Phase of a Fluid Using Partial Derivatives of Pressure, Volume, and Temperature without Reference to Saturation Properties: Applications in Phase Equilibria Calculations. Fluid Phase Equilib. 2011, 301 (2), 225. (16) Deiters, U. K. Calculation of Equilibria between Fluid and Solid Phases in Binary Mixtures at High Pressures from Equations of State. Fluid Phase Equilib. 1985, 20, 275. (17) Michelsen, M. L. The Isothermal Flash Problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, 1.

K


DENSITY MARCHING METHOD FOR DRAWING ... - ACS Publications

Recommend Documents