Ind. Eng. Chem. Res. 1999, 38, 1107-1113
1107
Phase-Boundary Calculations in Systems Involving More than Two Phases, with Application to Hydrocarbon Mixtures Niels Lindeloff,† Simon I. Andersen, and Erling H. Stenby Department of Chemical Engineering, Building 229, Technical University of Denmark, DK-2800 Lyngby, Denmark
Robert A. Heidemann* Department of Chemical and Petroleum Engineering, University of Calgary, 2500 University Drive, N.W., T2N 1N4 Calgary, AB Canada
An algorithm to track phase boundaries in systems with more than two phases is presented. The method follows the basic principles for phase diagram construction developed by Michelsen. The method is applied successfully to some normal alkane systems with solid-liquid-vapor (SLV) phase behavior. Some thermodynamic aspects of correlating the phase-boundary data for these systems are discussed, particularly the influence of equation of state parameters on the solid-liquid and SLV boundary calculation. Introduction Precipitation of solid waxes is a common problem in connection with the production and transport of hydrocarbons in cold environments. Most crude oils and even some gas condensates have an amount of heavy paraffinic molecules dissolved under reservoir conditions, which may precipitate as a result of cooling and/or depressurization in connection with production and pipeline transport. As drilling and production technologies have advanced, deep high-pressure offshore reservoirs have become subject to exploitation, which has led to a renewed interest in the effects of pressure on the precipitation of solid waxes. A way of getting an understanding of operating conditions leading to the precipitation of solid paraffins is to construct the pressure-temperature phase diagram for the mixture under consideration. The principles for calculation of vaporliquid-phase boundaries were outlined by Michelsen.1 Nghiem and Li2 applied the same general ideas to slightly more complex systems exhibiting three-phase (liquid-liquid-vapor) separation. In the present work, an algorithm following principles similar to those for calculating multiphase boundaries is applied for calculation of two- and three-phase boundaries for model mixtures consisting of methane, n-decane, and solidforming n-alkanes in the range of C18-C30. Daridon et al.3 have presented experimental data for several of these model systems. The algorithm is presented in a general form to track the boundary where any new phase may form from a multiphase equilibrium. The method includes the possibility of the new phase being a mixed solid. Formulation of the Algorithm In the following, the equations that may be used to determine the phase boundaries in a system of nc components in nph phases will be described. The problem is formulated as a set of equilibrium and † Fax: +45 45882258. E-mail:
[email protected]. * To whom correspondence should be addressed. Telephone: (403)220-8755.Fax: (403)284-4852.E-mail: heideman@ ucalgary.ca.
specification equations applying to nc components in nph phases. This set of algebraic equations is solved iteratively with respect to the independent variables by a Newtonian method, following the procedures outlined by Michelsen.1 The array of independent variables R contains the logarithm of the equilibrium constants, the phase fractions, and the logarithm of temperature and pressure: R ) [ln Kjn, ..., βn, ... ln T, ln P]. In the formulation of the equations, subscripts i, j, and k will be used as component identifiers, and superscripts l, m, and n will be used as phase identifiers. Superscript R identifies the reference phase. The equations in the form of an array of objective functions g are formulated in the following way: m nc(nph - 1) equilibria: gj ) ln Km i + ln φ i -
ln φ Ri ) 0, j ) nc(m - 1) + i (1) (nph - 1) summation equations: R gj ) (Xm k - Xk ) ) 0, j ) nc(nph - 1) + m (2)
∑k
1 constraint on the incipient phase: gj ) βincip ) 0, j ) nc(nph - 1) + nph (3) 1 phase fractions constraint: gj )
∑l βl - 1 ) 0,
j ) nc(nph - 1) + nph + 1 (4) 1 specification: gj ) Ri - S ) 0, j ) nc(nph - 1) + nph + 2 (5) The mole numbers X in eq 2 are defined from the mass balances: m Xm k ) zkKk /
∑l βlKlk
There are nc(nph - 1) + nph + 2 g equations in the nc(nph - 1) + nph + 2 unknowns r. Derivatives of g with respect to the independent variables r can be obtained to form a Jacobian J ) ∂g/∂r. The derivatives are presented in Appendix A. The set of equations can now be solved iteratively for the set of independent variables that satisfies the system of equations through a Newtonian scheme. The calculations are initialized in the first point using
10.1021/ie980415h CCC: $18.00 © 1999 American Chemical Society Published on Web 01/22/1999
1108
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
compositions from an isothermal flash calculation near the phase boundary. The general form of the iterative scheme is as follows:
J(k)∆(k) + g(r)(k) ) 0, r(k+1) ) r(k) + ∆(k)
(6)
The system is solved with respect to the correction vector ∆ by lower-upper decomposition and back substitution according to Crout’s method. Michelsen pointed out that, once a converged solution has been obtained, the sensitivity of a given solution r* to changes in the parameters, ∂R/∂RS; can be obtained very easily by solving J ∂R/∂RS + ∂R/∂RS ) 0. This only requires an additional call to the back-substitution subroutine. Thermodynamics The algorithm was applied to n-alkane systems forming solid, liquid, and vapor phases. In the following, the thermodynamic models used are discussed. Fugacity coefficients for the vapor and liquid phases were calculated using the Soave-Redlich-Kwong (SRK) equation of state with quadratic mixing rules. For evaluation of the fugacity coefficients, the Xk are treated formally as mole numbers. The derivatives of fugacity coefficients with respect to mole numbers pressure, and temperature required to form the Jacobian were obtained following the principles outlined by Michelsen and Mollerup.4 Becaise no experimental data for critical properties of the heavy normal alkanes are available, these had to be estimated using correlations. Various alternatives are available, and it was found that the choice of correlation has a profound effect on the model predictions. The correlations by Cavett,5 Twu;6 and Margoulas and Tassios7 were investigated. The Cavett expression yields average properties for crude oil cuts and thus should not be the optimum choice for mixtures of normal alkanes. The Twu correlations yield critical properties from perturbations with n-alkanes as the reference system. The properties for n-alkanes obtained with the Twu expressions are virtually identical to those obtained with the correlations by Margoulas and Tassios, but the former were preferred, because they may be applied to other systems as well. The acentric factors were in all cases estimated with the Edmister correlation. It was found that the parameters obtained with the Twu expressions gave the best representation of the shape of the experimental bubble point curves, but it was necessary to apply binary interaction coefficients in order to get a good match of the bubble point curves. The interaction coefficients were estimated using the expression by Chueh and Prausnitz:8
kij ) 1 -
(
1/6 2ν1/6 ci νcj 1/3 ν1/3 ci + νcj
)
Figure 1. Influence of critical parameters on predicted bubble points with the SRK equation. Daridon et al.3 mixture B: 43.8% methane, 45.9% n-decane, and 10.3% of the normal alkanes C18C30.
Figure 2. Influence of critical parameters on predicted bubble points with the SRK equation. Daridon et al.3 mixture C: 43.6% methane, 46.15% n-decane, and 10.25% of the normal alkanes C18-C30.
single mixed-solid phase with orthorhombic crystal structure. Electron diffraction studies9 have suggested that the orthorhombic structure is the preferred state in multicomponent paraffin waxes. For the solid phase, the fugacity coefficients were obtained in the following way. In general, the differential of the fugacity f of a pure component at temperature T and pressure P is
(
d ln f °i ) -
(7)
Figures 1 and 2 give examples of how the different approaches perform for a multicomponent mixture of normal alkanes and methane. The curves are calculated with our algorithm, and a comparison is made with the data of Daridon et al.3 Waxes may precipitate as multiple pure phases or as mixed-solid solutions. Paraffins that differ significantly in size, or with even and odd numbers of carbon atoms, tend to form separate phases. Mixing in the solid is possible with normal alkanes of similar sizes. In the present work, it is assumed that the paraffins form a
)
H°i - Hig i RT
2
dT +
νi dP RT
(8)
where H°l is the enthalpy of the pure substance, Hig l is the ideal gas enthalpy, and vi is the molar volume. This equation applies to any given phase. For the ratio of pure component fugacities in two phases, in this case a solid and a liquid, we obtain
d ln
f °iS
)f °iL
HSi - HLi RT2
dT +
νSi - νLi dP ) RT ∆Hfi RT
2
dT -
∆νi dP (9) RT
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 1109
We now do the line integration from the fusion temperature Tf to T and from the zero pressure limit to P. Under the assumption that the heat of fusion is a constant, the result for the ratio of the pure component fugacities at T and P is:
ln
(
f °iS
)
∆Hfi T ∆νP ) 1- f L RT RT f °i Ti
(10)
In this expression, we may use the approximation that the melting temperature is close to transition temperatures between different crystal structures and that heats of melting and transition may simply be added in cases where a low-temperature crystal form is considered. The refinement of adding a temperature variation in the heat of fusion would add very little to the correlating power of the model. Equation 10 gives the fugacity for the pure solid in terms of the fugacity of the pure liquid. The fugacity of a solid component in a mixture is related to that of the pure component through
f Si ) xSi γSi f °iS
Figure 3. Influence of critical parameters on predicted wax appearance points with the SRK equation. Daridon et al.3 mixture A: 43.7% methane, 46.1% n-decane, and 10.2% of the normal alkanes C20-C30.
(11)
By the introduction of fugacity coefficients (φi ) fi/xiP) and the combination of eqs 10 and 11, an expression for the solid fugacity coefficient is obtained:
ln φSi ) ln γSi + ln φ°iL -
(
)
∆Hfi ∆νiP T 1- f RT RT T i
(12)
As we use this equation, the fugacity coefficient for the pure liquid, φi°L is calculated from the equation of state for the pure substance at the mixture T and P. The equilibrium ratios for components in the solid phase are given by
ln KSi ) ln φRi - ln φSi
(13)
φRi
is the fugacity component for component i in where the reference phase, which would be one of the phases calculated from the equation of state model. For mixtures of normal alkanes with orthorhombic crystal structure, Coutinho and Stenby10 found that an adapted formulation of the Wilson equation gave a good representation of the solid-phase behavior. The Wilson activity-coefficient expression is
ln γi ) 1 - ln
∑j xjΛij - ∑k xkΛki/∑j xjΛkj
(14)
where the energy parameters Λij are given by
Λij )
(
)
νj λij - λji exp νi RT
(15)
In this expression, vi represents the molar volume. Coutinho and Stenby10 suggest that the λii energy parameters can be obtained a priori as
2 λii ) - (∆Hsubl - RT) Z
(16)
where the enthalpic term is the isothermal heat of
Figure 4. Influence of critical parameters on predicted wax appearance points with the SRK equation. Daridon et al.3 mixture D; 44% methane, 45.8% n-decane, 6.8% n-C18, and 3.4% n-C30.
sublimation. The heat of sublimation is calculated as the sum of the heat of fusion and the heat of vaporization, the latter being estimated using the correlations by Morgan and Kobayashi.11 The coordination number Z is 6 because only axial interactions are considered. For the λii cross energies, it is assumed that the interaction energy between two rods of different lengths placed alongside each other is characterized by that of the shorter molecule, because the contact area will be controlled by the shorter length. Thus, λij ) λii in the case where molecule i is shorter than molecule j, else λij ) λjj. The Coutinho and Stenby10 model was used in the calculations reported in this paper. In the calculations, methane and decane were essentially excluded from the solid phase. This was done by assigning a very high fugacity coefficient to these components in the solid. The derivatives of the solid fugacity coefficients needed to form the Jacobian are listed in Appendix B. Brief mention will be made of an alternative approach to the modeling of wax precipitation that was proposed by Won.12 In Won’s method, the liquid-vapor equilibrium is calculated from the equation of state, as is also done by us. The solid equilibrium ratio, however, is
1110
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
Table 1. Reduction in Hf and Tf (Keeping the Ratio Fixed) To Match Experimental Wax Appearance Temperatures mixture A mixture B mixture C mixture D overprediction (K) without tuning tuning (%)
7.2
7.2
5.0
4.7
2.3
2.4
1.7
1.3
calculated from
KSi ) xSi /xLi ) γLi f °iL/γSi f °iS
(17)
The ratio of the pure-liquid fugacity/pure-solid fugacity is calculated from an equation like eq 10 but with the terms included that account for a temperature-dependent heat of fusion. Won also evaluates activity coefficients for the solid wax-forming components in both the liquid phase and the solid phase from regular solution models with different solubility parameters in the two phases. Although a fugacity is calculated for a component in the liquid phase from the equation of state, it is not used in the solid-fluid equilibrium equations. In the method described in this paper, no separate model for liquid activity coefficients is super imposed on the equation of state calculations. Rather, we use the equation of state calculations to find the pure-liquid fugacity needed in eq 12. Discussion of Model Behavior The binary interaction coefficients in the mixing rule for the a parameter in the equation of state affect the temperature dependence of the vapor and liquid fugacities, and consequently, they will affect the model’s ability to predict the solid appearance line as well. Figures 3 and 4 demonstrate the model performance with and without the binary interaction coefficients. All lines are calculated using the proposed algorithm. The
figures show that the parameters obtained with the Twu correlations and without interaction coefficients can produce a reasonably good match of the solid appearance line in some instances, apart from the fact that the shape of the curve is wrong because of a poor prediction of the bubble points. The use of parameters from the Twu correlations and binary interaction coefficients fitted to the bubble point curve results in toohigh solid appearance temperatures. In conclusion, both phase boundaries cannot be modeled correctly without additional adjustment of the model. It is possible to compensate by tuning the heats of fusion and/or temperatures of fusion in eq 12. In doing so, it must be recognized that one is adjusting solid-model parameters to compensate for unphysical behavior of the liquid and vapor model contained in the equation of state. To simultaneously get a good correlation of solidliquid (SL), solid-liquid-vapor (SLV), and liquid-vapor (LV) phase boundaries for the n-alkane mixtures discussed, the critical parameters were estimated with the Twu correlations and binary interaction coefficients were applied. The heats of fusion and fusion temperatures for the higher alkanes were taken from the experimental work by Broadhurst.13 The ratio H f/Tf (i.e., the entropy of fusion) was fixed at the experimental values, and Hf and Tf were tuned to match the experimental solid appearance temperatures by multiplying all values by a single constant factor, generally resulting in a reduction of the fusion values by 2-3%. Table 1 summarizes the over prediction of the wax precipitation temperature without tuning and the percent reduction in the heats of fusion required to match the data of Daridon et al.3 The effect of including the Poynting-type correction term in eq 12 was found to be quite significant.14 A potential problem in that connection is to estimate the solid molar volumes. This was done by exploiting the experimental observation that normal alkanes shrink by approximately 10% upon solidification.15,16 The liquid
Figure 5. Calculation of wax appearance temperatures with and without the Poynting term. Daridon et al.3 mixture B.
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 1111
point curve. Without the solid, the algorithm goes from using two to just one thermodynamic model, and consequently, the derivatives ∂r/∂rS are more “well behaved”. Conclusions
Figure 6. Entire phase envelope showing the high pressure solid boundary. Daridon et al.3 mixture C.
In the present work, it was demonstrated that the phase-boundary calculation method by Michelsen could be extended successfully to systems involving more than two phases. The algorithm was applied to some n-alkane systems that exhibit solid wax precipitation at low temperatures. It was discussed how the choice of estimation method for the equation of state parameters significantly affects the model’s ability to predict correctly wax appearance temperatures. The results demonstrate that caution must be exerted when cubic equations of state are applied in hybrid models, as has previously been pointed out by Lundgaard and Mollerup.18 Another point is that some of the normal alkanes modeled in the present work are quite large molecules, certainly larger than those most cubic equations were originally designed to deal with. It was shown that experimental data could be correlated through adjustment of the fusion parameters in the solid model. Acknowledgment This research was supported, in part, by a grant from the Natural Sciences and Engineering Research Council of Canada. Nomenclature
Figure 7. Phase envelope of Figure 6 with emphasis on the vaporphase boundary.
molar volumes were estimated by a group contribution method17 and were treated as being pressure independent. Figure 5 shows an example of a calculation with and without the Poynting term in eq 12. The entire phase diagram for the systems can be constructed with the algorithm. The entire LV boundary may be calculated, including the critical point and the dew point line. In doing so, the curves are followed well beyond the cracking temperatures of the molecules, meaning that the high-temperature part of the curve is unphysical. However, this is the simplest way to construct the dew point line. Figures 6 and 7 contain an example of the entire phase diagram for one of the Daridon et al.3 mixtures. Figure 7 focuses on the vapor-liquid coexistence region. From a computational point of view, the distinct “kink” on the vapor boundary in the region where it crosses the solid line represents the greatest challenge to the algorithm. As the temperature is raised and the solid disappears, the system of equations must be reduced. If all of the equations are retained, the disappearance of a phase results in the corresponding phase fraction becoming negative. A small step size is required in the region where a phase boundary is crossed in order to get convergence. This is particularly the case when tracking the bubble
a ) equation of state parameter f ) fugacity g ) function to be zeroed at the solution g ) vector of functions to be zeroed H ) enthalpy J ) Jacobian matrix kij ) interaction parameter K ) equilibrium ratio m ) phase number nc ) number of components nph ) number of phases P ) pressure R ) universal gas constant S ) value of the specified variable T ) temperature v ) volume x ) mole fraction z ) feed mole fraction Z ) coordination number Greek Letters r ) vector of independent variables β ) phase amount ∆v ) volume difference ∆ ) vector, Newtonian iteration step γ ) activity coefficient λ, Λ ) Wilson equation parameters φ ) fugacity coefficient Subscripts c ) critical i, j, k ) component indices (k) ) iteration count j ) equation index S ) specified variable subl ) sublimation
1112
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999
Superscripts
∂gj
f ) fusion ig ) ideal gas l, m, incip ) phase indices L ) liquid ° ) pure substance R ) the reference phase S ) solid
∂βn
j ) nc(nph - 1) + nph + 1 ∂gj
∂ ln
Knj
) δijδnm + Xnj
nc
)
∂βn
ln φm i m ∂Xj
∂
∑k z
k
XRk
n
δnm -
∂ ln φRi
β
]
Xnj
∂XRk
∂βn
zj
Xm k
+
[ ]}
]
∂ ln φm i ∂Xm k
(
)
(
)
)1
Specification constraint gj ) ri - S ) 0, j ) nc(nph - 1) + nph + 2
∂ ln φRi βnXRj zj ∂XR
-
∂gj
j ) nc(m - 1) + i
{ [
[
Xnk
∂gj ∂gj ) )0 ∂ ln T ∂ ln P
m gj ) ln Km i + ln φi -
j
∂gj
)
Knj
∂ ln
ln φRi ) 0, ∂gj
∑l βl - 1 ) 0,
Phase fraction constraint gj )
Appendix A. Derivatives of g
Equilibrium equations
) δn,incip
∂gj ) 1, all others ) 0 ∂rS Appendix B. Derivatives of the Solid “Fugacity Coefficient”
∂ ln φSi
)
∂XSj
∂ ln γSi ∂XSj
XkΛkiΛkj
nc
1
)
+
nc
∑l Xl
∑k
∑l XlΛkl)2
(
(
∂gj ∂ ln φm ∂ ln φRi i )T ∂ ln T ∂T ∂T ∂ ln φRi ∂gj ∂ ln φm i )P ∂ ln P ∂P ∂P Mole number equations gj )
∑k (Xmk - XRk ) ) 0, j ) nc(nph - 1) + m
∂gj ∂ ln Knj
Xnj
)
∂gj ∂β
([
] [ ]}
βnXm βnXRj j δnm + zj zj
nc
)
n
Xnk
∑k z
(XRk - Xm k)
k
∂gj ∂gj ) )0 ∂ ln T ∂ ln P Incipient phase constraint gj ) βincip ) 0, j ) nc(nph - 1) + nph ∂gj ∂ ln
Knj
)
∂gj ∂gj ) )0 ∂ ln T ∂ ln P
-
nc
Λij nc
)
Λji
+ nc
∑l XlΛil ∑l XlΛjl
φSi
∂ ln ∂T
)
γSi
∂ ln ∂T nc
∂ ln γSi ∂T
)-
∂νSi P ∂ ln φ°i ∂T + + + ∂T RT RT RT
∂Λil
∑l Xl ∂T nc
∑l
L
XlΛil
Literature Cited
nc
-
{
∑k Xk
∆Hfi 2
νSi P 2
nc
∂Λki ∂T nc
∑l
-
nc
XlΛkl
∑l XlΛkl)2
(
}
∂Λkl
∑l Xl ∂T
Λki
(1) Michelsen, M. L. Calculation of Phase Envelopes and Critical Points for Multicomponent Mixtures. Fluid Phase Equilib. 1980, 4, 1. (2) Nghiem, L. X.; Li, Y.-K. Computation of Multiphase Equilibrium Phenomena with an Equation of State. Fluid Phase Equilib. 1984, 17, 77. (3) Daridon, J.-L.; Xans, P.; Montel, F. Phase Boundary Measurement on a Methane + Decane + Multiparaffins System. Fluid Phase Equilib. 1996, 117, 241. (4) Michelsen M. L.; Mollerup, J. Partial Derivatives of Thermodynamic Properties. AIChE J. 1986, 32, 1389. (5) Pedersen, K. S.; Fredenslund, Aa.; Thomassen, P. Properties of Oils and Natural Gases; Gulf Publishing: Houston, 1989. (6) Twu, C. H. An Internally Consistent Correlation for Predicting the Critical Properties and Molecular Weights of Petroleum and Coal-Tar Liquids. Fluid Phase Equilibria 1984, 16, 137. (7) Margoulas, K.; Tassios, D. Thermophysical Properties of n-Alkanes from C1 to C20 and Their Prediction for Higher Ones. Fluid Phase Equilibria 1990, 56, 119. (8) Chueh, P. L.; Prausnitz, J. M. Vapor-liquid equilibria at high pressures: Calculation of partial molar volumes in nonpolar liquid mixtures. AIChE J. 1967, 13, 1099.
Ind. Eng. Chem. Res., Vol. 38, No. 3, 1999 1113 (9) Dorset, D. L. The Crystal Structure of Waxes. Acta Crystallogr. 1995, B51, 1021. (10) Coutinho, J. A. P.; Stenby, E. H. Predictive Local Composition Models for Solid/Liquid Equilibrium in n-Alkane Systems: Wilson Equation for Multicomponent Systems. Ind. Eng. Chem. Res. 1996, 35, 918-925. (11) Morgan, D. L.; Kobayashi, R. Extension of Pitzer CSP Models for Vapor Pressures and Heats of Vaporization to Long Chain Hydrocarbons. Fluid Phase Equilibria 1994, 94, 51. (12) Won, K. W. Thermodynamics for Solid Solution-LiquidVapor Equilibria: Wax Phase Formation from Heavy Hydrocarbon Mixtures. Fluid Phase Equilibria 1986, 30, 265. (13) Broadhurst, M. G. An Analysis of the Solid-Phase Behavior of the Normal Paraffins. J. Res. Natl. Bur. Stand., Sect. A 1962, 66, 241. (14) Lindeloff, N. Formation of Solid Phases in Hydrocarbon Mixtures. M.Sc. Thesis, Technical University of Denmark, Lyngby, Denmark, 1995.
(15) Schaerer, A. A.; Busso, C. J.; Smith, A. E.; Skinner, L. B. Properties of Pure Normal Alkanes in the C17-C30 Range. J. Am. Chem. Soc. 1955, 77, 2017. (16) Templin, P. R. Coefficient of Volume Expansion for Petroleum Waxes and Pure n-Paraffins. Ind. Eng. Chem. 1956, 48, 154. (17) Elbro, H. S.; Fredenslund, Aa.; Rasmussen, P. Group Contribution Method for the Prediction of Liquid Densities as a Function of Temperature for Solvents, Oligomers, and Polymers. Ind. Eng. Chem. Res. 1991, 30, 2576. (18) Lundgaard, L.; Mollerup, J. The Influence of Gas-Phase Fugacity and Solubility on Correlation of Gas Hydrate Formation Pressure. Fluid Phase Equilibria 1991, 70, 199.
Received for review June 29, 1998 Revised manuscript received November 9, 1998 Accepted November 13, 1998 IE980415H
