VLL Equilibria and Critical End Points Calculation of Nitrogen

Jun 19, 2012 - VLL Equilibria and Critical End Points Calculation of Nitrogen-Containing LNG Systems: Application of SRK and PC-SAFT Equations of Stat...
1 downloads 0 Views 1MB Size
Article pubs.acs.org/IECR

VLL Equilibria and Critical End Points Calculation of NitrogenContaining LNG Systems: Application of SRK and PC-SAFT Equations of State Edgar Ramírez-Jiménez,† Daimler N. Justo-García,† Fernando García-Sánchez,‡,* and Roumiana P. Stateva§ †

Departamento de Ingeniería Química Petrolera, ESIQIE, Instituto Politécnico Nacional, Unidad Profesional Adolfo López Mateos, Zacatenco, 07738 México, D.F., México ‡ Laboratorio de Termodinámica, Programa de Investigación en Ingeniería Molecular, Instituto Mexicano del Petróleo. Eje Central Lázaro Cárdenas 152, 07730 México, D.F., México § Institute of Chemical Engineering, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria ABSTRACT: Two different numerical techniques, applying either the SRK or the PC-SAFT equations of state, were used to predict the multiphase behavior exhibited by the ternary systems nitrogen + methane + ethane, nitrogen + methane + propane, and nitrogen + methane + n-butane. The predictions with both equations of state were compared against experimental data in terms of nitrogen mole fractions of the two liquid phases and it was demonstrated that the agreement is good. Furthermore, a procedure similar to that reported by Gregorowicz and de Loos, in conjunction with the SRK equation of state, was introduced. It was applied to predict the K (L1−L2=V) and LCST (L1=L2−V) critical end points at constant temperature for the systems studied, and demonstrated a reliable and robust performance with good convergence characteristics.



INTRODUCTION In processing natural gas systems, the presence of a second liquid phase in a customarily liquid−vapor mixture can often cause problems and upset the expected design performance. Though only a limited number of immiscible binary systems (methane−n-hexane, methane−n-heptane, to name the most prominent ones) are relevant to natural gas processing, vapor− liquid−liquid (VLL) behavior can and does occur under certain conditions in ternary and higher liquefied natural gas (LNG) systems even when none of the constituent binaries themselves exhibit such a behavior. It is also known that the addition of nitrogen to miscible LNG systems can induce immiscibility and this necessarily affects the process design for these systems. The interest in the use of nitrogen gas to pressurize oil reservoirs to enhance recovery has resulted in natural gas process streams, rich in nitrogen, which are likely to display complex phase behavior. Investigators have studied experimentally ternary prototype LNG systems, containing nitrogen as a second solvent, and a lot of excellent data have been published.1−4 However, to help further understanding the possible occurrence of multiphase equilibria in LNG process systems, it is necessary to extend the knowledge of their phase behavior and of the variety of critical end point boundaries through an ability to correctly predict, model, and calculate them.5,6 In a previous work,7 the phase behavior of three systems of interest in the natural gas and oil industries, namely, nitrogen + methane + (n-pentane, or n-hexane, or n-heptane) was modeled by applying two techniques. The goal was, on their examples, to examine and analyze the challenges and difficulties encountered when predicting and modeling the complex phase behavior of LNG systems. Furthermore, two representatives of the © 2012 American Chemical Society

equations of state (EoSs)-type thermodynamic models, namely, SRK EoS8 and PC-SAFT EoS,9 were used, and their capabilities to represent the equilibrium fluid phases were compared. The aim of this study is 2-fold: The first aim is to exhaust the list of model nitrogen-containing ternary systems, where methane is the solvent, nitrogen is a second solvent, and any of the n-alkanes (from ethane to n-heptane) is the solute that display VLLE behavior and to predict and calculate their phase behavior applying two different computational procedures. A concise analysis of the numerical performance of the techniques is offered as well. The second aim is to extend the modeling objectives by including the calculation of K (upper critical end point) and/or LCST (lower critical solution temperature) critical end points of the systems studied. To realize the objectives, three systems, not studied previously, namely, nitrogen + methane + (ethane or propane or n-butane) were examined and their phase behavior predicted, and an algorithm, similar to that proposed by Gregorowicz and de Loos,10 was advocated and implemented to calculate the systems’ critical end points.



NUMERICAL PROCEDURES AND ALGORITHMS Phase Equilibrium Calculations. The numerical techniques used to predict and model the phase behavior and calculate the components’ distribution among the equilibrium phases for the three model systems were described in detail in

Received: Revised: Accepted: Published: 9409

February 10, 2012 June 14, 2012 June 19, 2012 June 19, 2012 dx.doi.org/10.1021/ie300372a | Ind. Eng. Chem. Res. 2012, 51, 9409−9418

Industrial & Engineering Chemistry Research

Article

our previous work.7 In view of this, a brief description only will be given here: Technique 1 employs the PC-SAFT EoS9 as the thermodynamic model, and uses an efficient computational procedure for solving the isothermal multiphase problem by assuming that the system is initially monophasic. A stability test allows verification of whether the system is stable or not. In the latter case, it provides an estimation of the composition of an additional phase; the number of phases is then increased by one, and equilibrium is achieved by minimizing the Gibbs energy. This approach, advocated as a stagewise procedure,11,12 is continued until a stable solution is found. In this technique, the stability analysis of a homogeneous system of composition z, based on the minimization of the distance separating the Gibbs energy surface from the tangent plane at z, is considered.13,14 In terms of fugacity coefficients, φi, this criterion for stability can be written as14

the quasi-Newton BFGS15 method which ensures the property of strict descent of the Gibbs energy surface. A detailed description of this approach for solving the isothermal multiphase problem can be found elsewhere.16 Technique 2 uses the SRK cubic EoS8 as the thermodynamic model and applies a rigorous thermodynamic stability analysis which is followed by a simple and effective procedure for identifying the phase configuration at equilibrium corresponding with the minimum Gibbs energy. The stability analysis is exercised once and on the initial system only; it is based on the well-known tangent plane criterion13,14 but uses a different objective function.17−19 The key point is to locate all zeros (y*) of a function Φ(y) given as N

