Open access peer-reviewed chapter - ONLINE FIRST

Perturbation Theory and Phase Behavior Calculations Using Equation of State Models

By Vassilis Gaganis

Submitted: February 22nd 2020Reviewed: August 25th 2020Published: October 12th 2020

DOI: 10.5772/intechopen.93736

Downloaded: 14


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.


  • statistical mechanics
  • perturbation theory
  • equation of state
  • phase behavior
  • phase stability
  • phase split
  • algorithms

1. 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 fpTvm=0which relate molar volume vmto 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.

2. The PR and SRK cubic EoS models

2.1 Development of the cubic EoS models

The ideal gas law pvm=RT, where the gas constant R=kBNAis 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 pvm=ZRTcan 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 aterm in Eq. (1) can be thought of as a term accounting for the attractive forces between molecules as it reduces pressure. Parameter baccounts for the molecules volume which becomes significant at high pressures (i.e. liquid state) as limpvm=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

PR1 + √21-√20.45724αR2Tc2/pc0.07780RTc/pc

Table 1.

Cubic EoS models constants.


where mis a function of the component acentric factor ωdefined by


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-C8 and pseudo-C11 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 aand b. For a mixture of known composition zi, they are given by


The Binary Interaction Parameters (BIP) kijaccount 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 Zcof 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=pvm/RTrather than for vmis preferred


where the dimensionless EoS constants are given by


2.2 Use of the cubic EoS models

As soon as the EoS constants have been defined, the compressibility factor Zcan 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 Mdenotes the fluid molar mass. Components fugacity coefficients φi, hence fugacity fi=φizip, can be computed by the following expressions


where Ai=aip/RT2, Bi=bip/RTand A/zi=j=1nzj1kijAiAj. Derivative properties such as the Joule-Thomson coefficient μJTcan be computed by differentiating the EoS and incorporating the derivatives in the rigorous thermodynamic definitions of the properties. For example,


2.3 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 vmEoSby some amount that depends on the fluid composition and its components properties. More specifically, the shifted volume is given by


Parameters ciare component specific and they are usually given as functions of the covolume parameters bi, that is


where values of sifor 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.

3. 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 urand the pair correlation function grrespectively, both functions of the distance raway 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 urand gris 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 urand grfunctions 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].

3.1 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 Eicorresponds to the energy of each possible microstate, β=1/kBTis the thermodynamic beta and kBis the Boltzmann constant. Note that Qis 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 QNVTat given volume and temperature is given by [3].




and ZNVTis known as the configuration integral. It is easy to show that if the system potential energy UNis assumed to be zero then the configuration integral ZNsimplifies to the system volume and the application of the related partition function leads simply to the ideal gas law. On the other hand, when UN0, it is often represented by a sum of pair-wise potentials, i.e.


Various pair potential models urhave 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, urand 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:


In the equations above σ is the sphere diameter, parameter γis used to scale the well width, lis 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.

3.2 Derivation of fluid properties for specific pair functions

By selecting the pair potential model urand incorporating it the configuration integral ZNand 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


Therefore, internal energy can be obtained as a function of the particle properties urand 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 ZN, 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


3.3 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


3.4 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, U0and Uprespectively, 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 U0and Upwhereas the original fluid energy is obtained for λ=1. By replacing that expression to the configuration integral we obtain


where the .operator 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




The treatment above has allowed the energy to be described by simpler expressions of the perturbation energy term based on the c1,c2,c3parameters. Therefore, to get the full energy expression one needs to choose the Upmodel, compute the values of the c1,c2,c3parameters 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 c1,c2,c3parameters, 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: Up=ij>iuprij.

  2. The pair potentials are equal between any pair of molecules: Up=ij>iuprij=NN12upr12.

By introducing those assumptions to Eq. (37) the calculation of coefficient c1simplifies to


To proceed we further need to introduce the following assumptions:

  1. The reference system to describe U0is the hard-sphere one

  2. The particles are uniformly distributed which implies that the pair correlation function g0ris equal to one at any distance beyond the limits of the particle and equal to zero inside that

  3. The free fluid volume is VNbwhere bis the particle volume that equals to b=2/3πσ3

