14 Phase Behavior Prediction and Modeling of LNG Systems with EoSs – What is Easy and What is Difficult ?

Blanca E. García-Flores1, Daimler N. Justo-García2, Roumiana P. Stateva3 and Fernando García-Sánchez1 1Laboratory of Thermodynamics, Research Program in Molecular Engineering, Mexican Petroleum Institute, Mexico, D.F., 2Department of Chemical and Petroleum Engineering, ESIQIE, National Polytechnic Institute, Mexico, D.F., 3Institute of Chemical Engineering, Bulgarian Academy of Sciences, Sofia, 1,2México, 3Bulgaria


Introduction
In the past decades, natural gas processing has increasingly encompassed low temperatures with the view to recover ethane and heavy components from natural gas.Liquefied natural gas (LNG) is a useful method for storing the gas in a small space at peak-shaving facilities of natural gas distribution companies.The increased density of LNG facilities allows transportation of natural gas via large ocean vessels from gas fields situated far from their potential markets.Also, in the past, natural gas has been an important source of chemical feed stocks like ethane, propane, butane, etc.In order to increase the recovery of these stocks, the gas processes have been shifting to lower temperatures where the recoveries are improved.In all of the above processes and applications of natural gas, knowledge of the LNG systems phase behavior and thermodynamic properties is required for the successful design and operation of LNG plants.
The question of how to predict and model phase equilibria behavior for natural gas systems is far from new and liquid-vapor equilibria problems have been successfully solved for many years now.At present, interest has shifted to systems containing not only species in the simple paraffinic homologous series, but also water, carbon dioxide, hydrogen sulfide, hydrogen, nitrogen to mention a few.Not unrelated to this growing variety of species in cryogenic process streams is the occurrence of multiphase phenomena.The situation of interest in natural gas processing is often a methane-rich stream for which there is the question whether a second solvent can alter dramatically the pattern of the phase behavior of the customarily liquid-vapor mixture and cause problems.
The active 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.However, in order to help further understanding the possible occurrence of multiphase equilibria in LNG process systems, it is necessary to acquire knowledge of their phase behavior and of the variety of critical end point boundaries through an ability to predict, model and calculate them.
There are several aspects of separation/refinement of natural and synthetic gases where multiphase equilibria play a role.The occurrence of liquid-liquid-vapor (LLV) behavior during the recovery of natural gas by low temperature distillation, especially from nitrogenrich gas mixtures, is just one such example.LLV behavior can also occur in the processing of gases containing high quantities of carbon monoxide that result from coal gasification, as well as in high pressure absorption processes for the removal of either desirable or undesirable components from natural and synthetic gas mixtures.In these latter processes, methanol can be used as the absorber to separate out feed stocks (ethane, ethylene, propane, propylene, etc.), harmful components such as hydrogen sulfide, carbon disulfide, and carbonyl sulfide, or simply unwanted species such as carbon dioxide.Obviously, the formation of a second liquid phase can upset the expected performance of these processes.
In the cryogenic processing of natural gas mixtures, species such as CO 2 and the heavier hydrocarbons can form solids and foul the gas processing equipment.For example, the formation of a solid phase can coat heat exchangers and foul expansion devices, leading to process shutdown and/or costly repairs.Knowledge of the precipitation conditions for gas streams is essential in minimizing downtime for cleanup and repairs.However, the appearance of a solid phase in a process is not always a liability.Off-gas (primary nitrogen) from power plants with light water reactors or from fuel processing plants can contain radioactive isotopes of krypton and xenon, which may be removed by solid precipitation.The phenomenological aspects of LLV/solid-liquid-vapor (SLV) behavior are also interesting because there is a need for a better understanding of the physical nature of the thermodynamic phase space, especially the non-analyticities of critical point region.
The success of the design and operation of separation processes in the oil and gas industry at low temperatures is critically dependent on accurate descriptions of the thermodynamic properties and phase behavior of the concerned multicomponent hydrocarbon mixtures with inorganic gases.Consequently, it is important to apply appropriate models within a thermodynamic modeling framework to predict, describe and validate the complex phase behavior of LNG mixtures.In this case, equations of state (EoS) are usually the primary choice.
The aim in this chapter is to examine and analyze the challenges and difficulties encountered when modeling the complex phase behavior of LNG systems, and to compare the capabilities of two numerical techniques advocated for phase behavior predictions and calculations of complex multicomponent systems.The chapter is organized as follows.In Section 2 the thermodynamics and topography of the phase behavior of LNG systems are presented.Section 3 describes two computational techniques, advocated by us, to predict and model complex phase behavior of nonideal mixtures, and their application to LNG systems.Section 4 is focused to the description of two numerical methods to directly calculate critical end points (K-and L-points).In Section 5 two representatives of the 361 equations of state (EoSs) type thermodynamic models, namely SRK EoS (Soave, 1972) and PC-SAFT EoS (Gross and Sadowski, 2001) are used to represent the equilibrium fluid phases.Results and discussion of the multiphase behavior modeling for selected ternary systems studied are presented in Section 6.Finally, conclusions derived from the work are given in Section 7.