Φ(y) =

where

∑ ξi[ln ξi + ln φi(ξ) − hi − 1] ≥ 0

ki(y) = ln φi(y) + ln yi − hi

i=1

∀ξ>0

(1)

i = 1, ..., N

(2)

Equation 1 requires that the tangent plane, at no point, lies above the Gibbs energy surface, and this is achieved when F(ξ) is positive in all its minima. The quasi-Newton BFGS minimization method15 is applied to eq 1 for determining the stability of a given system of composition z at specified temperature and pressure. Once instability is detected with the solution at π − 1 phases, the equilibrium calculation is solved by minimization of the following function π

Δg = min (ϕ) ni

N

⎛ x (ϕ)ϕ(ϕ)p ⎞ i i ⎟ ⎟ ° p ⎝ ⎠

∑ ∑ ni(ϕ) ln⎜⎜

ϕ=1 i=1

ki* = ln yi* + ln φ(y*) − hi

(3)

π−1

i = 1, ..., N (4)

ϕ=1

and ni(ϕ) ≥ 0

i = 1, ..., N ;

ϕ = 1, ..., π − 1

(7)

i = 1, ..., N

(8)

Furthermore, the number k*, which geometrically is the distance between two such hyperplanes, can be either positive or negative. When all calculated k* are positive, the initial system is stable; if at least one zero, to which a negative k* corresponds is located, it is unstable. The specific form of the modified Φ(y) (its zeros are its minima) and the fact that it is easily differentiated analytically allows the application of a local solver for determining its stationary points. Hence, the quasi-Newton BFGS method with a line search and superlinear order of convergence is used with a robust initialization strategy.17−19 The phase stability analysis performed provides also good approximations to the compositions to which the subsequent flash calculations of the “phase identification procedure” converge. Consequently, although local in nature, the flash routines neither diverge nor exhibit a tendency toward a strong attraction to the trivial solution. The numerical implementation, practical application, versatility, and robustness of Technique 1 and Technique 2 have been demonstrated over the years on a very large number of case studies, quite different in nature−from highly difficult theoretical problems to large practical problems in petroleum and chemical engineering. Some of the systems examined were chosen primarily for the high degree of complexity, which provided a rigorous test for the techniques’ capabilities. Further details on their performance can be found elsewhere.16,20−24 Calculation of K- and LCST-Points. Ternary systems which exhibit LLV behavior but do not exhibit such behavior in their constituent binaries have the immiscibility region bounded by a K (L1−L2=V)-point locus, a LCST (L1−L2=V) locus, and a Q (S−L1−L2−V)-point locus. Ternary systems which have immiscibility in a constituent binary can have boundaries similar to those mentioned above, besides the intrusion of the

subject to the inequality constraints given by

∑ ni(ϕ) ≤ zi

i = 1, ..., N

with hi given by eq 2 and assuming kN+1(y) = k1(y). From eqs 6 and 7, it follows that min Φ(y) = 0 when k1(y*) = k2(y*) = ... = kN(y*). The zeros of Φ(y) conform to points on the Gibbs energy hypersurface, where the local tangent hyperplane is parallel to that at z. To each zero y*, a number k* (equal for each yi* of a zero of the functional) corresponds, such that

where ξi are mole numbers with corresponding mole fractions as yi = ξi/∑j =N 1ξj, and hi = ln zi + ln φi(z)

(6)

i=1

N

F (ξ ) = 1 +

∑ [ki+ 1(y) − ki(y)]2

(5)

where Δg is the dimensionless Gibbs energy of the system, zi is the mole fraction of the component i in the system, n(ϕ) (i = 1, i ..., N ; ϕ = 1, ..., π − 1) is the mole number of component i in phase ϕ per mole of feed, x(ϕ) is the mole fraction of i component i in phase ϕ, T is the temperature, p is the pressure, and p° is the pressure at the standard state of 1 atm (101.325 (π) (π) kPa). In eq 3 the variables n(π) are considered i , xi , and φi functions of n(ϕ) . i Equation 3 is solved using an unconstrained minimization algorithm by keeping the variables n(ϕ) inside the convex i constraint domain given by eqs 4 and 5 during the search for the solution. In this case, a hybrid approach is used to minimize eq 3 starting with the steepest-descent method in conjunction with a robust initialization supplied from the stability test to ensure a certain progress from initializations, and ending with 9410

dx.doi.org/10.1021/ie300372a | Ind. Eng. Chem. Res. 2012, 51, 9409−9418

Industrial & Engineering Chemistry Research

Article

Another approach to calculate the K- and LCST-points was proposed by Mushrif and Phoenix.26 This approach utilizes an efficient critical point solver27 and a standard phase stability test14 within a nested-loop structure to directly locate K- and LCST-points. The algorithm consists of two nested inner loops to calculate a critical-point temperature and volume at fixed composition z. An outer loop uses the critical point as a test phase, searches for an incipient phase at a trial composition n̂, and updates the critical composition to iteratively decrease the tangent plane distance of the incipient phase to zero. The Newton−Raphson method with numerical derivatives was used in both the inner and outer loops. In this work, we have adopted the approach of Gregorowicz and de Loos10 to calculate the ternary K- and LCST-points rather than that of Mushrif and Phoenix26 because in the former the calculations are performed at chosen values of the temperature. This is the usual situation in practice where the experimental three-phase VLLE data are measured isothermally. We apply the Newton iteration technique to solve a set of only five nonlinear equations, namely, the conditions of criticality given by Baker and Luks,28

