Density Marching Method for Calculating Phase Envelopes. 2. Three

Article pubs.acs.org/IECR

Density Marching Method for Calculating Phase Envelopes. 2. ThreePhase Envelopes G. Venkatarathnam* Refrigeration and Air Conditioning Laboratory, Department of Mechanical Engineering, Indian Institute of Technology Madras, Chennai 600036, India ABSTRACT: There are very few automatic procedures in the literature to calculate three-phase envelopes (VLLE boundaries). The best existing procedures require selecting each step manually to complete the three-phase envelope. In this article, a density marching procedure and a saturation point calculation at specified density of one of the incipient phases developed recently for drawing two-phase envelopes is extended to three-phase envelopes. Procedures have been presented to determine three-phase bubble and dew points reliably. Two examples are presented to demonstrate the method. The method has direct application in process simulation and thermodynamic property programs. In the above equations, xl1 and xl2 refer to the composition of the first and second liquid phases and y to that of the vapor phase. The symbol ϕ refers to the fugacity coefficient, Kvi to the K value between the vapor phase and the second liquid (L2) phase of the ith component, and Kli to that between the first (L1) and the second (L2) liquid phases. βv, βl1, and βl2 are the phase fractions of the vapor, the first liquid, and the second liquid phases, respectively. The composition of the incipient phases is determined using the following expressions: zi xil 2 = v v 1 + β (K i − 1) + β l1(K il − 1) (6)

1. INTRODUCTION Many commonly used mixtures exhibit vapor−liquid−liquid equilibria (VLLE) at low temperatures. The identification of the VLLE boundaries is of importance to those who study global phase diagrams. The identification of the VLLE region is also of importance to engineers who work with these mixtures. For example, near constant temperature refrigeration can be provided when a mixture that exhibits VLLE at low temperatures is boiled to cool a load.1 In such cases, the choice of components of the mixture and their concentration is closely linked with identifying the region where the VLLE occurs. The term three-phase dew point has been used by some authors in the literature2,3 to describe the state where two liquid and one vapor phase are in equilibrium and the number of moles (phase fraction) of one of the liquid phases in equilibrium is zero. The line connecting these states is known as the three-phase dew line, or three-phase incipient liquid line. In this article, we prefer the former nomenclature. The threephase bubble point is similar to a two-phase bubble point except that a vapor phase is in equilibrium with two liquid phases, with the phase fraction of the vapor phase being zero. A three-phase bubble point or dew point calculation involves solving the following equilibrium relations and component mass balance equations: ln K iv + ln ϕi v − ln ϕi l 2 = 0, for i = 1 to nc

(1)

ln K il + ln ϕi l1 − ln ϕi l 2 = 0, for i = 1 to nc

(2)

1

− xil 2 = 0

(3)

i=1 nc

∑ yi − xil

2

(7)

yi = xil 2K iv

(8)

The variables involved in the saturation point calculation involving two liquid phases and one vapor phase (eqs 1−3) are ln Kvi , ln Kli, ln P, ln T, βv, and βl1 (for a three-phase dew point) with β being the phase fraction. The phase fraction of the second liquid phase can be determined from that of the other two (βl2 = 1 − βv − βl1). An additional specification equation is necessary to solve eqs 1−8 for ln Kvi , ln Kli, ln P, ln T, βv, and βl1. In the conventional approach of Michelsen,4 one of ln Kvi , ln l Ki, ln P, ln T, and one of βv or βl1 is specified to solve eqs 1−5 simultaneously. Equations 1, 4, and 5 alone need to be solved in the case of a two-phase saturation point calculation. Since the equations are very similar, the first thought is to use the procedure used for two-phase envelopes for three-phase envelopes as well. In the two-phase envelope calculation of Michelsen,5 the most sensitive variable in a previous saturation point calculation is used as the specification variable in the next step. Different variables are specified at different locations of the phase-

nc

∑ xil

xil1 = xil 2K il

=0 (4)

i=1

Received: May 5, 2014 Revised: June 26, 2014 Accepted: July 9, 2014

v

β = 0 (for three‐phase bubble point) or

β l1 = 0 (for three‐phase dew point) © XXXX American Chemical Society

(5) A