Thermodynamics and topography of the phase behavior of LNG systems
Though only a limited number of immiscible binary systems (methane-n-hexane, methanen-heptane, to name the most prominent ones) are relevant to natural gas processing, LLV behavior can and does occur under certain conditions in ternary and higher LNG systems even when none of the constituent binaries themselves exhibit such 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 qualitative classification of natural gas systems and the topography of the multiphase equilibrium behavior of the systems in the thermodynamic phase space and the nature of the phase boundaries will be addressed in this section.Kohn and Luks (1976, 1977, 1978), who carried out extensive experimental studies on the solubility of hydrocarbons in LNG and NGL (natural gas liquids) mixtures, qualitatively classify natural gas systems (or any system) as one of four types (Kohn and Luks, 1976).These types of binary systems show, for instance, that a first type system, methane-n-octane (Kohn and Bradish, 1964), has a solid-liquid-vapor (S-L-V) locus which starts at the triple point of the solute and terminates at a Q-point (S 1 -S 2 -L-V) near the triple point of the solvent with an S-V gap in the locus.The S-L-V branches terminate with a K-point where the liquid becomes critical with vapor in the presence of a solid.On the contrary, methane-n-heptane (Kohn, 1961), a second type system, has a Q-point (S-L 1 -L 2 -V) in the central portion of the SLV locus, at which point one loses the L 1 phase and gains the L 2 phase (solvent lean phase).These systems also have a L 1 -L 2 -V locus running from the Q-point to a K-point where the liquid L 2 becomes critical with vapor in the presence of L 1 .The other two types of systems have no discontinuities between the triple points of the solvents; however, one of these types, methane-n-hexane (Shim and Kohn, 1962), a typical third type system, has a L 1 -L 2 -V locus which is terminated by critical points, points where the two liquid phases, or a liquid and vapor, become critically identical.A typical system of the fourth type is ethane-n-octane (Kohn et al., 1976).The topographical evolution of multiphase equilibria behavior is discussed in detail by Luks and Kohn (1978).Ternary ( 3 N  ) or more complex systems exhibit similar types of phenomena, but the loci and exact points of a binary become 1 N  and 2 N  dimensional spaces in the more complex systems.
It is known that systems rich in methane can exhibit L 1 -L 2 -V behavior of the second and third type variety (e.g., methane-n-heptane and methane-n-hexane, respectively).However, investigations into ternary S-L-V and more complex systems revealed that L 1 -L 2 -V behavior can occur in systems whose binary pairs exhibit no immiscibility (Green et al., 1976;Orozco et al., 1977).Hottovy et al. (1981Hottovy et al. ( , 1982) ) observed that systems of the first type (methane-noctane) could be modified to behave like a second type system by adding a second solvent (e.g., methane-ethane-n-octane, methane-propane-n-octane, methane-n-butane-n-octane, and methane-carbon dioxide-n-octane).Merrill et al. (1983) reported the phase behavior of the ternary systems methane-n-pentane-n-octane and methane-n-hexane-n-octane, which also exhibit L 1 -L 2 -V immiscibility.Additionally, these authors studied the systems methane-nhexane-carbon dioxide, in which the carbon dioxide is added to the pair methane-n-hexane (third type) to induce L 1 -L 2 -V immiscibility.
In five of these seven ternary systems, methane-n-octane-(ethane, propane, n-butane, npentane, and carbon dioxide), the immiscibility region is bounded by loci of type K, type Q and LCST (lower critical solution temperature) points, whereas for the system methane-nhexane-n-octane the region of L 1 -L 2 -V immiscibility is bounded, apart from the K, Q, and LCST points, by the L 1 -L 2 -V locus of the methane-n-hexane binary system.In the case of the system methane-n-hexane-carbon dioxide, the immiscibility region is bounded by a locus of K and LCST points and by the L 1 -L 2 -V locus of the methane-n-hexane binary system.
It should be pointed out that the onset and evolution of LLV behavior in mixtures is related to the evolution of SLV behavior in those same systems.Thus, in natural gas it is often of interest to predict whether a methane-rich stream with one of the solutes (such as an nparaffin of carbon number four or higher; or benzene, or carbon dioxide) will form a solid phase.The reason behind that is that the presence of a second solvent can considerably change the solubility of the solid solute, if the solute is a hydrocarbon; this also occurs to a lesser extent if the solute is carbon dioxide.
On the other hand, the presence of nitrogen as a second solvent reduces the solubility of solids in methane, ethane, and mixtures thereof.Furthermore, the addition of nitrogen to miscible LNG systems can induce immiscibility.The number of nitrogen binary systems relevant to LNG that exhibit LLV immiscibility is few  nitrogen-ethane and nitrogenpropane, being among the most prominent ones.However, LLV phenomenon has been observed at certain conditions in many ternary and higher realistic nitrogen-rich LNG systems since LLV behavior can and does occur in multicomponent systems even when for none of the constituent binaries themselves an LLV locus is reported.
Another interesting example is the ternary mixture methane + ethane + n-octane; the species in the constituent binary pairs (methane + ethane) and (ethane + n-octane) are too similar in molecular nature to be LLV immiscible, while the pair methane + n-octane is too dissimilar to be immiscible.If a multicomponent mixture is considered to be solute plus solvent in the pseudo-component sense, it can be readily seen why the ternary mixture has a region of LLV behavior.The methane + ethane form a solvent background in which n-octane is immiscible.
The type of the LLV region displayed by a system depends on whether it contains an immiscible binary or not.For a ternary system with no constituent binary LLV behavior present the three-phase region is a "triangular" surface in the thermodynamic phase space with two degrees of freedom, while its boundaries have one.It is bounded from above by a K-point locus; from below by a LCST locus and at low temperatures by a Q-point locus.The systems methane + n-butane + nitrogen (Merrill et al., 1984a) and methane + n-pentane + nitrogen (Merrill et al., 1984b) belong to this class.
A K-point occurs when a liquid phase and a vapor phase become critical in the presence of a heavier noncritical equilibrium liquid phase, whereas an L-point occurs when two liquid phases becomes critical in the presence of a noncritical equilibrium vapor phase.These points, also called upper critical end points (UCEPs) and lower critical end points (LCEPs), are different depending on their location in the pressure-temperature space.That is, an UCEP is always located at a higher temperature and pressure than the LCEP, and when there exists only one critical end point, so that a K-point is always an UCEP point.On the contrary, an Lpoint can be an UCEP or an LCEP, depending on the global phase behavior of the system.
The system methane + ethane + n-octane does not exhibit immiscibility in any of its binary pairs.Although immiscibility has been reported in binary systems of methane + n-hexane and methane + n-heptane, solutes such as n-octane and higher normal paraffins crystallize as temperature decreases before any immiscibility occurs.On the other hand, with ethane as solvent, solutes beginning with n-C 19 and higher paraffins demonstrate LLV behavior.Apparently, the addition of modest amounts of ethane to methane creates a solvent mixture exhibiting immiscibility with n-octane.