binary LLV locus on the ternary LLV region. The K-point and LCST loci can intersect at a tricritical point where the three phases become critical; i.e., L1=L2=V. In a study on the modeling of the three-phase VLL region for ternary hydrocarbon mixtures with the SRK EoS, Gregorowicz and de Loos10 proposed a procedure for finding K- and LCSTpoints of ternary systems, based on the solution of thermodynamics conditions for the K- and LCST-points by using the Newton iteration technique and carefully chosen starting points. The procedure was applied to calculate the Kand LCST-point loci for two ternary systems, namely, C2 + C3 + C20 and C1 + C2 + C20, in which the constituent binary C2 + C20 exhibit immiscibility. Consequently, the extension of the three-phase LLV region of these systems is bounded by the binary L1−L2−V locus of the system C2 + C20 and the ternary K-point and L-point loci. Briefly, the strategy followed to find the critical end points was the following: (1) calculation of the critical line, the Kpoint, the three phase line, and the LCST-point for the system C2 + C20 using thermodynamic conditions, and (2) calculation of the K- and LCST-point loci for the ternary systems by using as starting points the coordinates of the K- and LCST-points found for the binary system C2 + C20. In this case, to obtain a K- and LCST-point for a ternary system, the following set of six nonlinear equations, AVV

AVx1

Dx1

D2 =

(9)

xα1 ,

xα2 ,

in seven variables, T, V , V , and have to be solved, where α designates either V for the L-point or L2 for the Kpoint, D and D* are the two determinants that must be satisfied at a critical point c, and μi is the chemical potential of component i. The ternary K- and LCST-points were calculated at chosen values of the temperature. The critical criteria given by eqs 1 and 2 are based on the Helmholtz energy, which can be expressed as,25 A=

∫V

∞⎛ ⎜



p−

N ⎛ V ⎞ nRT ⎞⎟ dV − RT ∑ ni ln⎜ ⎟ ⎠ V ⎝ niRT ⎠ i=1

N

+

∑ ni(ui0 − Tsi0) i=1

(14)

=0

Dn2

(15)

i = 1, 2, 3

xα1 ,

(16)

xα2 .

in five variables, p, and In eq 15, DV and Dni denote the derivatives of determinant D1 with respect to volume and number of moles of component i (i = 1, 2), respectively. The difference between the procedure of Gregorowicz and de Loos10 and the one reported here for finding K- and LCSTpoints of ternary systems at constant temperature, is in the manner of calculating the values of the two determinants that must be satisfied at the critical point. The method used by us to calculate the determinants D1 and D2 (eqs 14 and 15) was that proposed by Baker and Luks.28 In this method, the row ordering of the array of determinants D1 and D2 were chosen in such a way that the arrays are identical, except for the last row. The determinant D1 is then evaluated by computing the cofactors of the last row, which are, in turn, used to evaluate the determinant D2. In addition, Baker and Luks showed that the derivative of determinant D1 with respect to a given variable (e.g., V) is the sum of the element-by-element products of two matrices: one matrix contains the cofactors of the original matrix, whereas the other contains the derivative of the original matrix elements with respect to the variable. This method allows saving computer time, in particular for multicomponent systems. An additional important difference between the two procedures is that in the Gregorowicz and de Loos’ procedure a set of six nonlinear equations are solved in six variables (Vc, Vα, xc1,xc2, xα1 , and xα2 ), whereas in our procedure a set of five

(11)

xc2,

Dn1

xc1,

(12)

xc1,

AVn1 A n1n1 A n1n2

f ic − f iα = 0

(10)

pc − pα = 0 α

=0

and the equilibrium conditions in terms of fugacities, f i,

=0

i = 1, 2, 3

c

AVn1 A n1n1 A n1n2

DV

AVx2 A x1x2 A x2x2

μic − μiα = 0

AVn2

AVV AVn1 AVn2

=0

Dx2

AVx1 A x1x1 A x1x2

D* =

AVn1

AVn2 A n1n2 A n2n2

AVx2 A x1x2 A x2x2 DV

D1 =

AVx2

AVx1 A x1x1 A x1x2

D=

AVV

(13)

where p is the pressure, T is the temperature, V is the system volume, ni is the mole number of component i, u0i is the standard state molar internal energy, s0i is the standard state molar entropy, and R is the gas constant. Derivatives of the Helmholtz energy are denoted by a subscript in eqs 9 and 10 (e.g., AV, Axi) indicating the differentiation variable (volume V or mole fraction xi of component i). 9411

xc2,

dx.doi.org/10.1021/ie300372a | Ind. Eng. Chem. Res. 2012, 51, 9409−9418

Industrial & Engineering Chemistry Research

Article

nonlinear equations are solved in five variables (p, xc1,xc2, xα1 , and xα2 ), which makes the iterative process easier to perform and converge. In what follows a brief description of the algorithm realizing the Newton iteration technique employed is given: Once the derivatives of the Helmholtz energy with respect to volume and mole number are known, the two determinants D1 and D2 that must be satisfied at a critical point are evaluated at the critical phase composition for a given pressure and constant temperature, according to the procedure suggested by Baker and Luks.28 In this process, the volume of the desired phase (liquid or vapor, depending on the type of critical end point calculation) is calculated using the procedure of Topliss et al.29,30 at the same conditions of temperature and pressure. The pressure is then perturbed (with a value of 10−10, for example) and the determinants D1 and D2 are evaluated again at the same critical composition and temperature to determine numerically the derivative of D1 and D2 with respect to pressure. Similarly, the critical phase compositions xc1 and xc2 are perturbed and determinants D1 and D2 are evaluated at the same initial pressure and constant temperature to determine numerically the derivative of D1 and D2 with respect to xc1 and xc2. The derivatives of the fugacity with respect to pressure and the critical, x1c and x2c , and noncritical, x1α and x2α, phase compositions, are evaluated analytically using the volumes calculated with the procedure of Topliss et al. to evaluate the fugacities. The set of linear equations is then solved using LU decomposition.31 The corrections are added to the variables p, xc1,xc2, xα1 , and xα2 according to the Newton method. The process is iterative until convergence within the desired tolerance is achieved. Though, in principle, Newton methods are strongly initialization dependent and can either oscillate or diverge from the fixed point, our algorithm demonstrates good convergence characteristics owing to the fact that values of the variables p, xc1,xc2, xα1 , and xα2 obtained from the VLLE calculations, are used as initial estimates.