dx.doi.org/10.1021/ie501838y | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

Figure 1. Phase envelope of a mixture of nitrogen/methane/ethane/propane (28.7/17.1/24.2/30.0 mol %) determined using the master phase envelope utility of Aspen Hysys v 8.4 (equation of state: Peng−Robinson).

three-phase critical point is shown in the vapor region (outside the vapor boundary or dew line). It is evident from Figure 1 and the preceding discussion that the calculation of three-phase (VLLE) boundaries is not easy even for the simple mixture of nitrogen and light hydrocarbons in Figure 1. There is a need for procedures to calculate threephase envelopes automatically. This work addresses that need. In this paper, we extend a recently developed density marching method6 for tracing two-phase envelopes to threephase envelopes. Methods have also been developed to generate the initial estimate of the composition for liquid− liquid equilibrium and an estimation of three-phase dew points, the details of which are presented in the following sections.

envelope. The step size of the specified variable is normally assumed to be the same as in the previous step. Different heuristic approaches have been developed in the literature to find the most suitable specification variable in different regions.6 Some of the above heuristic approaches can also be used for a three-phase envelope. However, the complications in the automatic calculation of three-phase envelopes according to Michelsen and Møllerup7 are (a) initial estimates for the composition of the phases at liquid−liquid equilibrium are needed and (b) calculation of liquid−liquid equilibria is less well behaved than vapor−liquid equilibria. Michelsen4 states that automatic tracing of a three-phase envelope using internally generated specifications is not easily accomplished, and divergence occurs frequently in the vicinity of critical points. Hence, Michelsen4 chooses the variable to be specified and step length set manually after each step. Cismondi and Michelsen8 presented a method for drawing the three-phase envelopes of binary mixtures. In their procedure, the calculations start from the three-phase critical point (where two of the three incipient phases have identical composition) and use very small steps to calculate the threephase bubble and dew points to construct the three-phase boundaries. The main disadvantage with this method is the need to determine the three-phase critical point to start the phase envelope construction. Since many critical points can occur in close vicinity for multicomponent mixtures that exhibit VLLE at low temperatures, the extension of the method of Cismondi and Michelsen8 to multicomponent mixtures is not easy. Similarly, the composition of two liquid phases or that of one liquid phase and one vapor phase can be the same at the three-phase critical point depending on the type of mixture according to the Van Konynenburg and Scott classification.9 The determination of the three-phase critical point a priori is therefore not straightforward. One of the few commercial programs that allows calculation of the three-phase envelopes is the Aspen Hysys process simulation program.10 Figure 1 shows the phase-envelope calculated with the latest version of Aspen Hysys (version 8.4) for a mixture of nitrogen, methane, ethane, and propane with the master envelope utility applying the Peng−Robinson equation of state. It can be seen that the three-phase dew line (incipient liquid line) has not been calculated, and the three-phase bubble line crosses the open ended dew line (wrongly marked as bubble line). It can also be seen that the

2. FUNDAMENTAL BASIS OF THE DENSITY MARCHING METHOD FOR PHASE ENVELOPES Recently, Venkatarathnam6 showed for the first time that phase envelopes should be treated as space curves. Consider the space curve shown in Figure 2. The space curve is somewhat similar to a phase envelope of a simple type I mixture on a pressure− temperature plane. The relationship between y and x variables can only be described by two equations, one for the x variable and another for the y variable in terms of an independent variable α as follows:

Figure 2. Typical space curve, similar to the p-T phase envelope of a type I mixture.9 B


Industrial & Engineering Chemistry Research x = − 12α 2 + 24α + 400

y = − α 2 − 0.2α − 20

Article

The specification equation that needs to be solved to calculate the three-phase saturation point along with eqs 1 to 5 is therefore given as follows:

(9) (10)

The relationship between α and x,y is shown graphically in Figure 3. In order to trace the x−y curve, one needs to vary α

v ρ v = ρspecified

(11)