Modeling of the phase behavior of LNG systems
The success of the design and operation of separation processes in the oil and gas industry at low temperatures is critically dependent on the accurate descriptions of the thermodynamic properties and phase behavior of the concerned multicomponent hydrocarbon mixtures with inorganic gases.Phase-split calculations and phase stability analysis in natural gas systems simulation can take up as much as 50 % of the CPU time.In complicated problems it may take even more.Thus, it is important to develop a reliable thermodynamic modeling framework (TMF) that will be able to predict, describe and validate robustly and efficiently the complex phase behavior of LNG mixtures.
The TMF has three main elements: a library of thermodynamic parameters pertaining to pure-substances and binary interactions, thermodynamic models for mixture properties, and algorithms for solving the equilibrium relations.Reliable pure-component data for the main constituents of LNG systems are available experimentally; an equation of state is usually the primary choice for the thermodynamic model.Thus, the focal point of a TMF for phase behavior calculations of LNG systems is the availability of robust methods for thermodynamic stability analysis, and of reliable efficient and effective flash routines for three phase split calculations.

Computational technique 1
This first technique uses an efficient computational procedure for solving the isothermal multiphase problem by assuming that the system is initially monophasic.A stability test allows verifying 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 (Michelsen, 1982b;Nghiem and Li, 1984), 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 from the tangent plane at z , is considered (Baker et al., 1982;Michelsen, 1982a).In terms of fugacity coefficients, i  , this criterion for stability can be written as (Michelsen, 1982a)   where i  are mole numbers with corresponding mole fractions as , and 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.Consequently, a minimum of

 
F ξ should be considered in the interior of the permissible region test condition 1 for all trial compositions is not physically possible; it is thus sufficient to test the stability at all stationary points of   F ξ since this function is not negative at all stationary points.Here, the quasi-Newton BFGS minimization method (Fletcher, 1980) was 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 p  phases, the equilibrium calculation is solved by minimization of the following function subject to the inequality constraints given by 1 () 1 1,..., and where i z is the mole fraction of the component i in the system, () x  is the mole fraction of 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.325kPa).In eq 3 the variables () Equation 3 is solved using an unconstrained minimization algorithm by keeping the variables () i n  inside the convex constraint domain given by eqs 4 and 5 during the search for the solution.In this case, a hybrid approach to minimize eq 3 is used 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 the quasi-Newton BFGS 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 (Justo-García et al., 2008a).

www.intechopen.com
Phase Behavior Prediction and Modeling of LNG Systems with EoSs -What is Easy and What is Difficult?365