and we end up with


Finally, by utilizing first order approximation only (up to c1), replacing c1in the free energy Eq. (36) and differentiating over volume to obtain pressure (i.e. p=H/VT) the well known Van der Waals equation is obtained


Interestingly, the perturbed energy term uphas not been defined explicitly but it has been incorporated into the EoS aparameter. From a perturbation theory point of view, the accuracy of the Van der Waals equation of state can be improved by further considering the c2term, the calculation of which, however, is quite more complicated.

4. Conventional phase behavior calculations

4.1 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 ki=yi/zi, also known as distribution coefficients, rather than the bubble composition yiitself. The algorithm is as follows

  1. Compute fugacity of each component of the feed fiz)using the EoS model

  2. Initialize kiusing Wilson’s correlation [18]

  3. Assume feed is a liquid and look for a bubble, i.e. compute Yi=kizi

  4. Compute trial bubble composition sum SV=Yi

  5. Normalize composition yi=Yi/SVand compute its fugacity fiy)

  6. Compute correction factor Ri=1/SVfiz)/fiy)

  7. Check for convergence by evaluating Ri12<ε

  8. If convergence has not been achieved, update the equilibrium coefficients by applying kikiRi

  9. After convergence has been achieved check if the algorithm has arrived at a trivial solution by evaluating lnki2<δ

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 Xi=zi/ki

  2. Compute drop composition sum SL=Xi

  3. Normalize composition xi=Xi/SLand compute its fugacity fix)

  4. Compute correction factor Ri=1/SLfix)/fiz)

  5. Check for convergence by evaluating Ri12<ε

  6. If convergence has not been achieved, update the equilibrium coefficients by applying kikiRi

  7. After convergence has been achieved check if the algorithm has converged to a trivial solution by evaluating lnki2<δ

As soon as both calculations have been completed, Table 2 can be used to reckon on the phase state.

Vapor phase testLiquid phase testResult
Trivial solutionTrivial solutionStable
SV ≤ 1Trivial solutionStable
Trivial solutionSL ≤ 1Stable
SV ≤ 1SL ≤ 1Stable
SV > 1Trivial solutionUnstable
Trivial solutionSL > 1Unstable
SV > 1SL > 1Unstable
SV > 1SL ≤ 1Unstable
SV ≤ 1SL > 1Unstable

Table 2.

Stability test result selection.

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).

Figure 1.

The spinodal.

4.2 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, yiand xirespectively, 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+1equations, that is Eq. (41), (43) and (44), in 2n+1unknowns, i.e. yi, xiand β.

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/1ki, 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+1equations, that is Eq. (43) and (45), in n+1unknowns, i.e. kiand β. Of course, if the k-values are known by any means, the problem simplifies to the solution of the Rachford-Rice Eq. (45) to compute βand phase compositions are obtained from Eq. (46).

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+1branches and exhibits n1distinct 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, i.e. [β¯min=1/1kmax, β¯max=1/1kmin]. 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].

Figure 2.

The Rachford-Rice equation and its asymptotes.

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 kiunder the mass balance constraint in Eq. (45).

4.2.1 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 pk, 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

  1. Estimate convergence pressure pk

  2. Get kifrom convergence pressure-based correlations or Tables

  3. Solve the Rachford-Rice equation (Eq. (45)) for the vapor phase molar fraction.

  4. Compute phase compositions using Eq. (46).

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




4.2.2 Using composition dependent k-values from an EoS model

When an EoS model is available, components fugacity fi, hence fugacity coefficients φiand k-values ki=φix/φiycan be accurately computed rather than been read from charts. Apart from the nonlinearity of the Rachford-Rice equation (Eq. (45)), the complex formulae (Eq. (11)) relating phase composition to fugacity through the EoS model introduces additional nonlinearity to the calculation of the k-values thus imposing the need for iterative solution methods.

Computations may involve any one of the three methods available, i.e. Successive Substitution (SS), numerical solution of the systems of equations in Eq. (43) and (45) by means of the Newton Raphson method or direct minimization of the system Gibbs energy in Eq. (42) by means of optimization algorithms.

