*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 P *ln f* ð Þ*<sup>y</sup> <sup>i</sup> = f* ð Þ *x i* � �<sup>2</sup> becomes 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. *β* ¼ 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

$$\sum\_{i=1}^{n} z\_i k\_i - \mathbf{1} = \mathbf{0} \text{ for } p\_b \tag{55}$$

$$\sum\_{i=1}^{n} \mathbf{z}\_i / k\_i - \mathbf{1} = \mathbf{0} \quad \text{for } p\_d \,. \tag{56}$$

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) where *xi* ¼ *zi*. The *n* þ 1 unknowns 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* ¼ *zi* and 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 *f* ð Þ*y <sup>i</sup>* ¼ *f* ð Þ*z <sup>i</sup>* , hence *Yi* <sup>¼</sup> *ziki* <sup>¼</sup> *ziφ*ð Þ*<sup>z</sup> <sup>i</sup> <sup>=</sup>φ*ð Þ*<sup>y</sup> i* ¼ *yi f* ð Þ*z <sup>i</sup> = f* ð Þ*y <sup>i</sup>* which in turn implies <sup>P</sup>*Yi* <sup>¼</sup> <sup>P</sup>*yi* <sup>¼</sup> 1. When dealing with a dew point, i.e. Eq. (56), a similar result is obtained, <sup>P</sup>*Xi* <sup>¼</sup> <sup>P</sup>*xi* <sup>¼</sup> 1. Michelsen's algorithm [29] varies pressure until the following condition is met

$$Q(p, k\_i) = \mathbf{1} - \sum Y\_i = \mathbf{0} \text{ or } Q(p, k\_i) = \mathbf{1} - \sum X\_i = \mathbf{0}.\tag{57}$$

In detail, the algorithm is as follows


10.Check trivial solution by evaluating P *ln yi =zi* � �<sup>2</sup> <sup>&</sup>lt;*<sup>δ</sup>* or <sup>P</sup> *ln <sup>χ</sup>* ð Þ *<sup>i</sup>=zi* <sup>2</sup> < *δ*

The Newton-Raphson derivative is given by

$$\begin{aligned} \frac{\partial \mathbf{Q}}{\partial p} &= \sum\_{i=1}^{n} \mathcal{y}\_{i} \frac{f\_{i}^{(x)}}{f\_{i}^{(y)}} \left( \frac{\partial f\_{i}^{(y)}}{\partial p} \frac{\mathbf{1}\_{(y)}}{f\_{i}^{(y)}} - \frac{\partial\_{i}^{(x)}}{\partial p} \frac{\mathbf{1}\_{(x)}}{f\_{i}^{(x)}} \right) & \quad \text{for } p\_{b} \\\ \frac{\partial \mathbf{Q}}{\partial p} &= \sum\_{i=1}^{n} \mathcal{x}\_{i} \frac{f\_{i}^{(x)}}{f\_{i}^{(x)}} \left( \frac{\partial f\_{i}^{(x)}}{\partial p} \frac{\mathbf{1}\_{(x)}}{f\_{i}^{(x)}} - \frac{\partial f\_{i}^{(x)}}{\partial p} \frac{\mathbf{1}\_{(x)}}{f\_{i}^{(x)}} \right) & \quad \text{for } p\_{d} \end{aligned} \tag{58}$$

#### **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