Computational technique 2
The second approach to calculate multiphase equilibria used here applies a rigorous thermodynamic stability analysis and a simple and effective method for identifying the phase configuration at equilibrium with the minimum Gibbs energy.The rigorous stability analysis is exercised once and on the initial system only.It is based on the well-known tangent plane criterion (Baker et al., 1982;Michelsen, 1982a) but uses a different objective function (Stateva and Tsvetkov, 1994;Wakeham and Stateva, 2004).The key point is to locate all zeros (  y ) of a function    y given as where Furthermore, the number * k , which geometrically is the distance between two such hyperplanes, can be either positive or negative.A positive * k corresponds to a zero, which represents a more stable state of the system, in comparison to the initial one; a negative * k , a more unstable one.When all calculated * k are positive, the initial system is stable; otherwise, it is unstable.
It is widely acknowledged that the task to locate all zeros of the tangent-plane distance function (TPDF), (   y) in this particular case, is extremely challenging because a search over the entire composition space is required.The search is further complicated by the existence of a number of trivial solutions, corresponding with the number of equilibrium phases present (Zhang et al., 2011).The specific form of (  y) (its zeros are its minima) and the fact that it is easily differentiated analytically, allows the application of a non-linear minimisation technique for locating its stationary points, and in their works Stateva and Tsvetkov (1994) and Wakeham and Stateva (2004) used the BFGS method with a line-search (Fletcher, 1980).The implementation of any non-linear minimization technique requires a set of "good" initial estimates, and the BFGS method is no exception.All details of the organization and implementation of the initialization strategy employed by the stability analysis procedure are given elsewhere (Stateva and Tsvetkov 1994) and will not be discussed here.
Thus, as discussed by Wakeham and Stateva (2004), a method has been created which leads, in practice, to an "extensive" search in the multidimensional composition space.It has proved to be extremely reliable in locating almost all zeros of (  y) at a reasonable computational cost.The term "almost all" zeros is used because there is no theoreticallybased guarantee that the scheme will always find them all.If, however, a zero is missed, the method is self-recovering.Furthermore, the TPDF is minimized once only, which is a distinct difference from the approach that stage-wise methods generally adopt.Since stability analysis on its own cannot determine unequivocally which is the stable phase configuration of a system (identified as unstable at the given temperature and pressure), it is suggested to run a sequence of two-phase flash calculations to determine the correct number of the phases at equilibrium, and the distribution of the components among the phases.LLV   .

Calculation of K-and L-points
Needless to say that it is costly and time-consuming to determine the K and LCST loci in ternary systems experimentally; thus, the availability of appropriate algorithms and numerical routines in the third element of the TMF that will allow the prediction and reliable location in the thermodynamic phase space of such points, is indispensable in the study of the complex phase behavior of LNG model systems.Among the several such algorithms published in the open literature we have chosen to outline briefly and implement those of Gregorowics and de Loos (1996) and Mushrif and Phoenix (2006).In our choice we have been guided by the fact that the above algorithms can successfully predict the K-and L-points of binary and ternary systems.We will thus test their robustness and efficiency in the locating critical end points in LNG model systems.

Gregorowicz and de Loos' algorithm
In a study on the modeling of the three-phase LLV region for ternary hydrocarbon mixtures with the SRK EoS, Gregorowicz and de Loos (1996) proposed a procedure for finding K-and L-points of ternary systems, based on the solution of thermodynamics conditions for the Kand L-point using the Newton iteration technique and starting points carefully chosen.They applied their procedure to calculate the K-and L-point loci for two ternary systems, namely, C 2 + C 3 + C 20 and C 1 + C 2 + C 20 , in which the constituent binary C 2 + C 20 exhibit immiscibility.Consequently, the extension of the three-phase LLV region of these systems is bounded by the binary 12 LLV   locus of the system C 2 + C 20 and the ternary K-point and L-point loci.Briefly, the strategy followed by these authors to find the K-and L-points of the two ternary systems was the following: (1) calculation of the critical line, the K-point, the three phase line, and the L-point for the system C 2 + C 20 using thermodynamic conditions, and (2) calculation of the K-and L-point loci for the ternary systems by using as starting points the obtained coordinates of the K-and L-point calculations for the binary system C 2 + C 20 .In this case, to obtain a K-and L-point for a ternary system, the following set of six nonlinear equations, TV V x x x a n dx   have to be solved, where  designates V for the L-point or L 2 for the K-point, D and * D are the two determinants that must be satisfied at a critical point, and i  is the chemical potential of component i .
The critical criteria given by eqs 10 and 11 are based on the Helmholtz energy, which can be expressed as (Baker and Luks, 1980), where i n is the mole number of component i and V is the system volume.Derivatives of the Helmholtz energy are denoted by a subscript in eqs 10 and 11 (e.g., V A , i x A ) indicating the differentiation variable (volume V or mole fraction i x of component i ).Details to obtain the elements of determinants D and * D are given in Baker and Luks (1980).
It is worth mentioning that Gregorowicz and de Loos (1996) calculated the ternary K-and Lpoints at chosen values of the temperature, which it is important when the experiments are carried out isothermally.