The SS method starts with an estimation of the k-values, solves the Rachford-Rice equation for βto ensure mass balance and computes the phase composition and components fugacity. If phase fugacities are not equal, k-values are updated by the inverse fugacity coefficient ratio in Eq. (43). The algorithm is as follows

  1. Initialize ki

  2. Solve the Rachford-Rice equation (Eq. (45)) for the molar fraction β

  3. Compute phase compositions using Eq. (46)

  4. Solve the cubic polynomial of each phase (Eq. (8)) and compute components fugacity (Eq. (11))

  5. Check for convergence by evaluating lnfiy/fix2<ε

  6. If convergence has not been achieved, update the equilibrium coefficients by applying Eq. (43), i.e. kikiφix/φiy, and return to step 2

As mentioned above, the flash problem is governed by n+1equations in n+1unknowns. The problem can be further split to the solution of a system of nnonlinear equations (Eq. (43)) subject to one more nonlinear one (Eq. (45)). This way one needs to apply the Newton-Raphson method to solve the nnonlinear 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. gzk=0, by varying ki. 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 kkJ1gzkand return to step 2. The nxnJacobian 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


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.

4.2.3 k-value initialization

Flash equations are always satisfied by a “trivial” solution which simply implies that xi=yi=zi, hence ki=1. 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 SS iterations until the convergence criterion lnfiy/fix2becomes sufficiently small. Then the algorithm switches to any other method that converges rapidly to the solution. To initialize SS, Wilson’s correlation (Eq. (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].

4.3 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. β=0for a bubble point or β=1for 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+1equations that need to be solved for the case of a bubble point calculation are Eqs. (55) and (43) where xi=zi. The n+1unknowns are the bubble composition yi=ziki, 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 yi=ziand the drop composition is xi=zi/ki.

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 fiy=fiz, hence Yi=ziki=ziφiz/φiy=yifiz/fiywhich in turn implies Yi=yi=1. When dealing with a dew point, i.e. Eq. (56), a similar result is obtained, Xi=xi=1. Michelsen’s algorithm [29] varies pressure until the following condition is met


In detail, the algorithm is as follows

  1. Initialize psatto a pressure guaranteed to be in the two-phase region. This can be done by running a stability test at various pressures

  2. Initialize ki

  3. Compute Yi=zikifor a bubble point or Xi=zi/kifor a dew point

  4. Compute SV=Yior SL=Xi

  5. Normalize incipient phase composition using yi=Yi/SVor xi=Xi/SL

  6. Compute incipient phase fugacity fiyor fix

  7. Update incipient phase composition using Yi=yifiz/fiyor Xi=xifiz/fix

  8. Update pressure by running a Newton-Raphson iteration ppQQ/p

  9. Check for convergence by evaluating lnfiy/fiz2<εor lnfiχ/fiz2<ε

  10. Check trivial solution by evaluating lnyi/zi2<δor lnχi/zi2<δ

The Newton-Raphson derivative is given by


4.4 Negative flash calculations

Whitson and Michelsen [30] extended the regular phase split algorithm beyond the limits of the phase envelope to allow flash calculations at conditions where the 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 pk. At convergence pressure, the equilibrium coefficients become equal to unity whereas beyond pkthe 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.

Figure 3.

Regular phase envelope, shadow region and trivial solution region.

To interpret physically the results of a negative flash we firstly need to note that a molar fraction value of 0<β<1in a regular flash calculation implies that βmoles of gas of composition yineed to be added to 1βmoles of liquid of composition xito reconstruct the original feed composition zi. In a negative flash with β<0, β=βmoles of gas need to be removed from 1β=1+βmoles of liquid to reconstruct one mole of the original feed composition. Similarly, when β>1, β1moles 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.

4.5 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 CO2 content where two liquid phases (a CO2-rich and a CO2-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. xiand/or yiinstead of zi) 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 Fphases, the mass balance, equilibrium and composition consistency expressions generalize to


where βjdenotes the molar fraction of phase 1jFand yijdenotes the concentration of component 1inin phase 1jF. Michelsen [34] proposed varying yijand βjto 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 Fdenotes 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. kij=yij/yir, where kir=1. Let θjbe 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>0and θj=0for an existing phase whereas βj=0and θj>0for a nonexisting one.

To solve the phase split problem we need to determine all k-values kij, the molar fractions βjand the stability variables θjfor 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,,Fand kr. To satisfy the second condition, a minimum of the Gibbs energy is achieved when


subject to βj0, θj0for 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 θjare 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, θ=0and Eq. (65) simplifies to Eq. (45).

5. 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.

5.1 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 pseudo-components, 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 N2 with C1, CO2 with C2, nC4 with iC4 and nC5 with iC5. 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 C7+ composition that includes N2 and CO2, thus summing up to 11 components, reduces gradually the number of components to only 7 according to the lumping procedure shown in Table 3.

Original components set
Lumped components set
N2 + C1CO2 + C2C3 + iC4 + nC4iC5 + nC5 + C6F1F2F3

Table 3.

Components’ number reduction by splitting and lumping.

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].