3. DENSITY MARCHING ALGORITHM FOR THREE-PHASE ENVELOPES The density marching algorithm broadly consists of two steps: (a) determining saturation points at a specified density of an incipient phase and (b) increasing the specified density monotonically in equal steps from a low value to a high value to calculate the entire phase envelope automatically. The above procedure is identical to that adopted in ref 6 for the calculation of a two-phase envelope. However, there are many differences between calculating two and three-phase envelopes. The dew and bubble point temperatures can be calculated very easily at low pressures at a specified pressure when liquid and vapor phases are present by assuming that the vapor and liquid phases follow Dalton’s law and Raoult’s law, respectively. However, there is no easy method for determining the K values for a liquid−liquid split in the case of three-phase envelopes.7 Methods have therefore been developed to determine the three-phase bubble and dew points at low pressures reliably and are presented in the following section. One of the important steps in the calculation of the threephase envelope is to identify the phases that have zero phase fraction at the two boundaries (bubble and dew point) correctly. In this work, the two vanishing phases are designated as the V and L1 phase, respectively. The ln K values (see eqs 1 and 2) are defined with respect to the phase that does not vanish at the two boundaries (L2 phase). The identification of the L1 phase is not easy. For example, the heavier liquid phase of Mix-1 (Figure 2) vanishes at the three-phase dew point, while it is the lighter liquid phase that vanishes in the case of Mix-2 (Figure 8). In the present method, the first calculation to be performed is the calculation of the three-phase bubble point temperature at a specified low pressure. One of the liquid phases (L1 phase) vanishes on heating the saturated liquid feed to a certain vapor fraction. The heating of the saturated liquid feed is simulated by increasing the vapor fractions at a specified pressure (say 100 kPa) in steps until the calculations fail. The phase whose vapor fraction is less than zero when the calculations fail is designated the L1 phase. The last successful calculation that simulates heating of the three-phase mixture also serves as the initial value for a three-phase dew point calculation. More details of the determination of the three-phase bubble point are provided in the next section. The following procedure is recommended to calculate the entire three-phase boundary. 1. Calculate the three-phase bubble points at three specified pressures, say 80, 90, and 100 kPa, replacing eq 11 with p = pspecified 2. Identify the liquid phase that vanishes at the dew point and designate it as the L1 phase. 3. Increase the density of the y phase (or decrease that of the L1 phase) calculated in the previous iteration by a small value, say 0.1 to 1.0 mol/L, and use it in the specification equation (eq 11).

Figure 3. Relationship between the x and y variables of the space curve shown in Figure 2 and the independent variable α.

monotonically from a low value to a high value and determine the values of x and y independently using the above different equations. It is evident from the above discussion that the tracing of a phase envelope therefore requires varying an independent variable (similar to α in the above example) from a low value to a high value instead of varying one of the variables of the governing equations as in conventional methods. Many heuristics were invented to overcome the difficulties that arise at the many stationary points (∂T/∂p = 0, ∂p/∂T = 0) on the phase envelope because of the deficiency in not recognizing phase envelopes as space curves. In the density marching method of Venkatarathnam,6 the phase envelope is traced by varying the density of one of the phases, similar to varying α in Figure 3, to draw the space curve in Figure 2. The saturation point is itself calculated for the specified density of one of the incipient phases, and not at a specified pressure, temperature, or K value of one of the components as in conventional methods that follow Michelsen.5 This allows tracing the entire phase envelope by varying the density of one of the incipient phases in equal steps. This also makes it very easy to do any extrapolation to get the initial values for solving the governing equations with the Newton− Raphson method in the next step. Many examples were presented in ref 6 to demonstrate the ability of the density marching method to trace complex phase envelopes. The method has been implemented in the NIST Refprop program, version 9.1,11 released last year. A detailed description of the method is presented in ref 6 and is therefore not described again here. The density marching method can also be used for tracing three-phase envelopes. In the case of a three-phase envelope, the density of any of the incipient phases can be used as the independent variable. However, specifying the density of the vapor phase is more convenient when starting from the bubble end and proceeding toward the dew end and is therefore recommended. As in ref 6, the saturation point can is calculated at a specified density. C



Article