Mushrif and Phoenix's algorithm
The second approach to calculate K-and L-points was proposed by Mushrif and Phoenix (2006).This approach utilizes an efficient critical point solver and a standard phase stability test within a nested-loop structure to directly locate K-and L-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.
This algorithm is similar to that proposed by Gauter et al. (1999) to calculate critical end points (CEPs) for ternary systems of carbon dioxide as the near-critical solvent and two lowvolatile solutes (1-pentanol or 1-hexanol + n-tridecane) with the PR (Peng and Robinson, 1976) EoS.These authors used the approach of Heidemann and Khalil (1980) to follow a critical line in steps of the mole fraction of the carbon dioxide, searching for a CEP along this line; i.e., searching for the occurrence of an additional phase in zero amount by using the formulation of Michelsen (1982a).
The principal difference between the algorithm of Mushrif and Phoenix and that of Gauter et al. is that the former directly calculate the K-and L-points without following the critical lines.
In the algorithm proposed by Mushrif and Phoenix (2006), the critical criteria are based on the tangent-plane criteria developed by Michelsen and Heidemann (1988) at specified temperature and pressure.Michelsen and Heidemann (1988) also formulated the criticalpoint criteria in terms of the tangent plane distance based on the Helmholtz energy at fixed temperature and volume.The condition for the stability of a mixture, with respect to a trial phase n , is where F is the tangent plane distance from the Gibbs energy surface to the hyperplane tangent of the Gibbs energy surface at the composition z .
In eq 15 the sign of F will determine the stability of the test phase; i.e., if 0 F  the test phase is stable, if 0 F  the test phase is in equilibrium with some alternate phase, and if 0 F  the test phase is unstable.Michelsen (1984)


