Perturbation Theory and Phase Behavior Calculations Using Equation of State Models

Equations of State (EoS) live at the heart of all thermodynamic calculations in chemical engineering applications as they allow for the determination of all related fluid properties such as vapor pressure, density, enthalpy, specific heat, and speed of sound, in an accurate and consistent way. Both macroscopic EoS models such as the classic cubic EoS models as well as models based on statistical mechanics and developed by means of perturbation theory are available. Under suitable pressure and temperature conditions, fluids of known composition may split in more than one phases, usually vapor and liquid while solids may also be present, each one exhibiting its own composition. Therefore, computational methods are utilized to calculate the number and the composition of the equilibrium phases at which a feed composition will potentially split so as to estimate their thermodynamic properties by means of the EoS. This chapter focuses on two of the most pronounced EoS models, the cubic ones and those based on statistical mechanics incorporating perturbation analysis. Subsequently, it describes the existing algorithms to solve phase behavior problems that rely on the classic rigorous thermodynamics context as well as modern trends that aim at accelerating computations.


Introduction
Equations of State (EoS) have been widely used in the chemical engineering industry for the calculation of process fluids phase properties. EoS models are algebraic expressions of the form fp , T, v m ðÞ ¼ 0 which relate molar volume v m to pressure and temperature. Since the derivation of the ideal gas law and following the pioneering work of Van der Waals, dozens of EoS models of various complexity and thermodynamic considerations have been presented to accurately estimate thermophysical properties. Among them, basic and extended cubic equations of state, virial forms, EoS models with association terms and models based on statistical physics. Of them, the ones most widely used in the chemical engineering industry are the cubic ones [1] due to their simplicity and speed of calculations, thus minimizing the computing time required for flow simulations in processes, porous media and pipelines. Less simple but more accurate models incorporating associating theory are often used in midstream and downstream applications [2].
EoS models based on the application of statistical mechanics in conjunction to perturbation theory to describe the thermodynamic behavior of substances at a microscopic level are commonly used to estimate properties of liquids [3]. This approach is based on studying the microscopic behavior of a set of molecules by considering ensembles comprising of many instances of the set. Subsequently, the system energy and eventually all thermodynamic properties of interest are obtained by treating statistically the ensemble properties. However, as the derivation of a closed form of the energy function is usually intractable, perturbation theory greatly simplifies that task. A known closed form solution for a simple reference system is firstly adopted, and the additional energy terms required to improve the simple reference system to the complex one are considered as a perturbation of the original reference system. Perturbation theory utilizes linearization to lead to approximate closed form solutions of the combined complex system.
To obtain estimates of the thermophysical properties using EoS models, it is necessary that the composition of the fluid is known and that a reliable characterization of the mixture components, by means of specific components properties values, is available. Cubic EoS models though simple they are not predictive, and the reliability of their predictions can only be ensured by "tuning" the model, i.e. varying the components properties so that the model predictions match accurately the available experimental measurements.
Once a tuned EoS model is available, properties such as density, fugacity coefficient, enthalpy, heat capacity, Joule-Thomson coefficient and speed of sound can be easily computed by simple expressions. Calculations become more complex when the phase state of a mixture is not known a priori. As an example consider a control volume, i.e. a grid block, in a flow simulation model where the pressure and saturation change of each coexisting phase at current timestep need to be determined in order to get a description of the fluid state. The pressure change in the control volume is related to mass influx and outflux through fluids density and compressibility. When the control volume content is a single-phase fluid both properties can be easily computed by means of the EoS model. However, when the content is saturated, it will split into two or more phases, each one exhibiting its own properties, thus introducing the need to identify the number and composition of the equilibrating phases, hence their density and compressibility.
In such cases, a test to determine if the fluid appears in a single or two phases needs to be run, known as stability test [4]. If the test indicates the presence of two or more phases in equilibrium, the phase split problem further needs to be solved to compute the composition and the amount of the two coexisting properties [5]. By knowing their composition, all properties of the equilibrium phases can then be computed regularly.
In this chapter, the utilization of EoS models of the cubic form and those based on perturbation theory is discussed and their application to compute fluids thermophysical properties is presented. Algorithms to run phase stability and phase split in the classic context as well as in the reduced variables one are also discussed. Additionally, the chapter discusses the recent developments in the use of soft computing techniques to accelerate the solution of the stability and phase split problems in flow simulations.