4. Guess the initial values of ln Kvi , ln Kli, ln P, ln T, βv, and βl1 for the next calculation by extrapolating them with respect to the density at the specified incipient phase density. 5. Solve eqs 1−7 and 11 simultaneously using a Newton− Raphson method with the above initial estimates. 6. Repeat steps 2−4 until the pressure on the dew side is lower than the minimum specified pressure (say 80 kPa). 7. The L1 liquid phase and the vapor phase get switched when the calculations cross the three-phase (critical) point. This is essentially similar to the liquid and vapor phases getting switched during a phase envelope calculation with two phases. 8. Use the phase identification method of Venkatarathnam and Oellrich12 to identify the correct phase of the different phases as in ref 6 as follows: Liquid or liquid‐like vapor: Π > 1

(12)

Vapor: Π ≤ 1

(13)

Figure 4. Difference between the three-phase bubble point and the saturation temperature of the lowest boiler for Mix-1 and Mix-2 studied in this work.

where Π is a dimensionless parameter, called the phase identification parameter, and is given by the expression: ⎡⎛ ∂ 2p ⎞ ⎛ ∂p ⎞ ⎛ ∂ 2p ⎞ ⎛ ∂p ⎞ ⎤ ⎟/⎜ ⎟ − ⎜ 2 ⎟/⎜ ⎟ ⎥ Π = v ⎢⎜ ⎢⎣⎝ ∂v∂T ⎠ ⎝ ∂T ⎠v ⎝ ∂v ⎠ ⎝ ∂v ⎠T ⎥⎦

hydrochlorfluorocarbon, hydrofluorocarbon, and hydrocarbon fluids. Multicomponent mixtures that exhibit VLLE at low temperatures are used in cryogenic systems essentially because of this feature.1 For example, the bubble point temperature at 1 bar is less than 1 K above the saturation temperature of nitrogen in the case of Mix-1 (Figure 1) even though the concentration of nitrogen is less than 30% and that of propane is more than 30%. In our procedure, the guess values for the first three calculations (three-phase bubble and dew point temperatures) are obtained as follows: 1. Assume the three-phase bubble point temperature to be 1 K above the saturation temperature of the lowest boiler at sub-atmospheric pressures. In most cases, the difference is less than 2 K even at 5 bar. 2. Perform a tangent plane stability calculation at a temperature 5 to 10 K below the saturation temperature of the lowest boiler and estimate the Kl values. The ln Kl values serve as the guess values for the full three-phase bubble point calculation. Since the fugacity coefficient of the vapor phase is one for an ideal gas, at low pressures, the ln Kv values can be l assumed to be ln ϕi2. 3. Since the phase fractions do not vary markedly with temperature in the subcooled liquid region, the phase fraction βl1 and its complementary βl2 (= 1 − βl1) obtained in the previous step in the subcooled liquid state can be used as guess values at the saturated liquid state for the full three-phase bubble point calculation. βv is zero at the three-phase bubble point. 4. When βv = 0, eq 3 reduces to the regular Rachford−Rice equation:

(14)

Or, in terms of density, the phase identification parameter (Π) is given by the following expression: ⎡⎛ ∂ 2p ⎞ ⎛ ∂p ⎞ ⎛ ∂ 2p ⎞ ⎛ ∂p ⎞ ⎤ ⎟/⎜ ⎟ − ⎜ 2 ⎟/⎜ ⎟ ⎥ Π = 2 − ρ ⎢⎜ ⎢⎣⎝ ∂ρ∂T ⎠ ⎝ ∂T ⎠v ⎝ ∂ρ ⎠ ⎝ ∂ρ ⎠T ⎥⎦

(15)

The phase identification parameter (Π) is only dependent on the partial derivatives of pressure, volume/density, and temperature and is therefore applicable at all pressures and temperatures. The calculation of the phase identification parameter (Π) at a specified density and calculated temperature is straightforward for all equations of state. With the above method, there is no need to keep track of the phase of any incipient phase during the calculation. The phase identification method is available in Aspen Hysys10 as one of the phase identification options (Venkatarathnam and Oellrich method) and can also be calculated as a property in the NIST Refprop program11 (phase identification parameter).

4. THREE-PHASE BUBBLE POINT CALCULATION Michelsen and Møllerup7 suggest solving the following equation for temperature to guess the three-phase bubble point temperature: nc

∑ pisat (T ) = p(specified) i=1

(16)

nc

∑