The Nitrogen + Methane + Ethane System. The threephase VLL region displayed by this ternary system is bounded from above by a K-point locus and from below by a lower critical solution temperature (LCST) locus, while at low temperatures by a Q-point locus. Owing to the fact that the system contains a binary pair (nitrogen + ethane) which exhibits VLL behavior,37 its VLL space is truncated. In this case, the partially miscible pair nitrogen + ethane spans the VLL space from a position of the LCST locus to a position on the Q-point space. Because methane is of intermediate volatility compared with nitrogen and ethane, it creates a three-phase VLL space which extends from the binary VLL locus upward in temperature. The topographical nature of the regions of immiscibility for the system nitrogen + methane + ethane is shown in Figure 1. In this figure, it can be seen that the L−L=V and L=L−V critical end-point loci intersect at a tricritical point (L=L=V).

Figure 1. Experimental boundaries of the L1−L2−V immiscibility region for the system nitrogen + methane + ethane.1,37



RESULTS AND DISCUSSION The phase behavior of the three systems of interest was modeled applying computational techniques 1 and 2, and the results were compared with the experimental data reported by Llave et al.1 for the systems nitrogen + methane + ethane and nitrogen + methane + propane, and by Merrill et al.2 for the system nitrogen + methane + n-butane. The binary interaction parameters used in this work for the SRK EoS were taken from Knapp et al.32 and from Nagy and Shirkovskiy,33 and their values are: kN2−C1 = 0.0278, kN2−C2 = 0.0407, kN2−C3 = 0.0763, kN2−C4 = 0.07, kC1−C2 = −0.0078, kC1−C3 = 0.009, and kC1−C4 = 0.0056, while the corresponding binary interaction parameters for the PC-SAFT equation are: kN2−C1 = 0.0307, kN2−C2 = 0.0458, kN2−C3 = 0.0759, kN2−C4 = 0.057, kC1−C2 = 0.0039, kC1−C3 = 0.0019, and kC1−C4 = 0.0192. The latter values are those of ́ Garci a-Sá n chez et al.34 and Justo-Garciá et al.35 The components’ physical properties required for the calculations performed with the SRK EoS were taken from DIPPR36 while the three pure-component parameters (i.e., temperature independent segment diameter σ, depth of the potential ε, and number of segments per chain m) of these compounds for the PC-SAFT equation of state were taken from Gross and Sadowski.9

Figures 2 and 3 present the experimental and calculated L1− L2−V phase behavior (in terms of L1−L2 nitrogen mole fraction data) for the nitrogen + methane + ethane system at, respectively, 125 and 129 K and different pressures. These figures show a reasonable agreement between the experimental values of liquid phases L1 and L2 and those calculated by both models. Figure 1 shows the locus of K-points for this ternary system at higher temperatures; hence, it would be interesting to predict the K-point and the LCST loci with the SRK and PC-SAFT equations by using the same interaction parameters as those used in the VLLE calculations. An attempt to directly calculate the LCST points at T = (125 and 129) K for this ternary system was carried out by using the procedure described above. As pointed out before, the use of good initial guesses, obtained from the VLL calculations close to the LCST, promotes the convergence of the Newton procedure and correct values of the critical end points for the SRK equation were obtained. The predicted LCST points at 125 K (2.402 MPa and 0.5694 N2 mole fraction) and 129 K (2.837 MPa and 0.5760 N2 mole fraction), are shown in Figures 2 and 3, respectively. Figures 2 and 3 also show that at pressures away from the LCST point, the PC-SAFT model gives a better representation of the experimental compositions for liquid phase L1 while both 9412

dx.doi.org/10.1021/ie300372a | Ind. Eng. Chem. Res. 2012, 51, 9409−9418

Industrial & Engineering Chemistry Research

Article

temperatures are presented in Table 1. The table includes the pressures estimated and the critical and noncritical mole fractions of the three components. Figure 4 shows the experimental and calculated K-point and LCST loci for this system on a pressure−temperature phase diagram. As can be seen, the three-phase VLL region bounded by these two loci is predicted at lower pressures than the experimental ones. This figure also shows that the position of the tricritical point is found at a temperature and pressure higher than the “hypothetical” experimental one. Figure 5 shows on a temperature−composition phase diagram, the variation of the mole fractions of nitrogen, methane, and ethane as temperature increases: the mole fraction of nitrogen decreases as temperature increases for the K-point locus, whereas it increases as temperature increases for the LCST locus; for methane and ethane the respective mole fractions increase as temperature increases for the K-point locus and decrease as temperature increases for the LCST locus. It is also evident from Figure 5 and Table 1 that the mole fraction of each component at the K-point and LCST loci approaches the tricritical point as temperature increases: for example at T = 149.5 K the difference in pressure between the calculated Kpoint and LCST is 0.003 MPa. Although the deviations between the experimentally measured and the calculated K- and LCST-points with the SRK EoS can be attributed to the fact that binary interaction parameters used were determined from binary vapor−liquid equilibrium data, it can be said that, on the whole, this equation gives a satisfactory representation of the K-point and LCST loci for this ternary system. The Nitrogen + Methane + Propane System. This system is topographically similar to the nitrogen + methane + ethane system, sharing the same sort of boundaries; that is, Kpoint, Q-point, LCST, and binary (nitrogen + propane) VLL loci.37 However, its three-phase VLL region extends over a much larger area and upward in temperature and pressure from the binary VLL locus in comparison to the nitrogen + methane + ethane system, as shown in Figure 6. This figure shows that the LCST and K-point loci intersect at a tricritical point. Notice that the experimental Q-point locus for this system is not shown due to the limitations of the experimental apparatus used by Llave et al.1 for measuring these Q points. Figures 7−9 present, respectively, the experimental and calculated V−L1−L2 phase behavior (in terms of L1−L2 nitrogen mole fraction data) for the nitrogen + methane + propane system at the temperatures of 125, 127, and 134 K, and different pressures. Overall, there is a satisfactory agreement between the experimental values of liquid phases L1 and L2 and those predicted with both models. It should be noted that the VLL calculations were performed up to almost the position of the binary VLL boundary at T = (125 and 127) K, where methane mole fractions are practically zero in all three phases; that is, this ternary system exhibits LCST but not Kpoints at these temperatures. At pressures away from the LCST, the SRK EoS gives a better representation of the experimental compositions for liquid phase L1 while both models agree closely with each other for the liquid L2 phase and give an adequate representation of the experimental data (Figures 7 and 8). Conversely, Figure 9 shows that when pressure approaches the LCST point, the PCSAFT EoS gives a better representation of the experimental compositions for liquids L1 and L2 at 134 K. Since there are no experimental data below 2.216 MPa at 125 K and below 2.335