The PR and SRK cubic EoS models 2.1 Development of the cubic EoS models
The ideal gas law pv m ¼ RT, where the gas constant R ¼ k B N A is defined as the product the Boltzmann constant and the Avogadro number, only considers the elastic collision of molecules thus considering the thermodynamic behavior of the fluid as a purely kinetic process. As a result, it exhibits accurate predictions of the molar volume only when gases at pressures and temperatures close to the atmospheric ones are considered. On the other hand, the real gas law pv m ¼ ZRT can be used to describe accurately the properties of any fluid and at any conditions provided that the appropriate value of the compressibility factor Z (also known as deviation factor in the sense that it considers the deviation of the real gas law from the ideal gas one) can be computed. Clearly, the real gas law simplifies to the ideal one by simply setting Z ¼ 1.
Van der Waals was first to recognize the need to separately consider attractive and repulsive forces between the fluid molecules thus leading to the first cubic equation Indeed, the a term in Eq. (1) can be thought of as a term accounting for the attractive forces between molecules as it reduces pressure. Parameter b accounts for the molecules volume which becomes significant at high pressures (i.e. liquid state) as lim p!∞ v m ¼ b. Both parameters are functions of the properties of the component or mixture under consideration. Clearly, by setting both parameters to zero we revert back to the ideal gas law.
Ever since, various new cubic EoS models have been proposed with the Soave-Redlich-Kwong (SRK) and the Peng-Robinson (PR) ones [6] being by far the most commonly used ones in the chemical engineering industry. Both are pressure explicit and are defined by the following expression where the parameters values are given in Table 1. The temperature dependent term in that Table is given by where m is a function of the component acentric factor ω defined by m ¼ 0:48 þ 1:574ω À 0:176ω 2 SRK 0:37464 þ 1:54226ω À 0:26992ω 2 PR, ω ≤ 0:49 0:3796 þ 1:485ω À 0:1644ω 2 þ 0:01667ω 3 PR, ω > 0:49 The required properties of pure components can be found in any standard petroleum thermodynamics textbook [7]. When pseudo-components are used to describe the fluid composition, such as such as pseudo-C 8 and pseudo-C 11 in petroleum mixtures, average values can also be obtained from the literature. Custom pseudo-components such as petroleum mixtures heavy end need to be treated by means of suitable correlations which utilize molar mass and density to provide estimates of the critical properties and the accentric factor or other required properties [8]. When it comes to mixtures, parameters mixing rules need to be utilized to estimate a and b. For a mixture of known composition z i , they are given by The Binary Interaction Parameters (BIP) k ij account for the interaction between different constituents and are usually initialized either to zero or by the Prausnitz [9] rule where the critical molar volume is obtained by solving the EoS at critical conditions and the critical value Z c of the compressibility factor for the PR EoS equals to 0.3074. Parameter θ is user dependent and is usually set to 1.2. Note Eq. (6) is only used to determine BIPs between hydrocarbon components. BIPs between nonhydrocarbons or between hydrocarbon and nonhydrocarbon components are taken from Tables [6].
Once all parameters have been estimated for a mixture of known composition at fixed pressure and temperature, the EoS can be solved for volume. Usually a dimensionless form that can be solved for Z ¼ pv m =RT rather than for v m is preferred where the dimensionless EoS constants are given by

Use of the cubic EoS models
As soon as the EoS constants have been defined, the compressibility factor Z can be obtained by solving the cubic polynomial Eq. (8) [10]. When more than one real positive roots are obtained, the smallest one is selected when the fluid is a liquid whereas the largest one is used for a gas. Molar volume and density can be easily computed by where M denotes the fluid molar mass. Components fugacity coefficients φ i , hence fugacity f i ¼ φ i z i p, can be computed by the following expressions Derivative properties such as the Joule-Thomson coefficient μ JT can be computed by differentiating the EoS and incorporating the derivatives in the rigorous thermodynamic definitions of the properties. For example,

Volume translation
Cubic EoS models are notoriously known for their deficiency in estimating liquid density. A simple modification, known as volume shifting or volume translation, originally proposed by Peneloux [11], can greatly improve the capabilities of cubic EoS. The idea lies in "shifting" the predicted phase molar volumes v EoS m by some amount that depends on the fluid composition and its components properties. More specifically, the shifted volume is given by Parameters c i are component specific and they are usually given as functions of the covolume parameters b i , that is where values of s i for common pure components are available in Tables [6]. It should be noted that "shifting" (or "translating") the volume also affects the Z factor which needs to be updated to ensure calculations consistency It can be shown that when applying volume translation to two phases that equilibrate, the fugacities of the components do change but they do in the same amount so that they remain equal, thus not disturbing the equilibrium. As a result, volume translation does not affect phase compositions in flash calculations or saturation conditions but only phase density.

EoS models in the thermodynamic perturbation theory context
Unlike macroscopic EoS models such as those described in the previous section, major efforts have been oriented toward the development of microscopic approaches based on statistical mechanics where the individual behavior of each particle in a fluid substance is considered. The repulsive and attractive forces are handled separately and combined to provide a description of the thermodynamic properties of fluids through methods based on statistical physics. The basic idea is to study the microscopic behavior of a set of molecules by considering many instances of the set, each one corresponding to one possible state. This ensemble is described through the statistical properties averaged over all possible states. The basic components for this task is the pair potential function ur ðÞand the pair correlation function gr ðÞrespectively, both functions of the distance r away from the center of some molecule. By defining them one can generate expressions to compute the system free energy and eventually all thermodynamic properties of interest [3].
Arriving to the energy expression while starting from ur ðÞand gr ðÞis a very complex task from the mathematical treatment point of view. Complex expressions of the two functions might correspond to more accurate description of the molecules dynamics but they also lead to intractable mathematical expressions. For this task perturbation theory has greatly enhanced the derivation of EoS models by firstly utilizing known closed form solutions for simple reference functions. Subsequently, the small changes between the accurate ur ðÞand gr ðÞfunctions and the reference ones are treated in a very elegant way by means of perturbation theory thus leading to approximate closed form solutions for complex pair functions [12].