u are determined by inverse iteration (Wilkinson, 1965), then min 2 b   . If 0 b  , the system is at the limit of intrinsic stability.At a critical point coefficients b and c in eq 16 are zero for a given eigenvector u of B corresponding to the smallest eigenvalue min  .For the evaluation of coefficient c , Michelsen (1984) showed that this can be determined efficiently from information already available of i u and i g .
The solution procedure to calculate a critical point is as follows: since coefficient b is to be a zero eigenvalue of B , then it must be a singular matrix with a zero determinant; i.e., The criterion of 0 b  is met when the matrix B is singular and is used to find the critical temperature at a fixed composition and volume.The determinant of matrix B is calculated through a LU decomposition of B (  LU B where L is lower triangular and U is upper triangular); i.e., the determinant of B is the product of the diagonal elements of the LU decomposed matrix.Once the iteration to find the stability limit has converged, the vector u is determined by inverse iteration technique.
The implementation of the equation-of-state approach for calculating critical points using this procedure requires that temperature T and volume V are iterated in a nested way.That is, based on an initial guess of V , the temperature is determined in a inner loop until the determinant of the matrix B becomes equal to zero; then the convergence criterion for the coefficient c is checked.If this coefficient, evaluated at the stability limit, is equal to zero, the calculation ends; otherwise, a new estimate for the volume is generated in an outer loop and the iteration on T is evaluated again.Once T and V have been obtained, the pressure P is evaluated from the equation of state.
After having calculated the critical point, a noncritical equilibrium phase is searched for at constant temperature and pressure conditions.Mushrif and Phoenix (2006) used the stability test implemented by Michelsen (1982a)  The tangent plane to the Gibbs energy surface at the trial composition is parallel to the reference-phase tangent plane (critical phase) when crit crit ln ln ln( ) ln 0 where 1 ln( ) is the dimensionless distance between the two tangent planes and i  is the fugacity coefficient evaluated at composition y .
If 0   , the trial phase is in equilibrium with the reference phase; if 0   , the trial phase is an incipient phase, and if 0   , the reference phase is unstable.By combining eq 17 with the definition of  , the set of N equations to solve for a stationary point can be written as crit exp ln ln ( ) ln 0 1,..., which can efficiently be solved by Newton iteration or through a minimization method.
A K-or L-point is found if 0   .When  does not meet the convergence criterion (e.g., where the derivative 1 dd z  is approximated by perturbing the critical composition, recalculating  from eq 17 and calculating the finite difference analogue of the derivative.Mushrif and Phoenix (2006) have pointed out that failure of the algorithm can occur when a critical composition is updated to a value where no stationary point exist other than the trivial solution.
To calculate a K-or L-point using this algorithm, it is necessary to provide appropriate initial estimates of composition, temperature, and volume.In this case, good initial guesses for critical temperatures and volumes were, depending on the type of calculation, the same as those used by Heidemann and Khalil (1980).The success of the algorithm to locate a K-or L-point strongly depend on (1) the binary interaction parameters used in the equation of state and (2) the initial critical composition (0)  z .However, it would seem that the value of the initial critical composition significantly affects the successful convergence of the method to locate a K-or L-point.

Thermodynamic models
Modeling of the complex phase behavior of LNG systems requires a suitable thermodynamic model and a robust and efficient computational algorithm for performing phase stability and multiphase flash calculations interwoven in the second element of the TMF.Regarding the thermodynamic models, the SRK EoS and the PC-SAFT EoS have received wide acceptance in the industry because of their ability to predict accurately the phase behavior of oil-gas systems.

The SRK equation of state
The explicit form of the SRK equation of state (Soave, 1972) For mixtures, constants a and b are given by 11 1 ; and ij a is defined as where ij k is an adjustable interaction parameter characterizing the binary formed by components i and j .
Eq 20 can be written in terms of compressibility factor, ZP v R T  , as where Aa PR T  and

 
Bb PR T  .
The expression for the fugacity coefficient, ii i f yP

 
, is given by

The PC-SAFT equation of state
In the PC-SAFT EoS (Gross and Sadowski, 2001), the molecules are conceived to be chains composed of spherical segments, in which the pair potential for the segment of a chain is given by a modified square-well potential (Chen and Kreglewski, 1977).Non-associating molecules are characterized by three pure component parameters: the temperatureindependent segment diameter  , the depth of the potential  , and the number of segments per chain m .
The PC-SAFT EoS written in terms of the Helmholtz energy for an N-component mixture of non-associating chains consists of a hard-chain reference contribution and a perturbation contribution to account for the attractive interactions.In terms of reduced quantities, this equation can be expressed as The hard-chain reference contribution is given by 1 (1 ) l n ( ) where is the mean segment number in the mixture The Helmholtz energy of the hard-sphere fluid is given on a per-segment basis as where k is the Boltzmann constant and T the absolute temperature.
The dispersion contribution to the Helmholtz energy is given by where hc Z is the compressibility factor of the hard-chain reference contribution, and 23 3 11 The parameters for a pair of unlike segments are obtained by using conventional combining rules where ij k is a binary interaction parameter, which is introduced to correct the segmentsegment interactions of unlike chains.
where the coefficients i a and i b depend on the chain length as given in Gross and Sadowski (2001).
The density to a given system pressure sys P is determined iteratively by adjusting the reduced density  until sys calc PP  .For a converged value of  , the number density of molecules  , given in Å -3 , is calculated from Using Avogadro's number and appropriate conversion factors,  produces the molar density in different units such as The pressure can be calculated in units of 2 Pa N m   by applying the relation from which the compressibility factor Z , can be derived.The expression for the fugacity coefficient is given by   In eq 40, the partial derivatives with respect to mole fractions are calculated regardless of the summation relation 1 1

Results and discussion
Experimental data reported by Llave et al. (1987) for the system nitrogen + methane + ethane, by Hottovy et al. (1981) for the system methane + ethane + n-octane, and by Fall and Luks (1988) for the system carbon dioxide + nitrogen + n-nonadecane, were used to test and compare the robustness, efficiency and reliability of the two computational techniques and the SRK and PC-SAFT EoS thermodynamic models embedded in the respective elements of the TMF.The prediction and modeling of the phase behavior of these systems demonstrates in a clear-cut way the usual numerical difficulties encountered in the process.
The binary interaction parameters used with the PC-SAFT equation were taken from García-Sánchez et al. ( 2004) and from Justo-García et al. (2008b), while those used with the SRK EoS were taken from Knapp et al. (1982).Some of the interaction parameters were also obtained from the minimization of the sum of squares of the differences between experimental and calculated bubble-point pressures.
The binary interaction parameters employed are: k   for the PC-SAFT equation, respectively.The components' physical properties required for the calculations performed with the SRK EoS were taken from DIPPR (Rowley et al., 2006) 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 (2001).

The nitrogen + methane + ethane system
The three-phase VLL region displayed by this ternary system is bounded from above by a K-point locus, from below by a lower critical solution temperature LCST locus, at low temperatures by a Q-point locus, and, due to the fact that this system contains a binary pair (nitrogen + ethane) which exhibits LLV behavior, its LLV space is truncated.In this case, the partially miscible pair nitrogen + ethane spans the LLV 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 LLV space which extends from the binary LLV locus upward in temperature.The topographical nature of the regions of immiscibility for the system nitrogen + methane + ethane is shown in Fig. 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).Fig. 2 presents the experimental and calculated L 1 -L 2 -V phase behavior (in terms of L 1 -L 2 nitrogen mole fraction data) for the nitrogen + methane + ethane system at 135 K and different pressures.This figure shows a reasonable agreement between the experimental values of liquid phases L 1 and L 2 and those predicted with both models.Notwithstanding, although the LLV calculations performed with both equations up to a position near the Kpoint (about 41.25 bar with both models), this point is away from the experimental one (43.05bar at 135 K).
An attempt to directly calculate either the K-or L-point for this ternary system using the algorithm of Mushrif and Phoenix (2006) was carried out.However, the algorithm was not able to give correct values of these critical end points.This is because the algorithm is strongly initialization dependent and hence gives different values of these points, depending on the initial guess of the critical composition, which meets the convergence criterion.Our preliminary results show that the algorithm advocated by Gregorowicz and de Loos (1996) is more stable than that of Mushrif and Phoenix, even if it also depends on the initial values of the critical composition.LCST (L1=L2-V) Fig. 1.P-T space of the boundaries of the three-phase L 1 -L 2 -V regions for the system nitrogen + methane + ethane.Experimental data from Llave et al. (1985Llave et al. ( , 1987)).models agree with each other for the liquid L 2 phase but represent the experimental data very closely.Since there are not experimental data below 22.16 bar, 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 LCST point with the PC-SAFT model (37.20 bar) seems to be closer to the "hypothetical" experimental LCST point (38.37 bar) in comparison with the LCST point estimated from the SRK model (35.45 bar).Of course, these discrepancies can be due to the fact of using binary interaction parameters determined from VL equilibrium data, which, apparently, led to less accurate results.