5.2 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=1kij, rather than the number nof components used. Michelsen [41] derived the first reduced variables algorithm for cubic EoS models with zero BIPs (kij=0) by showing that the equations to be solved could be reduced to only 3. Simply speaking, although the phase composition, e.g. yi, is a vector with ncomponents, it is incorporated in the mixing rules only through its scalar projections to the components’ aiand biconstants vectors, thus forming only two variables, i.e. amix=aiyi2and bmix=biyi. By further considering the molar fraction β, the number of variables to be determined reduces to only 3.

In general, the n+1original variables (i.e. the k-values and the molar fraction) are replaced by a set of m+2reduced ones, with mn, 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 Γ=1kijbe decomposed to a set of eigenvalues λiand eigenvectors tiby use of the Singular Value Decomposition method [28], so that Γ=i=1mλititiT, where mdenotes the rank of Γ. For the vapor phase, we define the projection vectors qi=tiaand the reduced variables hV=QTy, where a=aiis the vector containing the components energy parameters, Q=q1qmand operator denotes the Hadamard vector product (by element multiplication). The phase energy parameter aVand its derivative (required for the computation of phase fugacity) can be computed as functions of the reduced variables, that is aV=hVTΛhVand aV/y=2hV, where Λ=diagλ1λm. By further considering the vapor phase volume parameter bVas an unknown variable all required quantities (i.e. compressibility factor ZVfrom Eq. (7) and fugacity coefficients from Eq. (10)) can now be completed as functions of hVand bV.

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. hL=QTzβhV/1βand bL=bTzβbV/1β, thus allowing for the computation of the liquid phase properties as well.

To summarize, hV, bVand β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 yiand β

  2. Compute hV=QTy

  3. Compute aV=hVTΛhV, aV/y=2hVand bV=bTy

  4. Compute hL=QTzβhV/1βand bL=bTzβbV/1β

  5. Compute aL=hLTΛhL, aL/x=2hLand bL=bTx

  6. Solve the cubic polynomial for both phases (Eq. (8))

  7. Compute fugacity coefficients for both phases using (Eq. (11))

  8. Compute ki=φix/φiyand phase compositions using (Eq. (46))

  9. Check convergence by evaluating if Eq. (66) are satisfied

  10. If convergence has not been achieved, update hV, bVand βby means of a Newton-Raphson step and return to step 3.

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 (mn) which implies that the number of equations that need to be solved is significantly reduced. Moreover, reduced variables hicorresponding to very low eigenvalues λican 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].

5.3 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 k-values 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 two-phase when lying well inside the phase envelope. If a certain answer cannot be obtained, a regular stability algorithm is invoked.

Figure 4.

The SVM output equals to zero at the phase boundary.

6. 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.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Vassilis Gaganis (October 12th 2020). Perturbation Theory and Phase Behavior Calculations Using Equation of State Models [Online First], IntechOpen, DOI: 10.5772/intechopen.93736. Available from:

chapter statistics

14total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us