Figure 2. Experimental1 and calculated N2 mole fractions for the coexisting L1 and L2 phases of the system nitrogen (1) + methane (2) + ethane (3) at 125 K.

Figure 3. Experimental1 and calculated N2 mole fractions for the coexisting L1 and L2 phases of the system nitrogen (1) + methane (2) + ethane (3) at 129 K.

models agree with each other for the liquid L2 phase and represent the experimental data closely. Since there are not experimental data below p = (2.615 and 3.167) MPa at T = (125 and 129) K, respectively, comparisons of the models with experiment in the region of the LCST are not possible. However, according to the predictions with both models at these temperatures, the “estimated” LCST points with the PCSAFT model (2.46 and 2.94 MPa) seems to be in a better agreement with the “hypothetical” experimental LCST points (2.52 and 2.92 MPa) in comparison with the LCST points calculated from the SRK model (2.402 and 2.837 MPa). Of course, these discrepancies could be due to the fact that binary interaction parameters determined from VL equilibrium data were used, which, apparently, led to less accurate results. Further, the capability of the procedure to predict the Kpoint and LCST loci of ternary systems was tested using the SRK equation along with the interaction parameters used in the VLLE calculations. Results of the K- and LCST-points calculations for the nitrogen + methane + ethane at different 9413

dx.doi.org/10.1021/ie300372a | Ind. Eng. Chem. Res. 2012, 51, 9409−9418

Industrial & Engineering Chemistry Research

Article

Table 1. Predicted (SRK EoS) K-Point and LCST Loci for the N2 (1) + CH4 (2) + C2H6 (3) System critical phase mole fractions T (K)

p (MPa)

x1

x2

133.0 135.0 135.7 140.0 142.5 145.0 146.0 146.5 147.0 147.3 147.6 148.0 148.5 149.5

3.982 4.138 4.343 4.558 4.786 5.027 5.128 5.179 5.232 5.264 5.296 5.339 5.395 5.505

0.9462 0.9223 0.8920 0.8611 0.8289 0.7943 0.7794 0.7716 0.7634 0.7583 0.7530 0.7456 0.7356 0.6819

0.0431 0.0653 0.0928 0.1203 0.1481 0.1766 0.1884 0.1944 0.2006 0.2044 0.2083 0.2136 0.2206 0.2514

125.0 129.0 133.0 135.0 137.5 140.0 142.5 145.0 146.0 146.5 147.0 147.3 147.6 148.0 148.5 149.5

2.402 2.837 3.309 3.551 3.879 4.211 4.550 4.893 5.031 5.100 5.168 5.209 5.250 5.305 5.372 5.502

0.5684 0.5759 0.5838 0.5880 0.5936 0.6001 0.6079 0.6182 0.6234 0.6264 0.6298 0.6320 0.6345 0.6380 0.6433 0.6683

0.2523 0.2526 0.2545 0.2562 0.2585 0.2612 0.2640 0.2663 0.2669 0.2671 0.2671 0.2670 0.2668 0.2665 0.2657 0.2578

noncritical phase mole fractions x3

x1

x2

x3

0.0107 0.0124 0.0152 0.0186 0.0231 0.0291 0.0322 0.0340 0.0360 0.0373 0.0387 0.0408 0.0438 0.0667

0.3372 0.3507 0.3821 0.4165 0.4549 0.4997 0.5204 0.5316 0.5436 0.5513 0.5593 0.5708 0.5868 0.6373

0.1065 0.1508 0.1972 0.2342 0.2621 0.2805 0.2847 0.2861 0.2868 0.2668 0.2866 0.2859 0.2838 0.2707

0.5663 0.4985 0.4207 0.3493 0.2830 0.2198 0.1949 0.1823 0.1696 0.1619 0.1541 0.1433 0.1294 0.0920

0.1793 0.1715 0.1617 0.1558 0.1479 0.1387 0.1281 0.1155 0.1097 0.1065 0.1031 0.1010 0.0987 0.0955 0.0910 0.0739

0.9441 0.9328 0.9182 0.9093 0.8959 0.8796 0.8590 0.8318 0.8181 0.8103 0.8018 0.7961 0.7903 0.7813 0.7689 0.7425

0.0552 0.0660 0.0797 0.0879 0.1000 0.1145 0.1322 0.1546 0.1654 0.1714 0.1778 0.1821 0.1863 0.1928 0.2014 0.2192

0.0007 0.0012 0.0021 0.0028 0.0041 0.0059 0.0088 0.0136 0.0165 0.0183 0.0204 0.0218 0.0234 0.0259 0.0297 0.0383

K-Point

LCST

Figure 4. Experimental1 and calculated boundaries of the L1−L2−V immiscibility region for the system nitrogen + methane + ethane.

Figure 5. Mole fractions of nitrogen, methane, and ethane at the predicted (SRK EoS) K-point and LCST loci of the system nitrogen (1) + methane (2) + ethane (3).