The methane + ethane + n-octane system
The three-phase LLV region displayed by the methane + ethane + n-octane ternary system (a surface in the thermodynamic phase space with two degree of freedom) is bounded from above by a K-point locus (L-L=V), from below by a lower critical solution temperature LCST locus (L=L-V), and at low temperatures by a Q-point locus (S-L-L-V).For the three components in this system, there is no binary immiscibility.The topographical nature of the regions of immiscibility for this system is shown in Fig. 3, where symbols are the experimental data given by Hottovy et al. (1981) identifying the boundaries of the three phase LLV region for this ternary system.This Figure shows also that the L-L=V and L=L-V critical end-point loci intersect at a tricritical point (L=L=V) at the upper temperature limit.
Fig. 3. P-T space of the boundaries of the three-phase L 1 -L 2 -V regions for the system methane + ethane + n-octane.Experimental data from Hottovy et al. (1981).
Following the immiscibility region, a single temperature was chosen to test the capabilities of the PC-SAFT with computational technique 1, and the SRK EoS with computational technique 2, to predict the phase behavior for the system methane + ethane + n-octane.Fig. 4 compares the performance of the two methods at 210 K, and at different pressures on the basis of the experimentally measured and calculated nitrogen mole fractions for the liquid L 1 and L 2 phases by the two thermodynamic models employed.In this case, the predictions of the PC-SAFT model are closer to the experimental composition values than those of the SRK model.However, it should be mentioned that it was not possible to continue the calculations with this model to approach either the K-or L-point because the three-phase LLV triangles become so very narrow as pressure decreases or increases that it is extremely difficult to determine an appropriate global composition able to separate this mixture into three-phase LLV equilibria.On the other hand, because the three-phase LLV triangles predicted with the SRK model are wider than those predicted with the PC-SAFT model, it was easier to get a good initial global composition to calculate the three-phase LLV equilibria from the LCST point (52.80 bar) to the K-point (59.09 bar) applying technique 2.
An inspection of this figure shows that the predictions of both EoSs don't follow the behavior of the liquid phase L 1 as well as the variation of liquid phase L 2 as pressure decreases.Also, it is interesting to note that although there is not a true experimental value of the LCST at the temperature considered, Fig. 4   Fig. 4 also shows that at the different pressures the predicted L 1 -L 2 -V region (in terms of L 1 and L 2 methane mole fraction data) with the SRK model deviates considerably from the experimental one, while the PC-SAFT model predictions are more reasonable.However, since the interaction parameters for the SRK and PC-SAFT models were determined from binary vapor-liquid equilibrium data, the rather poor fit in this region with either model is not unexpected.

The carbon dioxide + nitrogen + n-nonadecane system
As mentioned in Section 2, the type of the LLV region displayed by a ternary system depends on whether it contains an immiscible pair or not.In this context, the system carbon dioxide + nitrogen + n-nonadecane exhibits immiscibility in the carbon dioxide + nnonadecane binary pair (Fall et al., 1985), so that its three-phase region is similar to that exhibited by the system nitrogen + methane + n-pentane (Merrill et al., 1984b).Therefore, the LLV region is "triangular" and is bounded from above by a K-point locus (L-L=V), at low temperatures by a Q-point locus (S-L-L-V), and, from a position of the Q-point locus to a position on the K-point space, by a binary carbon dioxide + n-nonadecane LLV locus.
Fig. 5 presents the experimental pressure-temperature diagram of the LLV space displayed by the system (Fall and Luks, 1986).An examination of the figure shows that this system does not have a LCST (L=L-V) locus and that the Q-point locus terminates at an invariant V-L1-L2 (Binary CO2-C19) Fig. 5. P-T space of the boundaries of the three-phase L 1 -L 2 -V regions for the system carbon dioxide + nitrogen + n-nonadecane.Experimental data from Fall et al. (1985) and Fall and Luks (1986).379 point of the S-L-L=V type.Thus, due to the fact that the carbon dioxide + nitrogen + nnonadecane system contains a binary pair (carbon dioxide + n-nonadecane) exhibiting LLV behavior, its "triangular" LLV region is a three-sided space without a tricritical point.
A temperature of 297 K was chosen to study this ternary system and the results obtained are presented in Fig. 6.This figure shows that the PC-SAFT model predicts well the experimental carbon dioxide compositions of the liquid L 2 and vapor phases down to the lowest measured pressure of 62.63 bar (i.e., the binary carbon dioxide + n-nonadecane data) for this isotherm.However, the SRK model is superior to the PC-SAFT model in predicting the three phases in equilibrium for this ternary system.In this case, the calculated L 1 , L 2 , and V phases are close to the experimental ones.Furthermore, though the SRK EoS overpredicts the "experimental" K-point (87.30 bar) by 4.39 bar the performance of the PC-SAFT EoS in this particular case is inferior as it overpredicts by 14.6 bar the "experimental" point.Fig. 6.Comparison of L 1 and L 2 compositional data of carbon dioxide at 297 K for the system carbon dioxide + nitrogen + n-nonadecane.Experimental data from Fall and Luks (1986).
Nevertheless, it should be recalled that all calculations were performed using binary interaction parameters obtained from regression of binary experimental VL data, many of them measured at temperatures higher than that studied here.Furthermore, we are confident that the performance of the corresponding thermodynamic models could be improved considerably provided the interaction parameters were obtained from the regression of the experimental data at three phases.Still, if those sets of interaction parameters are used to predict the phase behavior of a given system at conditions different from the original ones then there is the risk that the equilibria predictions and calculations could either give physical meaningless results or fail altogether.
We also tried to directly calculate the K-points at the temperature considered for this ternary system by using the algorithm of Mushrif and Phoenix (2006); however, once again, the algorithm predicted different values of the critical end points, depending on the initial critical composition.In this case, the strategy to find a critical temperature close to the temperature of study was to use a series of initial compositions.Unfortunately, none of the compositions produced a K-point similar to those obtained from the VLL calculations with both models.
Finally, it should be pointed out that for the sake of comparison both equations of state were interwoven into computational procedure 1 and 2, respectively, and that a series of multiphase flash calculations were carried out at different temperatures and pressures for the three ternary systems studied obtaining the same results for the specific EoS, irrespectively of the computational procedure utilized.
The results obtained showed that there are not any essential differences between, or particular advantages of any of two computational procedures, either in their efficiency, effectiveness, robustness or in their convergence behavior.Thus, it can be said that both procedures 1 and 2 can be used to predict the phase behavior of a wide variety of multicomponent nonideal systems over wide ranges of temperature and pressure.