The correlation function formalism to derive EoS models
Consider a thermodynamically large system comprising of a fixed number of molecules, at fixed temperature and volume, which is allowed to exchange heat with the environment. Subsequently, consider a collection of many such probable systems forming what is known as a canonical ensemble. The aggregate thermodynamic properties of such systems can be described as functions of the statistic properties of the ensemble. For this task the canonical partition function is defined by where E i corresponds to the energy of each possible microstate, β ¼ 1=k B T is the thermodynamic beta and k B is the Boltzmann constant. Note that Q is dimensionless and as it will be shown later it relates macroscopic thermodynamic properties of the system to the energy of the microscopic systems forming the ensemble.
For a system comprising of N identical molecules the partition function Q N V, T ðÞ at given volume and temperature is given by [3]. where and Z N V, T ðÞ is known as the configuration integral. It is easy to show that if the system potential energy U N is assumed to be zero then the configuration integral Z N simplifies to the system volume and the application of the related partition function leads simply to the ideal gas law. On the other hand, when U N 6 ¼ 0, it is often represented by a sum of pair-wise potentials, i.e.
Various pair potential models ur ðÞhave been presented with the hard-sphere, the square well and the Lennard-Jones being the most pronounced ones [12]. The hard sphere model assumes that the particles are perfect spheres of diameter σ, the potential at distances less than the sphere diameter is equal to infinite (hence the "hard" sphere) and zero beyond that. Therefore, ur ðÞand the corresponding Boltzmann factor are given by The square-well model [13] further allows for a negative value at some distance beyond the hard sphere diameter: The Lennard-Jones model [14] offers the advantage of being defined by a continuous function of the distance r: ur ðÞ¼4ε l=r ðÞ 12 À l=r ðÞ 6 hi : In the equations above σ is the sphere diameter, parameter γ is used to scale the well width, l is the length parameter and ε is the energy parameter. A detailed description on how to use the hard-sphere model pair potential function to develop an EoS model is given in Section 3.3.

Derivation of fluid properties for specific pair functions
By selecting the pair potential model ur ðÞand incorporating it the configuration integral Z N and eventually to the canonical partition function expression, internal energy can be obtained by noting that By utilizing the pair-wise potential energy model of Eq. (19) it can be shown that ur ðÞ gr ðÞ r 2 dr: Therefore, internal energy can be obtained as a function of the particle properties ur ðÞand gr ðÞ . Clearly, the first term corresponds to the kinetic energy of the particles, that is the ideal gas contribution of the system.
Using similar arguments, pressure can be obtained by as the volume derivative of the configuration integral Z N , that is Again, the first term corresponds to the ideal gas pressure term. The system Helmholtz energy is defined by The chemical potential corresponds to the energy required to add one more particle in the collection it is given by and eventually Given the expression above, entropy can be obtained by

The hard-sphere model
The generic fluid properties expressions derived in the previous section are now applied to the hard sphere model for the pair potential energy. By noting that the derivative of the Boltzmann factor of the hard-sphere model is simply the Dirac delta function [15] and replacing it to the generic properties' expressions of the previous paragraph, it follows for pressure that where the pressure is now split into the ideal gas and the excess part and the packing function η which corresponds to the ratio of the particles volume over the total one is given by The expressions for the other properties of interest are obtained similarly and they are given by

Thermodynamic perturbation theory
Although the hard-sphere pair potential model allows for an explicit calculation of thermodynamic properties of interest, its results are not that accurate mostly due to the inherent simplicity of the hard-sphere model itself. Nevertheless, many researchers have pointed out that the comparison between the experimental structure factor and the one obtained computationally from the hard-sphere model indicates that the two curves are quite close to each other. To get a better match a more complex pair potential model could be sought which, however, would inevitably lead to mathematically intractable expressions of the properties. Alternatively, perturbation methods can be applied to the original simple hard-sphere model to add thermodynamic complexity under controlled extra computational burden.
The idea, firstly presented by Zwanzig [16], is to divide the total potential energy into two terms, U 0 and U p respectively, where the first term corresponds to a reference system and the second one corresponds to the perturbation, which needs to be significantly smaller than the reference one for the perturbation method to be applied successfully. The total energy is then given by The perturbation parameter λ allows for various mixtures of U 0 and U p whereas the original fluid energy is obtained for λ ¼ 1. By replacing that expression to the configuration integral we obtain where the : hioperator denotes the statistical average of the reference system. Replacing Eq. (34) to the expression for the Helmholtz energy we obtain where the first term corresponds to a multiple of the Helmholtz energy of the reference system and the second term accounts for the energy of the perturbation. By combining the Taylor expansion forms of the exponential term and of the logarithmic term we obtain where The treatment above has allowed the energy to be described by simpler expressions of the perturbation energy term based on the c 1 , c 2 , c 3 parameters. Therefore, to get the full energy expression one needs to choose the U p model, compute the values of the c 1 , c 2 , c 3 parameters and replace then in Eq. (36) while setting λ ¼ 1. All thermodynamic properties of interest can then be computed as functions of the energy as shown in Section 3.2.
The beauty of the perturbation theory is that although the calculation of the c 1 , c 2 , c 3 parameters, which have been introduced by the application of perturbation theory and the assumption of a simple reference system, is not an easy task it still is significantly easier than replacing a complex pair potential and pair correlation function and running mathematical operations in Eq. (18).
As an example application of perturbation theory in statistical mechanics based thermodynamics, consider the use of the model obtained by perturbation theory to generate the Van der Waals EoS, this time from a statistical mechanics point of view instead from the classic macroscopic one. Firstly, let us state the following assumptions: 1. The potential energy consists of the sum of all pair potentials: 2. The pair potentials are equal between any pair of molecules: By introducing those assumptions to Eq. (37) the calculation of coefficient c 1 simplifies to To proceed we further need to introduce the following assumptions: 1. The reference system to describe U 0 is the hard-sphere one 2. The particles are uniformly distributed which implies that the pair correlation function g 0 r ðÞis equal to one at any distance beyond the limits of the particle and equal to zero inside that 3. The free fluid volume is V À Nb where b is the particle volume that equals to b ¼ 2=3πσ 3 and we end up with Finally, by utilizing first order approximation only (up to c 1 ), replacing c 1 in the free energy Eq. (36) and differentiating over volume to obtain pressure (i.e. p ¼À∂H=∂Vj T ) the well known Van der Waals equation is obtained Interestingly, the perturbed energy term u p has not been defined explicitly but it has been incorporated into the EoS a parameter. From a perturbation theory point of view, the accuracy of the Van der Waals equation of state can be improved by further considering the c 2 term, the calculation of which, however, is quite more complicated.