MPa at 127 K, comparisons of the models with experiment in the region of the LCST are not possible. However, according to the predictions with both models, the estimated (by interpolation) LCST with the PC-SAFT, namely (1.877, 2.055, and 2.757) MPa for the temperatures of (125, 127, and 134) K, respectively, seems to be closer to the

“hypothetical” experimental LCST point in comparison with the LCST points calculated from the SRK (1.806, 1.964, and 2.572 MPa) for the temperatures of (125, 127, and 134) K. Figure 9 also shows the K-point calculated with the SRK; its position is consistent with the predicted liquid L2 and vapor V loci obtained from the VLLE calculations. Hence, the 9414

dx.doi.org/10.1021/ie300372a | Ind. Eng. Chem. Res. 2012, 51, 9409−9418

Industrial & Engineering Chemistry Research

Article

Figure 6. Experimental boundaries of the L1−L2−V immiscibility region for the system nitrogen + methane + propane.1,37.

Figure 8. Experimental1 and calculated N2 mole fractions for the coexisting L1 and L2 phases of the system nitrogen (1) + methane (2) + propane (3) at 127 K.

Figure 7. Experimental1 and calculated N2 mole fractions for the coexisting L1 and L2 phases of the system nitrogen (1) + methane (2) + propane (3) at 125 K.

Figure 9. Experimental1 and calculated N2 mole fractions for the coexisting L1 and L2 phases of the system nitrogen (1) + methane (2) + propane (3) at 134 K.

procedure we have introduced can be used with confidence to predict critical end points at constant temperature. The Nitrogen + Methane + n-Butane System. As mentioned above, the type of the VLL region displayed by the ternary systems depends on whether they contain an immiscible pair or not. In this context, the system nitrogen + methane + n-butane does not exhibit immiscibility in any of their binary pairs. Consequently, the three-phase region displayed by this system is “triangular”, which is bounded from above by a K-point locus (L−L=V), from below by a LCST (L=L−V) locus, and by a Q-point locus (S−L−L−V) at low temperatures. These three loci intersect at invariant points for this system: the K-point and LCST loci, at a tricritical point, whereas the Q-point locus terminates at a point of the type S− L−L=V from above and at a point of the type S−L=L−V from below, respectively. Figure 10 presents the experimental pressure−temperature diagram of the three-phase VLL space developed by this ternary system.2 This is a very interesting system because it is in contrast to the earlier systems, which did not exhibit L1−L2−V behavior outside the boundary loci when it is projected on pressure−temperature space. That is, above

150 K the critical solution temperature (CST) locus is composed of LCST points, similar to those observed in other nitrogen-containing ternary systems.1,3,4 Below 150 K the CST boundary changes to an UCST (upper critical solution temperature) locus in a manner similar to that exhibited by the systems methane +2-methylpentane +2-ethyl-1-butene and methane +3,3-dimethylpentane +2-methylhexane reported by van Konynenburg.38 Figure 10 presents the topographical nature of the regions of immiscibility for the system nitrogen + methane + n-butane. As can be seen, the L−L=V and L=L−V critical end-point loci intersect at a tricritical point (L=L=V). For this system, VLL calculations were performed at three temperatures T = (128, 139, and 149) K and the results obtained are presented in Figures 11−13. Figure 11 shows that both thermodynamic models predict well the experimental nitrogen compositions of liquid L2 phase down to the lowest experimentally measured pressure of 1.658 MPa; however, the SRK model is superior to the PC-SAFT model. Below this pressure, the nitrogen compositions predicted in both the L1 and L2 phases with PC-SAFT differ 9415

dx.doi.org/10.1021/ie300372a | Ind. Eng. Chem. Res. 2012, 51, 9409−9418

Industrial & Engineering Chemistry Research

Article

Figure 10. Experimental boundaries of the L1−L2−V immiscibility region for the system nitrogen + methane + n-butane.2

Figure 12. Experimental2 and calculated N2 mole fractions for the coexisting L1 and L2 phases of the system nitrogen (1) + methane (2) + n-butane (3) at 139 K.

Figure 11. Experimental2 and calculated N2 mole fractions for the coexisting L1 and L2 phases of the system nitrogen (1) + methane (2) + n-butane (3) at 128 K.

Figure 13. Experimental2 and calculated N2 mole fractions for the coexisting L1 and L2 phases of the system nitrogen (1) + methane (2) + n-butane (3) at 149 K.

considerably from the experimental data, while the SRK predicts the experimental values very closely. These unexpected results may be attributed again to the interaction parameter used for the nitrogen + n-butane binary by the PC −SAFT. Because at 128 K there is not any K-point, only the UCST was calculated with the SRK model, which gave a value of 1.712 MPa. No attempt was made to estimate (by interpolation, for example) this critical end point with the PC-SAFT EoS because this point is far away from the “hypothetical” experimental UCST. Figures 12 and 13 present the VLL nitrogen compositions predicted at 139 and 149 K, respectively. Both thermodynamic models give adequate predictions with respect to the experimentally measured nitrogen compositions of the liquid L2 phase above to the lowest pressure measured (p = 2.464 MPa for the isotherm of 139 K). Both models predict similar values for nitrogen compositions of the liquid L1 phase, but these compositions differ from the experimental data down to the lowest measured pressure of 2.512 MPa. This figure also shows that the liquid L1 and L2 nitrogen compositions

predicted with the PC-SAFT follow a trend similar to that shown in Figure 11, though closer to the experimental data. Since there are not any experimental UCST data for this isotherm, it was possible to estimate it by interpolation. The “hypothetical” experimental UCST value found is 2.44 MPa while the calculated UCST with the SRK model is 2.499 MPa, which is only 0.06 MPa above the “experimental” value. The estimated (by interpolation) UCST with the PC-SAFT is 2.21 MPa, still away from the estimated “experimental” value. The K-point calculated with the SRK model is 4.159 MPa, but because there are not any experimental data for the nitrogen compositions of the vapor phase, a comparison with the experiment in the region of the K-point at this temperature cannot be realized. A similar value for this critical end point seems to be obtained with the PC-SAFT. At T = 149 K the liquid L1 and L2 nitrogen mole fractions predicted with the SRK are in a better agreement with the experimental liquid L1 data than those predicted with the PCSAFT (Figure 13). However, the liquid L2 nitrogen mole 9416