Conclusions
Though there has been much progress and advance in two-phase stability and two-phase split calculations with EoSs, there is still not much progress in three-phase split calculations in natural gas systems despite the large number of publications devoted to the subject.The reason behind that is that the difficulties and challenges are dominating over the easy to perform calculations, if any.Thus, it can be accepted that the algorithms and numerical methods advocated are not robust enough for incorporation in a process simulator.In view of this, the further development of a reliable, robust and efficient TMF and its subsequent approbation on model systems, typical representatives of LNG, is of considerable interest both to scientists and engineers.
On the example of three ternary systems (nitrogen-methane-ethane, methane-ethane-noctane, and carbon dioxide-nitrogen-n-nonadecane) the capabilities of a TMF, advocated by us, to predict and model complex phase behavior of systems of importance to LNG processing are demonstrated.The TMF employs two numerical techniques which embed different thermodynamic models and stability analysis routines.
The results obtained show that there are many and different challenges and difficulties that are not always possible to overcome completely.For example, the two techniques for multiphase flash calculations cannot always assure steady and non-oscillatory convergence with no tendency towards a strong attraction to the trivial solution, particularly in cases close to the critical lines.Besides, it is known that the phase equilibrium equations are often difficult to converge in the critical region and that the use of inappropriate initial estimates can lead to the trivial solution.In view of this, the availability in a TMF of a robust stability analysis routine that will provide good set of initial estimates for the compositions of possible equilibrium phases and will guarantee steady convergence of the flash routines is of great importance.Regarding the prediction of K-and L-points with the Gregorowicz and de Loos and Mushrif and Phoenix's methods, it is clear that the both algorithms are strongly initialization dependent.To overcome this problem, for example, Mushrif and Phoenix suggest that the critical phase composition is updated based on the values of parameter  , calculated from equilibrium phase calculations using the Newton-Raphson iteration.However, during this process of upgrading, the composition may change to a value where there is no phase in equilibrium, particularly when  differs significantly from iteration to iteration.
The procedure of Gregorowicz and de Loos to calculate K-and L-points seems to be a better fitted method to carry out this task.However, the evaluation of the determinants for solving the conditions of criticality requires the second derivatives of the Helmholtz energy with respect to volume and composition, which makes difficult their evaluation, particularly when these derivatives have to be found analytically applying the PC-SAFT EoS.Still, of course, a possible solution to this problem is to evaluate these derivatives numerically.
Finally, both the PC-SAFT EoS and the RKS cubic EoS are capable of representing with a reasonable accuracy the experimentally observed phase behavior of the ternary systems studied.
using the critical composition z as the reference phase

Fig. 2
Fig.2also shows that at pressures away from the LCST point, the PC-SAFT model gives a better representation of the experimental compositions for liquid phase L 1 while both

Fig. 2 .
Fig.2.Comparison of L 1 and L 2 compositional data of nitrogen at 135 K for the system nitrogen + methane + ethane.Experimental data fromLlave et al. (1987).

Fig. 4 .
Fig. 4. Comparison of L 1 and L 2 compositional data of nitrogen at 210 K for the system methane + ethane + n-octane.Experimental data fromHottovy et al. (1981).
www.intechopen.comPhase Behavior Prediction and Modeling of LNG Systems with EoSs -What is Easy and What is Difficult?381 developed the criteria for critical points by expanding the tangent plane distance function in a Taylor series around the test point as that defines the distance in composition space from the test point at 0 s  .As the sign of the tangent plane distance function F determines the stability of the test phase, it is necessary to find the minimum of this function using scaled mole numbers as s dF ds   hold at the test point; s being a parameter min , can be written as rr TT indicates that the estimated LCST point with the PC-SAFT model (56.35 bar) is closer to the "experimental" one (56.64bar) than that obtained with the SRK model (52.80 bar).Nonetheless, the "experimental" K-point (60.19 bar) is closer to the one calculated with the SRK model.
Phase Behavior Prediction and Modeling of LNG Systems with EoSs -What is Easy and What is Difficult? www.intechopen.com