A simpler procedure, however, is to assume that the threephase bubble point is only 1 to 2 K above the saturation temperature of the lowest boiler at the specified temperature. Figure 4 shows the difference between the three-phase bubble points of the two mixtures (Mix-1 and Mix-2) studied in this work. It can be seen that the difference between the three-phase bubble point and the saturation temperature of nitrogen in the case of Mix-1 and that of methane in the case of Mix-2 is less than 2 K even at a pressure of about 5 bar in both cases. The difference is normally less than 1 K for most mixtures at atmospheric pressures as shown experimentally by Lavrenchenko et al.13 for different mixtures of tetrafluorocarbon,

i=1

zi(K il − 1) 1 + β l1(K il − 1)

=0 (17)

Since is known, solve the above equation for β . 5. Similarly, eq 4 reduces to the regular Rachford−Rice equation: Kli

l1

nc

f (T ):

∑ i=1

zi(K iv − 1) 1 + β l1(K il − 1)

=0 (18)

6. Assuming that β is independent of temperature but Kv and Kl are dependent on temperature (eqs 1 and 2), obtain the next temperature estimate by solving the above equation using l1

D



Article

a Newton−Raphson method once. The main assumption made here is that βl1 is independent of temperature. 7. Determine the composition of the different phases at the estimated temperature and K values using eqs 6−8. 8. Update the Kl and Kv values using the fugacity coefficients determined at the new phase compositions (eqs 1 and 2). 9. Continue the temperature iterations (steps 4 to 8) until the difference in temperature between successive iterations is less than a desired value, say 0.01 K. 10. The phase fraction βl1, converged temperature T, Kv, and l K obtained in step 9 serve as the guess values for a full saturation point calculation that solves eqs 1−5 simultaneously using a Newton−Raphson method. 11. Alternately, one can additionally solve eq 17 for the phase fraction of the first liquid phase (βl1) after solving eq 18 after step 8 using updated K values (eqs 1 and 2). The calculations are continued until eqs 1 and 2 are satisfied at the estimated temperature T and phase fraction βl1. This successive substitution procedure is slower than a full Newton−Raphson method in step 10 but is generally guaranteed to converge.

Figure 5. Dew line determined using the two-phase envelope method of ref 6 and three-phase envelope determined using the method presented in this article (equation of state: Peng−Robinson).