dx.doi.org/10.1021/ie300372a | Ind. Eng. Chem. Res. 2012, 51, 9409−9418

Industrial & Engineering Chemistry Research

Article

Nitrogen + Methane + Propane Systems. J. Chem. Eng. Data 1987, 32, 14−17. (2) Merrill, R. C.; Luks, K. D.; Kohn, J. P. Three-Phase Liquid− Liquid−Vapor Equilibria in the Methane + n-Butane + Nitrogen System. Adv. Cryog. Eng. 1984, 29, 949−955. (3) Merrill, R. C.; Luks, K. D.; Kohn, J. P. Three-Phase Liquid− Liquid−Vapor Equilibria in the Methane + n-Hexane + Nitrogen + and Methane + n-Pentane + Nitrogen Systems. J. Chem. Eng. Data 1984, 29, 272−276. (4) Chen, W.-L.; Luks, K. D.; Kohn, J. P. Three-Phase Liquid− Liquid−Vapor Equilibria in the Nitrogen + Methane + n-Heptane System. J. Chem. Eng. Data 1989, 34, 312−314. (5) Luks, K. D.; Merrill, R. C.; Kohn, J. P. Partial Miscibility Behavior in Cryogenic Natural Gas Systems. Fluid Phase Equilib. 1983, 14, 193− 201. (6) Luks, K. D. The Occurrence and Measurement of Multiphase Equilibria Behavior. Fluid Phase Equilib. 1986, 29, 209−224. (7) Justo-García, D. N.; García-Sánchez, F.; Stateva, R. P.; GarcíaFlores, B. E. Modeling of the Multiphase Behavior of NitrogenContaining Systems at Low Temperatures with Equations of State. J. Chem. Eng. Data 2009, 54, 2689−2695. (8) Soave, G. Equilibrium Constants from a Modified Redlich− Kwong Equation of State. Chem. Eng. Sci. 1972, 27, 1107−1203. (9) Gross, J.; Sadowski, G. Perturbed-Chain SAFT: An Equation of State Based on Perturbation Theory for Chain Molecules. Ind. Eng. Chem. Res. 2001, 40, 1244−1260. (10) Gregorowicz, J.; de Loos, Th. W. Modelling of the Three Phase LLV Region for Ternary Hydrocarbon Mixtures with the Soave− Redlich−Kwong. Fluid Phase Equilib. 1996, 118, 121−132. (11) Michelsen, M. L. The Isothermal Flash Problem. Part II. PhaseSplit Calculation. Fluid Phase Equilib. 1982, 9, 21−40. (12) Nghiem, L. X.; Li, Y.-K. Computation of Multiphase Equilibrium Phenomena with an Equation of State. Fluid Phase Equilib. 1984, 17, 77−95. (13) Baker, L. E.; Pierce, A. C.; Luks, K. D. Gibbs Energy Analysis of Phase Equilibria. Soc. Pet. Eng. J. 1982, 22, 731−742. (14) Michelsen, M. L. The Isothermal Flash Problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, 1−19. (15) Fletcher, R. Practical Methods for Optimization, 2nd ed.; Wiley Interscience: New York, 1987. (16) Justo-García, D. N.; García-Sánchez, F.; Romero-Martínez, A. Isothermal Multiphase Flash Calculations with the PC-SAFT Equation of State. Am. Inst. Phys. Conf. Proc. 2008, 979, 195−214. (17) Stateva, R. P.; Tsvetkov, G. A Rigorous Approach to Stability Analysis as a First Step When Solving the Isothermal Multiphase Flash Problem. Technol. Today 1991, 4, 233−240. (18) Stateva, R. P.; Tsvetkov, S. G. A Diverse Approach for the Solution of the Isothermal Multiphase Flash Problem. Application to Vapor−Liquid−Liquid Systems. Can. J. Chem. Eng. 1994, 72, 722− 734. (19) Wakeham, W. A.; Stateva, R. P. Numerical Solution of the Isothermal Multiphase Flash Problem. Rev. Chem. Eng. J. 2004, 20, 1− 56. (20) García-Sánchez, F.; Schwartzentruber, J.; Ammar, N. M.; Renon, H. Modeling of Multiphase Liquid Equilibria for Multicomponent Mixtures. Fluid Phase Equilib. 1996, 121, 207−225. (21) García-Sánchez, F.; Eliosa-Jiménez, G.; Salas-Padrón, A.; Hernández-Garduza, O.; Á pam-Martínez, D. Modeling of Microemulsion Phase Diagrams from Excess Gibbs Energy Models. Chem. Eng. J. 2001, 84, 257−274. (22) Stateva, R. P. Predicting and Calculating Complex Phase Equilibrium in Supercritical Fluid Systems with a New Technique. Collect. Czech. Chem. Commun. 1995, 60, 188−210. (23) Tsvetkov, G., St.; Stateva, R. P. A Computationally Efficient Algorithm for Simultaneous Phase and Chemical Equilibrium Calculations. Collect. Czech. Chem. Commun. 1997, 62, 558−574. (24) Stateva, R. P.; Cholakov, G. S.; Galushko, A. A.; Wakeham, W. A. A Powerful Algorithm for Liquid−Liquid−Liquid Equilibria Predictions and Calculations. Chem. Eng. Sci. 2004, 55, 2121−2129.