The stability test
The question answered by a stability test is whether a mixture of given composition, at given pressure and temperature, will appear as a single phase or as a multi-phase one. Clearly, the question can be answered if the bubble point/upper dew point pressure and the lower dew point pressure (if any) of the mixture at operating temperature is known. Any fluid above its bubble point pressure will appear as single-phase liquid whereas when above its upper dew point or below its lower dew point pressure it will appear as single-phase gas.
As the saturation pressure calculation is very costly, a brilliant approach by Michelsen [4] is most preferably used. The idea lies in the fact that if a mixture is unstable, i.e. if it splits into two or more phases when in equilibrium, there exists at least one composition which when forms a second phase in an infinitesimal quantity leads to a reduction of the system's Gibbs energy. Therefore, one should try a bubble/ drop of any possible composition, consider that as a second phase that coexists with the original fluid and examine whether the system Gibbs energy is reduced compared to that of the original single phase fluid. To avoid looking over all possible compositions, Michelsen suggested that one should only look for compositions that minimize the mixture's Gibbs energy rather than simply reduce it. If all minima lie above the single-phase fluid Gibbs energy, then there is no composition that allows for an energy reduction, hence the fluid is single phase, and otherwise it lies in two-phase equilibrium. The Gibbs energy difference, the sign of the minimum of which is used to determine the fluid phase state, is referred to as the Tangent Plane Distance (TPD).
Locating the minima of the TPD is not an easy task as any optimization algorithm may be trapped in a local positive rather than the global negative minimum, thus leading to wrong conclusion about the number of phases present. Additionally, the stability problem has a natural "trivial solution", the one corresponding to a second phase composition same to that of the original fluid. This solution leads to a zero TPD value and it may attract any optimization algorithm, thus misleading the stability algorithm away from the true TPD minimum.
To overcome those issues two approaches can be envisaged. Firstly, one might use global minimization algorithms which ensure that the minimum found is the global one [17]. Such algorithms take significant time to run hence they can only be applied to single calculations rather than batch ones, as is the case in fluid flow simulation. The second approach considers the repeated run of simple optimization algorithms, each time with appropriate initial values so that the global minimum will be located by at least one of those tries.
Based on the above observations, Michelsen [4] presented an algorithm which constitutes the standard approach to treat phase stability. To simplify calculations, it is recommended to optimize TPD by varying the equilibrium coefficients k i ¼ y i =z i , also known as distribution coefficients, rather than the bubble composition y i itself. The algorithm is as follows The algorithm needs to be repeated, this time by assuming that the feed is a gas and the second phase is a drop. In that case, the algorithm is as follows 1. Assume feed is a liquid and look for a drop, i.e. compute As soon as both calculations have been completed, Table 2 can be used to reckon on the phase state.
The algorithm described above is known as the "two-sided" stability test as the trial phase is tested both from the bubble as well as from the drop side. The bubble test converges to nontrivial negative solutions only when the test pressure and temperature conditions lie within the phase envelope in the range where the feed is predominantly liquid. Similarly, the drop test converges to nontrivial negative solutions only when the feed is predominantly gas. The two ranges overlap in a region known as "the spinodal" [19] where both tests converge to solutions with a negative TPD value indicating instability of the feed (Figure 1).

Vapor phase test Liquid phase test Result
Trivial solution Trivial solution Stable Trivial solution S L > 1 Unstable

The phase split
Once the stability test has indicated that the feed is split into two or more phases, a phase split algorithm, also known as flash, needs to be run to determine the composition and relative amount of each phase present in equilibrium. For the simple case of vapor-liquid equilibrium (VLE) which is most commonly encountered in petroleum engineering applications, the phase split algorithm will provide the compositions of the gas and liquid phase, y i and x i respectively, as well as the vapor phase molar fraction β. At equilibrium, the two phases should satisfy two conditions, namely the mass balance and the minimization of the system Gibbs energy. The first condition simply requires that the mass of each component in the feed should equal to sum of their mass in the resulting two phases in equilibrium, i.e.
The second condition additionally requires the phase compositions to be so that the two-phase system's Gibbs energy, defined by is at its minimum. It is easy to show that setting the Gibbs energy gradient equal to zero, an equivalent condition is obtained which requires that the fugacity of each component in the vapor phase is equal to its fugacity in the liquid phase, i.e.
Finally, we need to ensure that the composition of each equilibrium phase is consistent by summing up to unity. Equivalently, we may require that Summarizing, the solution of the flash problem can be seen as the solution of a system of 2n þ 1 equations, that is Eq. (41), (43) and (44), in 2n þ 1unknowns,i.e.y i , x i and β. Note that the mass balance equations are linear in the phase compositions, hence they can be solved and replaced in Eq. (44) thus allowing the flash problem to be reformulated in terms of the k-values and the molar fraction β. Indeed, by incorporating Eq. (41) and the k-values definition in Eq. (43) to Eq. (44), the famous Rachford-Rice equation is obtained where β i ¼ 1= 1 À k i ðÞ , which can be solved for the molar fraction β. Given β and the k-values, the equilibrium phase compositions can then be obtained by Therefore, the phase split problem can be treated as the solution of a system of n þ 1 equations, that is Eq. (43) and (45) From Eq. (45) it can be seen that the Rachford-Rice equation is a monotonically decreasing one, as its derivative is always negative, and that it is nonlinear in the molar fraction. In fact, as shown in example Figure 2, it is a sum of many decreasing hyperbolas each one defined by its own asymptote β i , hence it comprises of n þ 1 branches and exhibits n À 1 distinct roots. The only physically sound one is bounded in the [0, 1] range and it can be proved that the asymptotes which enclose that range are the ones corresponding to the maximum and to the minimum k-values, Beyond the obvious option of using the Newton-Raphson method to find the root, various alternative methods have been presented taking advantage of its special form to ensure safe and rapid convergence to the desired root [20,21].
Alternatively, the phase split problem can be treated as a constrained optimization problem where the system Gibbs energy in Eq. (42) needs to be minimized by varying k i under the mass balance constraint in Eq. (45).