for the calculation. The three-phase region is very narrow, and the three-phase envelope cannot be seen well on a p−T diagram. The three-phase envelope can be better appreciated on temperature-enthalpy or temperature-entropy planes (Figures 6

5. THREE-PHASE DEW POINT The procedure described in the previous section can be modified to generate a three-phase dew point. However, since the three-phase dew point can be higher than the saturation temperature of the lowest boiler anywhere from 1 to 30 K, depending on the mixture, a good initial temperature value is not available for the three-phase dew point calculation. Convergence is therefore not guaranteed with the procedure followed for the bubble point. The following procedure is therefore used in our program. 1. Calculate the three-phase bubble point temperature using the procedure described in the previous section. 2. Solve eqs 3 and 4 simultaneously by a Newton−Raphson method for temperature T and the phase fraction βl1 for βv = 0.05. Update the compositions using eqs 6−8 and determine fugacity coefficients at the new compositions. Update the K values using eqs 1 and 2, and iterate until eqs 1 and 2 are satisfied. 3. Increase the value of βv by, say, 0.05 and repeat step 2 until no convergence is obtained. The phase fraction of one of the liquid phases is likely to be less than zero under this condition. Designate that phase as the x1 phase. Use the T, βv, Kl, and Kv of the last successful calculation and solve eqs 1−5 simultaneously for the dew point temperature (βl1 = 0). The above procedure has been found to converge even for mixtures in which the three-phase dew point temperature is more than 30 K above the three-phase bubble point temperature, for example Mix-2 in the next section. Since only three calculations are performed using the above approach, the calculation time is not large compared to the total calculation time.

Figure 6. Three-phase bubble and dew lines calculated using the method presented in this work and represented on a T−h plane.

and 7). It can be seen that the enthalpy change between the dew point (where the phase fraction of the heavy liquid phase is zero) and that at the bubble point (where the phase fraction of

6. RESULTS AND DISCUSSION Consider a mixture of nitrogen/methane/ethane/propane (28.7/17.1/24.2/30.0 mol %), Mix-1, as presented in ref 1 (Table 4.2, Case 2). This mixture exhibits VLLE at low temperatures and is therefore used in a J-T refrigerator. Figure 5 shows the dew line determined using the two-phase envelope method of ref 6 and the three-phase envelope determined using the method presented in this article. The volume translated Peng−Robinson equation of state of NIST Refprop11 was used

Figure 7. Three-phase bubble and dew lines calculated using the method presented in this work and represented on a T−s plane. E



Article

pressure. The calculations that were repeated from the bubble end showed the existence of more than one three-phase region (Figure 9) as predicted by Kohse and Heidemann. The three-

the vapor phase is zero) is over 1000 J/mol at 80 kPa. The corresponding difference in entropy is about 14 J/mol K at 80 kPa. Both the dew and bubble lines have a positive slope on the temperature−enthalpy and temperature−entropy planes. The entire phase envelope cannot be traced by varying either enthalpy or entropy monotonically. Thus, enthalpy and entropy are unsuitable as the independent specification variable (eq 11) of the problem. It is possible to not extrapolate at all and just use the values of the problem variables from the previous iteration for the next iteration but with a different specified density (eq 11). However, the calculations take a few more iterations, and convergence close to the three-phase point may be difficult in some cases. No problem was faced while crossing the threephase critical point with extrapolation. Since extrapolation is straightforward, a second order extrapolation of all the problem variables with respect to the phase whose density is specified is recommended. The density marching method presented in this article can also be used to locate the approximate three-phase critical point. A change of sign of ln Kvi or ln Klt, or the ln Kvi − ln Kli value between two successive saturation point calculations, indicates the existence of a critical point. The interpolated temperature, pressure, and density values serve as the guess values for a full three-phase critical point calculation (l1 = l2 − v or l1 − l2 = v or l2 − l1 = v). In general, the three-phase critical point calculation takes only a few iterations to converge with this approach. Calculations can begin on the dew end and proceed toward the bubble end or vice versa. Since initial three-phase dew point calculations (see section 5) require the calculation of threephase bubble points, it is advantageous to start the calculations on the bubble end and proceed toward the dew end. However, calculations do fail sometimes with this approach. Figure 8

Figure 9. Closer view of the three-phase boundaries of a mixture of CH4/C6H14/CO2/H2S (50/2.63/5.74/41.63 mol %) determined using the GERG equation of state.11 (1) Saturation point calculations started from the bubble end. (2) Saturation point calculations started on the dew end.

phase envelope procedure described in this work is not capable of tracing multiple three-phase regions since the guess values for the next saturation point calculation are generated from the previous three calculations by extrapolation. The calculations that were started from the bubble end, however, do give an indication of the complex phase behavior. The example also shows the necessity of having a method to generate the dew point calculations in mixtures that exhibit complex multiphase equilibria.

7. CONCLUSIONS The following conclusions can be drawn from this work: 1. The method presented in this paper allows automatic tracing of three-phase envelopes except in very complex mixtures exhibiting more than one three-phase region. Even in such cases, a good idea of the complex behavior can be obtained from calculations started from the bubble and the dew ends. The calculations from the opposite end can be started when the saturation pressure in the failed iteration is greater than the minimum pressure specified. 2. The density marching algorithm of ref 6 can be used to generate three-phase envelopes as well. However, an elaborate procedure needs to be followed to generate the initial values required for the Newton−Raphson method used to solve the governing equations (eqs 1−5) simultaneously. 3. It can be seen from Figure 6 that the phase envelope represented on a temperature−enthalpy diagram gives better insight into the three-phase regions compared to the traditional p−T plane used in most studies. 4. The three-phase bubble point temperature method takes advantage of the small difference between the three-phase bubble point and the saturation temperature of the low boiler in the mixture (normally less than 2 K at low pressures), and the small variation in the phase fractions of the two liquid phases with temperature, to reliably predict the three-phase bubble point. 5. The density marching method can pass through the threephase critical point easily, and only equal density steps are

Figure 8. Two- and three-phase envelopes of a mixture of CH4/ C6H14/CO2/H2S (50/2.63/5.74/41.63 mol %) determined using the GERG equation of state.11