fractions predicted with both thermodynamic models do not agree well with the experimental data; that is, the SRK underpredicts the experimental data whereas the PC-SAFT overpredicted them. In fact, Figure 13 shows that the PC-SAFT overpredicts both liquid phases, which indicates that this thermodynamic model tends to predict better the experimental LLV behavior as temperature increases. This could be a result of the fact that the interaction parameter of the nitrogen + nbutane binary system for this model was determined at temperatures higher than those studied here. Figure 13 also shows that the UCST calculated with the SRK (3.295 MPa) is only slightly closer to the “hypothetical” experimental UCST (3.36 MPa) than that estimated with the PC-SAFT (3.44 MPa), whereas the calculated K-point with the SRK (4.72 MPa) is very similar to that estimated with the PCSAFT (3.69 MPa).



CONCLUSIONS The prediction and calculation of the phase behavior, exhibited at low temperatures, by three nitrogen-containing systems (nitrogen + methane + ethane, nitrogen + methane + propane, and nitrogen + methane + n-butane), of interest to the natural gas and oil industry, with two numerical techniques applying the PC-SAFT and SRK EoSs, was investigated in this work. The results demonstrate that both thermodynamic models perform well and predict with reasonable accuracy the experimentally observed phase behavior of the systems considered. This study complements the results involving the prediction of the threephase VLL behavior for other three ternary systems (nitrogen + methane + n-pentane, nitrogen + methane + n-hexane, and nitrogen + methane + n-heptane) studied in a previous work by the authors. Furthermore, a procedure similar to that presented by Gregorowicz and de Loos, was advocated and applied to predict K- and/or LCST (UCST for the nitrogen + methane + nbutane) points for the systems studied at constant temperature. The thermodynamic model used was the SRK EoS. The procedure was also utilized to predict the K- and LCST-point loci of the nitrogen + methane + ethane system up to 149.5 K; a temperature close to the tricritical point and in all cases studied demonstrated a reliable and robust performance with good convergence characteristics. Therefore, it can be used with confidence to predict the K-point and LCST of ternary systems using any thermodynamic model, as will be demonstrated in the future for PC-SAFT EoS.



AUTHOR INFORMATION

Corresponding Author

*Tel.: +52 55 9175 6574. E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was partially supported by the Mexican Petroleum Institute under Project D.00406. One of the authors (D. N. J.G.) gratefully acknowledges the National Polytechnic Institute of Mexico for their financial support through the Project SIP20110150.



REFERENCES

(1) Llave, F. M.; Luks, K. D.; Kohn, J. P. Three-Phase Liquid− Liquid−Vapor Equilibria in the Nitrogen + Methane + Ethane and 9417

dx.doi.org/10.1021/ie300372a | Ind. Eng. Chem. Res. 2012, 51, 9409−9418

Industrial & Engineering Chemistry Research

Article

(25) Prausnitz, J. M.; Lichtenthaler, R. N.; Gomes de Acevedo, E. Molecular Thermodynamics of Fluid-Phase Equilibria, 3rd ed.; Prentice Hall PTR: Upper Saddle River, NJ, 1999. (26) Mushrif, S. H.; Phoenix, A. V. An Algorithm to Calculate K- and L-Points. Ind. Eng. Chem. Res. 2006, 45, 9161−9170. (27) Heidemann, R. A.; Khalil, A. M. The Calculation of Critical Points. AIChE J. 1980, 26, 769−779. (28) Baker, L. E.; Luks, K. D. Critical Point and Saturation Pressure Calculations for Multicomponent Systems. Soc. Pet. Eng. J. 1980, 20, 15−24. (29) Topliss, R. J.; Dimitrelis, D.; Prausnitz, J. M. Computational Aspects of a Noncubic Equation of State for Phase Equilibrium Calculations. Effect of Density-Dependent Mixing Rules. Comput. Chem. Eng. 1988, 12, 483−489. (30) Topliss, R. J. Techniques to Facilitate the Use of Equations of State for Complex Fluid-Phase Equilibria. Ph.D. Dissertation, University of California, Berkeley, CA, 1985. (31) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in FORTRAN: The Art of scientific Computing; Cambridge University Press: New York, 1992. (32) Knapp, H. R.; Döring, R.; Oellrich, L.; Plöcker, U.; Prausnitz, J. M. Vapor-Liquid Equilibria for Mixtures of Low Boiling Substances; DECHEMA Chemistry Data Series: Frankfurt, Germany, 1982; Vol. VI. (33) Nagy, Z.; Shirkovskiy, I. A. Mathematical Simulation of Natural Gas Condensation Processes Using the Peng-Robinson Equation of State. Paper SPE 10982 presented at the 57th Annual Fall Technical Conference and Exhibition of the Society of Petroleum Engineers of AIME, New Orleans, LA, September 26−29, 1982. (34) García-Sánchez, F.; Eliosa-Jiménez, G.; Silva-Oliver, G.; Vázquez-Román, R. Vapor−Liquid Equilibria of Nitrogen-Hydrocarbon Systems Using the PC-SAFT Equation of State. Fluid Phase Equilib. 2004, 217, 241−253. (35) Justo-García, D. N.; García-Sánchez, F.; Díaz-Ramírez, N. L.; Romero-Martínez, A. Calculation of Critical Points for Multicomponent Mixtures Containing Hydrocarbon and Nonhydrocarbon Components with the PC-SAFT Equation of State. Fluid Phase Equilib. 2008, 265, 192−204. (36) Rowley, R. L.; Wilding, W. V.; Oscarson, J. L.; Yang, Y.; Zundel, N. A. DIPPR Data Compilation of Pure Chemical Properties: Design Institute for Physical Properties; Brigham Young University: Provo, UT, 2006. http//dippr.byu.edu. (37) Llave, F. M.; Luks, K. D.; Kohn, J. P. Three-Phase Liquid− Liquid−Vapor Equilibria in the Binary Systems Nitrogen + Ethane and Nitrogen + Propane. J. Chem. Eng. Data 1985, 30, 435−438. (38) Konynenburg, P. H. Critical Lines and Phase Equilibria in Binary Mixtures. Ph.D. Dissertation, University of California: Los Angeles, CA, 1968.

9418

dx.doi.org/10.1021/ie300372a | Ind. Eng. Chem. Res. 2012, 51, 9409−9418