Using k-values from correlations and charts
Equilibrium coefficients are functions of pressure, temperature and composition. However, at low pressures and temperatures, such as those prevailing at surface separators, the dependency on composition is very loose, thus allowing for the derivation of k-values correlations which only utilize pressure and temperature such as the one by Wilson [18].
. Similar correlations by Standing [22] and Whitson and Torp [23] have also been presented. An alternative approach is based on the utilization of charts which provide k-values at various pressures and temperatures. The generation of those charts is based on the observation that at high pressures k-values approach unity. In fact, there exists a composition dependent pressure value, known as the convergence pressure p k , at which all k-values become equal to unity. Charts for various convergence pressure values and system temperatures provide plots of the k-values as functions of pressure [24]. To utilize them in flash calculations, the user needs to determine the convergence pressure by means of any of the available methods [23,25,26] and select the appropriate chart where from the prevailing k-values can be obtained.
The solution algorithm is as follows The Rachford-Rice equation needs to be solved by means of any iterative function-solving method such as the Newton-Raphson one and the molar fraction update is given by where

Using composition dependent k-values from an EoS model
When an EoS model is available, components fugacity f i , hence fugacity coefficients φ i and k-values k i ¼ φ As mentioned above, the flash problem is governed by n þ 1 equations in n þ 1 unknowns. The problem can be further split to the solution of a system of n nonlinear equations (Eq. (43)) subject to one more nonlinear one (Eq. (45)). This way one needs to apply the Newton-Raphson method to solve the n nonlinear thermodynamic equilibrium equations and at each iteration compute β to ensure mass balance and composition consistency. To describe this algorithm, we define or equivalently, in a vector format: which needs to be driven to zero, i.e. gz , k ðÞ ¼ 0, by varying k i . The algorithm is identical to the SS one except step 6 which now reads.
6. If convergence has not been achieved, update the equilibrium coefficients by the Newton-Raphson method k k À J À1 gz , k ðÞ and return to step 2. The nxn Jacobian matrix is defined by The optimization approach uses any optimization method to minimize Gibbs energy subject to mass balance. Quasi-Newton methods such as the BFGS [27] only require computation of the Gibbs energy gradient with respect to the k-values, whereas a Newton method also requires the Hessian [27]. Hence, step 6 now reads.
6a. If convergence has not been achieved, compute the Gibbs energy gradient, update k-values by means of the BFGS method and return to step 2 or. 6b. If convergence has not been achieved, compute the Gibbs energy gradient and Hessian, update k-values by means of the Newton method and go to step 2.
The gradient and Hessian are defined by : (54) Although the Jacobian, gradient and Hessian formulae are rather complex to compute they allow for the very quick convergence of the optimization algorithm to its solution.

k-value initialization
Flash equations are always satisfied by a "trivial" solution which simply implies that Clearly, that solution satisfies mass balance and equilibrium conditions (Eq. (41) and (43)) and it also satisfies composition consistency (Eq. (44)) for any vapor phase molar fraction value. Converging to the physically sound rather than the trivial solution can only be ensured by utilizing appropriate initial estimates of the equilibrium coefficients. SS has proved to be more robust, yet slow, when initialized away from the true solution, as opposed to the Newton-Raphson, BFGS and Newton methods which perform rapidly only provided that they are initialized close to the solution.
To benefit from the advantages of each method most flash algorithms run a few  (47)) might be used. If a stability test has been run before the phase split, the k-values obtained can be used as a very good estimate of the final solution.
In flow applications where physical properties are obtained by EoS models, k-values are often initialized to the values they exhibited at the same point in the previous timestep, thus taking advantage of the fact that flow in petroleum engineering applications is a slow varying process with time. Even more accurate estimations can be obtained by extrapolating the converged k-values obtained in the previous 2 or 3 timesteps using linear or quadratic interpolation respectively [28].

