**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. *xi* and/or *yi* instead 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 *F* phases, the mass balance, equilibrium and composition consistency expressions generalize to

$$\begin{aligned} &\sum\_{j=1}^{F} \beta\_j \mathbf{y}\_i^j = \mathbf{z}\_i \ \mathbf{i} = \mathbf{1}, \dots, n\\ &f\_i^{(\mathbf{y}\_1)} = f\_i^{(\mathbf{y}\_2)} = \dots = f\_i^{(\mathbf{y}\_F)} \Leftrightarrow \mathbf{y}\_i^1 \boldsymbol{\phi}\_i^{(\mathbf{y}\_1)} = \mathbf{y}\_i^2 \boldsymbol{\phi}\_i^{(\mathbf{y}\_2)} = \dots = \mathbf{y}\_i^F \boldsymbol{\phi}\_i^{(\mathbf{y}\_F)} \ i = \mathbf{1}, \dots, n \end{aligned} \tag{59}$$

where *<sup>β</sup> <sup>j</sup>* denotes the molar fraction of phase 1≤*j*<sup>≤</sup> *<sup>F</sup>* and *<sup>y</sup> <sup>j</sup> <sup>i</sup>* denotes the concentration of component 1≤*i* ≤*n* in phase 1≤*j*≤ *F*. Michelsen [34] proposed varying *y <sup>j</sup> <sup>i</sup>* and *β <sup>j</sup>* to minimize the objective function given by

$$Q = \sum\_{j=1}^{F} \beta\_j - \sum\_{i=1}^{n} z\_i \sum\_{k=1}^{F} \frac{\beta\_k}{\phi\_i^{(\mathbf{y}\_k)}},\tag{60}$$

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 <sup>j</sup> <sup>i</sup>* <sup>¼</sup> *<sup>y</sup> <sup>j</sup> <sup>i</sup> =y<sup>r</sup> <sup>i</sup>* , where *kr <sup>i</sup>* ¼ 1. Let *θ <sup>j</sup>* be the stability variable of a phase, defined so that it is equal to zero when the phase is present (hence *β <sup>j</sup>* >0) or

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

exhibits a positive value when the phase does not exist (i.e. when *β <sup>j</sup>* ¼ 0). Therefore, *β <sup>j</sup>* >0 and *θ <sup>j</sup>* ¼ 0 for an existing phase whereas *β <sup>j</sup>* ¼ 0 and *θ <sup>j</sup>* >0 for a nonexisting one.

To solve the phase split problem we need to determine all k-values *k <sup>j</sup> <sup>i</sup>* , the molar fractions *β <sup>j</sup>* and the stability variables *θ <sup>j</sup>* for all phases but the reference one. Indeed, once those variables have been determined, the composition of any equilibrium phase can be computed by

$$y\_i^j = \frac{z\_i}{1 + \sum\_{\substack{j=1 \\ j \neq r}}^r \beta\_j \left(k\_i^j \mathbf{e}^{\theta\_j} - 1\right)}, \ i = 1, \dots, n, \ j = 1, \dots, F. \tag{61}$$

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

$$r\_k(\mathfrak{g}, \mathfrak{g}) = \sum\_{i=1}^n \frac{z\_i \left(k\_i^k \mathbf{e}^{\boldsymbol{\theta}\_k} - \mathbf{1}\right)}{\mathbf{1} + \sum\_{j=1}^F \beta\_j \left(k\_j^k \mathbf{e}^{\boldsymbol{\theta}\_j} - \mathbf{1}\right)}, \ j = \mathbf{1}, \dots, F. \tag{62}$$

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

$$
\beta\_j \theta\_j = \mathbf{0},
\tag{63}
$$

subject to *β <sup>j</sup>* ≥ 0, *θ <sup>j</sup>* ≥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 *θ<sup>r</sup>* ¼ 0.

To solve the numerical problem it is initially assumed that all phases are present, hence all *θ <sup>j</sup>* 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

$$k\_i^j = \phi\_i^{(r)} / \phi\_i^{(j)},\tag{64}$$

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

$$r(\boldsymbol{\beta}, \theta) = \frac{z\_i(k\_i \mathbf{e}^{\theta} - \mathbf{1})}{\mathbf{1} + \beta(k\_i \mathbf{e}^{\theta} - \mathbf{1})}. \tag{65}$$

When both phases are present, *θ* ¼ 0 and Eq. (65) simplifies to Eq. (45).