shows one such example. The mixture of CH4/C6H14/CO2/ H2S (50/2.63/5.74/41.63 mol %), Mix-2, was studied by Kohse and Heidemann14 using the Soave−Redlich−Kwong equation of state. In this work, the GERG equation of state implemented in REFPROP11 was used to generate the three-phase envelope. The two-phase dew point calculations were limited to about 20 MPa. No effort was made to determine solid boundaries. The calculations that started on the dew end at a pressure of 5 bar stopped on the dew end before reaching the same starting F



Article

(8) Cismondi, M.; Michelsen, M. L. Global Phase Equilibrium Calculations: Critical Lines, Critical End Points and Liquid−liquid− vapour Equilibrium in Binary Mixtures. J. Supercrit. Fluids 2007, 39, 287. (9) 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, 495. (10) Aspen Hysys Process Simulation Program V 8.4; Aspen Technologies: Burlington, MA, 2013. (11) 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. (12) 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, 225− 233. (13) Lavrenchenko, G. K.; Trotsenko, A. B.; Ruvinsky, G. Y.; Tabachnik, E. U.; Balkin, B. N.; Balavachenko, B. U. Investigation of Phase Equilibria in Binary Mixtures of Refrigerants with Limits of Mixing. Kholod. Tekh. Tekhnol. 1983, 37, 38. (14) Kohse, B. F.; Heidemann, R. A. Tricritical Lines and Multiphase Equilibria in Quaternary Mixtures. Fluid Phase Equilib. 1992, 75, 11.

required. If large density steps are chosen, it may be advantageous to restrict the density variation of the vapor phase at low pressures to a smaller value, say, 0.1 or 0.2 mol/L. The method presented has direct application in process simulation and thermodynamic property programs.

■

AUTHOR INFORMATION

Corresponding Author

*Phone: +91 44 2257 4685. Fax: +91 44 2257 0509. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The author is grateful to Dr. E.W. Lemmon of the National Institute of Standards and Technology, Boulder and Prof. L.R. Oellrich, Karlsruhe Institute of Technology, Germany for reviewing this work.

■

NOMENCLATURE = K value of ith component between first liquid (L1) and the second liquid (L2) phase (see eq 2) Kvi = K value of ith component between vapor and the second liquid (L2) phase (see eq 1) nc = number of components P = pressure (MPa) T = temperature (K) V = specific volume (L/mol) l xi1 = mole fraction of ith component in the first liquid phase l xi2 = mole fraction of ith component in the second liquid phase yi = mole fraction of ith component in the vapor phase zi = mole fraction of ith component of feed Kli

Greek Letters

β = phase fraction (mol/mol) ϕ = fugacity coefficient ρ = density (mol/L) Π = phase identification parameter (eqs 12−15) Superscripts

l1 = first liquid phase l2 = second liquid phase v = vapor phase Subscript

i = ith component

■

REFERENCES

(1) Venkatarathnam, G. Cryogenic Mixed Refrigerant Processes; Springer: New York, 2008. (2) Huang, S. S.-S.; Leu, A.-D.; Ng, H.-J.; Robinson, D. B. The Phase Behavior of Two Mixtures of Methane, Carbon Dioxide, Hydrogen Sulfide, and Water. Fluid Phase Equilib. 1985, 19, 21. (3) Sampath, V. R.; Leipziger, S. Vapor-Liquid-Liquid Equilibria Computations. Ind. Eng. Chem. Process Des. Dev. 1985, 24, 652. (4) Michelsen, M. L. Some Aspects of Multiphase Calculations. Fluid Phase Equilib. 1986, 30, 15. (5) Michelsen, M. Calculation of Phase Envelopes and Critical Points for Multicomponent Mixtures. Fluid Phase Equilib. 1980, 4, 1. (6) Venkatarathnam, G. Density Marching Method for Calculating Phase Envelopes. Ind. Eng. Chem. Res. 2014, 53, 3723. (7) Michelsen, M. L.; Møllerup, J. M. Thermodynamic Models: Fundamentals & Computational Aspects; Tie-Line Publications: Denmark, 2004. G


Density Marching Method for Calculating Phase Envelopes. 2. Three

Recommend Documents