Saturation condition calculations
The estimation of saturation pressure or temperature can be considered as a special case of a flash calculation where the molar ratio is known, i.e. β ¼ 0 for a bubble point or β ¼ 1 for a dew one, whereas pressure or temperature needs to estimated. The bubble or drop composition, known as incipient phase, needs to be estimated as well. At saturation conditions, the Rachford-Rice equation reads Equilibrium, i.e. equality of fugacity between the feed and the incipient phase, needs to be respected. Therefore, the n þ 1 equations that need to be solved for the case of a bubble point calculation are Eqs. (55) and (43) are the bubble composition y i ¼ z i k i , or equivalently the prevailing k-values, and the saturation pressure or temperature. For the case of a dew point calculation, the equations are (56) and (43) where y i ¼ z i and the drop composition is An alternative, more elegant approach is based on the fact that the TPD at a saturation point needs to be equal to zero. In other words, forming a bubble with the incipient phase composition, different than the feed one, retains the system Gibbs energy. A zero TPD value implies f y ðÞ When dealing with a dew point, i.e. Eq. (56), a similar result is obtained, Michelsen's algorithm [29] varies pressure until the following condition is met In detail, the algorithm is as follows 1. Initialize p sat to a pressure guaranteed to be in the two-phase region. This can be done by running a stability test at various pressures The Newton-Raphson derivative is given by fluid is physically single phase. They showed that the phase split equations can still be satisfied, this time with negative β values at pressures above the bubble point or with β values above unity at pressures above the upper or below the lower dew point. The more is the distance from the phase boundary the more is the absolute value of the molar fraction, eventually approaching À∞ and þ∞ at the convergence pressure p k . At convergence pressure, the equilibrium coefficients become equal to unity whereas beyond p k the flash equations have only one solution, the trivial one. Algorithms to compute the locus of the convergence pressure over a temperature range, known as "convergence locus" (CL), have been developed [31]. The negative flash area between the regular phase envelope and the CL is often referred to as the "shadow region" [32]. They also showed that stability tests can also be interpreted outside the phase envelope. Each of the two trial phases converges to a nontrivial solution (i.e. the TPD distance is positive) up to a locus in the shadow region, known as "stability test limit locus", STLL) which is enclosed by the CL. Such stability test results can be used to initialize negative flash calculations. Beyond STLL, the stability test only converges to the trivial solution. The regions discussed are shown in Figure 3 for a black oil, where the phase envelope interior is shown in red, cyan and yellow color and the latter corresponds to the spinodal. The shadow regions above the bubble point and the dew point lines are shown in pink and blue color respectively. Green color indicates the area outside the CL where the trivial solution is the only one to the phase split problem.
To interpret physically the results of a negative flash we firstly need to note that a molar fraction value of 0 < β < 1 in a regular flash calculation implies that β moles of gas of composition y i need to be added to 1 À β moles of liquid of composition x i to reconstruct the original feed composition z i . In a negative flash with β < 0, β jj¼ Àβ moles of gas need to be removed from 1 À β ¼ 1 þ β jjmoles of liquid to reconstruct one mole of the original feed composition. Similarly, when β > 1, β À 1 moles of liquid need to be removed from β moles of gas.
Clearly, negative flash solutions are not of any direct use in fluid flow calculations. However, they can significantly improve the convergence properties of the regular flash calculations close to the phase boundary by allowing the solution at some iteration to escape temporarily outside the phase envelope while trying to arrive to the exact solution.

Multiphase calculations
The need for multiphase calculations varies depending on the chemical engineering field. For example, in the upstream petroleum industry it is not that intense as multiphase equilibrium very rarely occurs in the reservoir and only when special studies in the wellbore and pipeline flow are considered. A case that is possible to happen in the reservoir is the presence of oil with high CO 2 content where two liquid phases (a CO 2 -rich and a CO 2 -poor one) and a vapor one could be formed. Things become more complicated when solids are considered as is the case with asphaltenes, waxes or hydrates. In the latter case, the phases that need to considered as possible to form are the solid one which may correspond to more than one hydrate structures (i.e. sI, sII and sH [33]), the aqueous phase which can be in liquid of solid form (ice) and the hydrocarbons phase (liquid, vapor or both). Nevertheless, multiphase equilibrium appears very often in chemical engineering processes taking place in process plants.
To identify such situations the standard approach is to repeatedly use the conventional two-phase Michelsen's stability test. Firstly, the test is run and if instability is detected then the vapor-liquid flash problem is solved. Subsequently, the equilibrium phase compositions are used as feeds (i.e. x i and/or y i instead of z i ) with suitable initial k-values to further detect if indeed they are stable or if one of them (e.g. the liquid one) will further split to two liquids.
Although many multiphase flash algorithms have been presented, the one developed by Michelsen is still considered as the most elegant one. By directly extending the two-phase flash requirements to a total of F phases, the mass balance, equilibrium and composition consistency expressions generalize to where β j denotes the molar fraction of phase 1 ≤ j ≤ F and y j i denotes the concentration of component 1 ≤ i ≤ n in phase 1 ≤ j ≤ F. Michelsen [34] proposed varying y j i and β j to minimize the objective function given by which satisfies Eq. (57) at its minimum. An alternative approach that combines stability and flash calculations in a single algorithm [35] at the cost of an increased set of variables that need to be determined, has also been presented. Unlike the previous algorithms, in the one presented here F denotes the maximum number of phases that might be present in equilibrium rather than the actual number of them. Upon convergence, this algorithm will also provide information about the presence or absence of each one of the potential phases.
The algorithm requires that one phase, surely known to be present in the mixture, is considered as the reference one, say phase r. This way, the equilibrium coefficients of any other potential phase can be defined with respect to the reference one, i.e. k j i ¼ y j i =y r i , where k r i ¼ 1. Let θ j be the stability variable of a phase, defined so that it is equal to zero when the phase is present (hence β j > 0) or exhibits a positive value when the phase does not exist (i.e. when β j ¼ 0). Therefore, β j > 0 and θ j ¼ 0 for an existing phase whereas β j ¼ 0 and θ j > 0 for a nonexisting one.
To solve the phase split problem we need to determine all k-values k j i , the molar fractions β j and the stability variables θ j for all phases but the reference one. Indeed, once those variables have been determined, the composition of any equilibrium phase can be computed by At the solution the mass balance and thermodynamic equilibrium conditions need to be simultaneously satisfied. For the first condition, the two-phase Rachford-Rice equation is extended to multiphase calculations as follows Note that the above equation needs to be satisfied for all k ¼ 1, ⋯, F and k 6 ¼ r. To satisfy the second condition, a minimum of the Gibbs energy is achieved when subject to β j ≥ 0, θ j ≥ 0 for all phases. Note that Eq. (63) is satisfied by definition for the reference phase, i.e. j ¼ r, as that phase is known to exist, hence θ r ¼ 0.
To solve the numerical problem it is initially assumed that all phases are present, hence all θ j are set to zero, the k-values are initialized using appropriate correlations or expected equilibrium phase compositions and molar fractions are equally spaced. Firstly, the mass balance and equilibrium equations are solved for the molar fractions and the stability variables using the currently estimates of the k-values. Subsequently, phase compositions and fugacities are computed using Eq. (61). Finally, k-values are updated in an inner loop by and calculations are repeated until convergence. It is interesting to note that for the case of VLE phase split calculations, by defining the liquid phase to be the reference one, Eq. (61) simplifies to Eq. (46). Furthermore, the extended Rachford-Rice equation reduces to When both phases are present, θ ¼ 0 and Eq. (65) simplifies to Eq. (45).

