*Perturbation Theory and Phase Behavior Calculations Using Equation of State Models DOI: http://dx.doi.org/10.5772/intechopen.93736*

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 *yi* itself. The algorithm is as follows

1.Compute fugacity of each component of the feed *f* ð*z*Þ *<sup>i</sup>* using the EoS model


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


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


**Table 2.** *Stability test result selection.* *Perturbation Theory and Phase Behavior Calculations Using Equation of State Models DOI: http://dx.doi.org/10.5772/intechopen.93736*

**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, *yi* and *xi* 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.

$$\mathbf{z}\_{i} = (\mathbf{1} - \boldsymbol{\beta})\mathbf{x}\_{i} + \boldsymbol{\beta}\mathbf{y}\_{i} \quad i = \mathbf{1}, \ldots, n. \tag{41}$$

The second condition additionally requires the phase compositions to be so that the two-phase system's Gibbs energy, defined by

$$G = (\mathbf{1} - \boldsymbol{\beta}) \sum\_{i=1}^{n} \mathbf{x}\_i \ln f\_i^{(\mathbf{x})} + \boldsymbol{\beta} \sum\_{i=1}^{n} \mathbf{y}\_i \ln f\_i^{(\mathbf{y})},\tag{42}$$

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.

$$f\_i^{(x)} - f\_i^{(y)} = 0 \Rightarrow \phi\_i^{(y)} y\_i p - \phi\_i^{(x)} \ge\_i p = 0 \Rightarrow \frac{\phi\_i^{(x)}}{\phi\_i^{(y)}} = \frac{y\_i}{\varkappa\_i} = k\_i, \quad i = 1, \ldots, n. \tag{43}$$

Finally, we need to ensure that the composition of each equilibrium phase is consistent by summing up to unity. Equivalently, we may require that

$$\sum\_{i=1}^{n} (\mathbf{x}\_i - \mathbf{y}\_i) = \mathbf{0}.\tag{44}$$

Summarizing, the solution of the flash problem can be seen as the solution of a system of 2*n* þ 1 equations, that is Eq. (41), (43) and (44), in 2*n* þ 1 unknowns, i.e. *yi* , *xi* and *β*.

### *A Collection of Papers on Chaos Theory and Its Applications*

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

$$r(\beta) = \sum\_{i=1}^{n} \frac{z\_i}{\beta - \overline{\beta}\_i} = \mathbf{0},\tag{45}$$

where *β<sup>i</sup>* ¼ 1*=*ð Þ 1 � *ki* , which can be solved for the molar fraction *β*. Given *β* and the k-values, the equilibrium phase compositions can then be obtained by

$$\mathbf{x}\_{i} = \frac{1}{k\_{i} - 1} \frac{x\_{i}}{\beta - \overline{\beta}\_{i}}, \quad \mathbf{y}\_{i} = k\_{i} \mathbf{x}\_{i}, \quad i = 1, \ldots, n \tag{46}$$

Therefore, the phase split problem can be treated as the solution of a system of *n* þ 1 equations, that is Eq. (43) and (45), in *n* þ 1 unknowns, i.e. *ki* and *β*. Of course, if the kvalues 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* þ 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, i.e. [*βmin* ¼ 1*=*ð Þ 1 � *kmax* , *βmax* ¼ 1*=*ð Þ 1 � *kmin* ]. 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 *ki* under 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

**Figure 2.** *The Rachford-Rice equation and its asymptotes.*

*Perturbation Theory and Phase Behavior Calculations Using Equation of State Models DOI: http://dx.doi.org/10.5772/intechopen.93736*

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

$$k\_i(p, T) = \frac{\exp\left(5.37(1 + \alpha\_l)(1 - T\_{c\_i}/T)\right)}{p/p\_{c\_i}}.\tag{47}$$

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


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

$$
\beta \leftarrow \beta - \frac{dr}{d\beta}\Big|\_{\beta} r(\beta),
\tag{48}
$$

where

$$\frac{dr}{d\beta} = -\sum\_{i=1}^{n} \frac{z\_i}{\left(\beta - \overline{\beta}\_i\right)^2}.\tag{49}$$

#### *4.2.2 Using composition dependent k-values from an EoS model*

When an EoS model is available, components fugacity *fi* , hence fugacity coefficients *<sup>φ</sup><sup>i</sup>* and k-values *ki* <sup>¼</sup> *<sup>φ</sup>*ð Þ *<sup>x</sup> <sup>i</sup> <sup>=</sup>φ*ð Þ*<sup>y</sup> <sup>i</sup>* can 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*


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

$$\mathbf{g}\_i = f\_i^{(y)} - f\_i^{(x)},\tag{50}$$

or equivalently, in a vector format:

$$\mathbf{g(z,k)} = \mathbf{f^{(y)}} - \mathbf{f^{(y)}},\tag{51}$$

which needs to be driven to zero, i.e. **g z**ð Þ¼ , **k 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 **k k** � **J** �1 **g z**ð Þ , **k** and return to step 2. The *nxn* Jacobian matrix is defined by

$$\mathbf{J} = \frac{\partial \mathbf{g}(\mathbf{z}, \mathbf{k})}{\partial \mathbf{k}} = \left\{ \frac{\partial \mathbf{g}\_i}{\partial k\_j} \right\} = \left\{ \frac{\partial f\_i^{(x)}}{\partial k\_j} - \frac{\partial f\_i^{(y)}}{\partial k\_j} \right\}. \tag{52}$$

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.


*Perturbation Theory and Phase Behavior Calculations Using Equation of State Models DOI: http://dx.doi.org/10.5772/intechopen.93736*

The gradient and Hessian are defined by

$$\frac{\partial G(\mathbf{z}, \mathbf{k})}{\partial \mathbf{k}} = \begin{Bmatrix} \partial G \\ \frac{\partial k\_i}{\partial k\_i} \end{Bmatrix} \tag{53}$$

$$\frac{\partial^2 G(\mathbf{z}, \mathbf{k})}{\partial \mathbf{k} \partial \mathbf{k}^T} = \left\{ \frac{\partial^2 G}{\partial k\_i \partial k\_j} \right\} = \left\{ \frac{\partial^2}{\partial k\_i \partial k\_j} \left( (\mathbf{1} - \boldsymbol{\beta}) \sum\_{i=1}^n \boldsymbol{x}\_i \ln f\_i^{(\mathbf{x})} + \boldsymbol{\beta} \sum\_{i=1}^n \boldsymbol{y}\_i \ln f\_i^{(\mathbf{y})} \right) \right\}. \tag{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.