Accelerated phase behavior calculations
When flow simulations are considered, reliability undoubtedly comes first as lack of convergence or obtaining unrealistic results during the calculations at any grid block would lead to a general failure of the reservoir simulation run. However, some tolerance can be shown to the accuracy of the EoS model produced results due to the latter's inherent simplicity, to the nonexhaustive fluid's compositional analysis available and to questionable tuning procedures. In fact, small inaccuracies in the fluid behavior calculations that might be introduced can be partially remediated by the history matching procedure of the field model.
On the other hand, the ever increasing demand for complex flow domain models in terms of both grid and fluid models complexity has rendered nowadays the speed of phase behavior calculations as one of the most critical issues of flow simulation, especially for cases of complex thermodynamic phenomena such as near critical phase behavior and multiphase equilibrium in the presence of solids. As a result, speeding up phase behavior calculations is considered as a major issue, even if this involves some sacrifice in the calculations accuracy.

Rigorous methods
Reducing the number of components used to describe the fluid composition through a splitting and lumping procedure is the standard way to obtain simpler, hence faster EoS models. Firstly, the heavy end, usually corresponding to a limited carbon number, needs to be replaced with a large number of pseudo-components defined by means of computational methods. This way the flexibility during the EoS model tuning increases. The most pronounced method is the one developed by Whitson that utilizes the Gamma distribution [6]. Subsequently, the extended number of components is reduced (lumped) to a small number of pseudocomponents, usually 3 to 5, by means of algorithms which aim at preserving the EoS model's performance [7]. Finally, pure components are grouped together to minimize the composition vector size. Typical selections are N 2 with C 1 ,CO 2 with C 2 , nC 4 with iC 4 and nC 5 with iC 5 . When two or more components are lumped together, the new group's properties need to be rematched against the available PVT measurements. A very illustrative example is given my Ahmed [7] where a full C 7+ composition that includes N 2 and CO 2 , thus summing up to 11 components, reduces gradually the number of components to only 7 according to the lumping procedure shown in Table 3.
Other accelerating methods include different treatments of the mathematical form of the problem or of its variables [36,37] and utilizing solution acceleration techniques such as the GDEM update one [5]. Rasmussen et al. [32] provided criteria to completely skip phase behavior calculations during a simulation run when the prevailing equilibrium conditions fall within specific regions of the fluid's phase diagram. Simply speaking, if the fluid is a single phase one, most probably it will keep so if its distance to the phase boundary is large enough. So is the case with fluids lying well inside the two-phase region. In both cases, the stability test can be skipped whereas in the former one the phase split can be skipped as well.
Finally, efforts have been concentrated on utilizing advanced code optimization [38] and High Performance Computing (HPC) techniques which take advantage of  the parallel computing capabilities of modern computer architectures [39]. Despite the difficulties in distributing the work load and in optimizing memory transfer between clusters, impressive acceleration factors have been reported [40].

The reduced variable framework
Reduced variables methods are based on the fact that the intrinsic dimensionality of the stability and phase split calculations, hence the number of equations to be solved, is related to the rank of the complementary BIP matrix Γ, defined by γ ij ¼ 1 À k ij , rather than the number n of components used. Michelsen [41] derived the first reduced variables algorithm for cubic EoS models with zero BIPs (k ij ¼ 0) by showing that the equations to be solved could be reduced to only 3. Simply speaking, although the phase composition, e.g. y i , is a vector with n components, it is incorporated in the mixing rules only through its scalar projections to the components' a i and b i constants vectors, thus forming only two variables, i.e. a mix ¼ P ffiffiffiffi a i p y i ÀÁ 2 and b mix ¼ P b i y i . By further considering the molar fraction β, the number of variables to be determined reduces to only 3.
In general, the n þ 1 original variables (i.e. the k-values and the molar fraction) are replaced by a set of m þ 2 reduced ones, with m ≪ n, thus significantly reducing the phase behavior problem dimensionality. Several authors extended Michelsen's idea to calculations with nonzero BIP [42][43][44] by applying Singular Value Decomposition to the BIP matrix so as to split it in a sum of rank-1 matrices. The less is the number of rank-1 matrices required to reconstruct accurately the original BIP matrix, the less is the number of reduced variables that need to be utilized, hence the less is m. Nichita and Graciaa [45] presented an alternative reduced variables set which allows for an easier Hessian matrix computation procedure and faster convergence while Gaganis and Varotsis [46] proposed a new procedure for generating improved reduced variables.
More specifically, let the complementary BIP matrix Γ ¼ 1 À k ij ÈÉ be decomposed to a set of eigenvalues λ i and eigenvectors t i by use of the Singular Value Decomposition method [28], so that Γ ¼ P m i¼1 λ i t i t T i , where m denotes the rank of Γ. For the vapor phase, we define the projection vectors q i ¼ t i ∘ ffiffiffi a p and the reduced variables h V ¼ Q T y, where a ¼ a i fg is the vector containing the components energy parameters, Q ¼ q 1 ⋯ q m ÂÃ and operator ∘ denotes the Hadamard vector product (by element multiplication). The phase energy parameter a V and its derivative (required for the computation of phase fugacity) can be computed as functions of the reduced variables, that is fg . By further considering the vapor phase volume parameter b V as an unknown variable all required quantities (i.e. compressibility factor Z V from Eq. (7) and fugacity coefficients from Eq. (10)) can now be completed as functions of h V and b V .
The corresponding variables of the liquid phase can be easily computed by considering the vapor phase molar fraction β as an unknown variable and applying mass balance, i.e. h L ¼ Q T z À βh V ÀÁ = 1 À β ðÞ and b L ¼ b T z À βb V ÀÁ = 1 À β ðÞ , thus allowing for the computation of the liquid phase properties as well.
To summarize, h V , b V and β form an alternative set of variables in terms of which the phase split problem can be cast. The constraining equations that need to be satisfied are The solution algorithm is as follows 1. Initialize y i and β Solve the cubic polynomial for both phases (Eq. (8)) 7. Compute fugacity coefficients for both phases using (Eq. (11) Eqs. (66) and (67) guarantee thermodynamic equilibrium whereas Eq. (68) ensures mass balance. Clearly, the equations are nonlinear and their solution still requires the utilization of iterative function solving methods. Nevertheless, the benefit of the reduced variables approach lies in the cardinality of the variables set which is usually smaller than that of the conventional approach as it equals to m þ 2. When the BIP matrix contains many small or even zero values, as it is commonly the case with the EoS modeling of multicomponent fluids, the rank of matrix Γ is much smaller than its size (m ≪ n) which implies that the number of equations that need to be solved is significantly reduced. Moreover, reduced variables h i corresponding to very low eigenvalues λ i can also be neglected at the cost of the truncation error of matrix Γ. For the extreme case where all BIPs are equal to zero, m ¼ 1, Eq. (66) simplifies to a scalar one and only three nonlinear equations need to be solved regardless of the number of the mixture components [41]. Nevertheless, there has been some questioning about the real benefit of reduction methods as modern computers architecture has significantly reduced their computing time gain against the conventional ones [47,48].

Soft computing methods
Soft computing methods aim at solving phase equilibrium problems by utilizing data points rather than solving the thermodynamically rigorous equations discussed in the previous sections. Simply speaking, data related to the stability and phase split problems are generated and subsequently used to build correlations which provide directly the variables of interest such as the TPD value and the prevailing kvalues for the stability and phase split problems respectively. Such flow-specific and fluid-specific soft computing models are case dependent as they are generated using data obtained either prior to the specific simulation of interest or during that.
The benefit lies in that the generated correlations consist of simple, noniterative calculations which are by orders of magnitude faster than the conventional iterative ones. Although the numerical treatment of the datapoints involves purely numerical techniques such as regression, classification and clustering [49], thermodynamics are still incorporated indirectly in the soft computing based models as the data points used to build the models have been generated in advance by conventional rigorous methods.
Composition independent correlations to estimate the equilibrium coefficients (k-values), such as those of Standing and Whitson as well as the convergence pressure method, all discussed in 3.2.1, can be thought of as the simplest soft computing method to treat the phase split problem as they provide k-values estimates without being based on a rigorous EoS model, hence avoiding the iterative solution of the fugacity equations or the minimization of the Gibbs energy.
Voskov and Tchelepi [50] proposed the generation and storage of the encountered tie-lines in Tables "on the fly". Initially, for each feed composition encountered during the simulation, the phase split problem is solved conventionally and the equilibrium compositions (i.e. the tie line endpoints) are stored. For each subsequent feed the algorithm searches quickly the Tables to identify the closest stored tie-lines and interpolate them linearly to get the equilibrium compositions. If no close enough tie lines can be found, the phase split problem is solved conventionally, and the table is enriched. Stability is determined by using the negative flash approach [30]. To reduce the computing time cost for accessing and further building-up the tie line Table, Belkadi et al. [51] proposed the Tie-line Distance Based Approximation which further accelerates the search procedure.
Gaganis and Varotsis [52,53] presented the methodology to develop proxy models for treating both the phase stability and phase split problems using machine learning tools. Their approach aims at solving conventionally the phase behavior problem for a set of sampled operating points and using the obtained data to generate explicit proxy models using multivariate regression models such as neural networks to directly predict the prevailing equilibrium coefficients values given feed composition, pressure and temperature (for nonisothermal runs). For the phase stability problem, their model outputs a positive nonlinear transformation of the conventional TPD value that exhibits the same sign as the former (Figure 4). Their model utilizes Support Vector Machines, SVM [54] to provide the same binary stable/unstable answers anywhere in the operating space even outside the stability test limit locus [31]. An improved stability test method has been presented by Gaganis [55] which reliefs the need to model accurately the phase boundary thus allowing for even simpler and faster to evaluate stability models. His approach develops two classifiers which only identify whether the point under question lies "far enough" from the phase boundary or not. If it lies far enough outside of the phase envelope, then the fluid is surely single phase whereas it is certainly at twophase when lying well inside the phase envelope. If a certain answer cannot be obtained, a regular stability algorithm is invoked.

Conclusions
Equations of State of varying complexity and accuracy are nowadays available to describe the thermodynamic behavior of almost all types of fluids. Beyond the classic and easy-to-implement cubic EoS models, recent advances in perturbation theory have allowed its application to the derivation of models that describe accurately in a microscopic level the behavior of fluids.
Phase behavior calculations by means of EoS models are massively required during all types of flow simulations, thus rendering the availability of robust, thermodynamic rigorous algorithms as of major importance. However, as the required computational load can be very heavy, various accelerating methods have been developed, and they have been proved to perform very well.

Author details
Vassilis Gaganis 1,2 1 National Technical University of Athens, Hellas, Athens, Greece 2 Foundation for Research and Technology, Hellas, Athens, Greece *Address all correspondence to: vgaganis@metal.ntua.gr © 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.