Published on *PNG 520: Phase Behavior of Natural Gas and Condensate Fluids* (https://www.e-education.psu.edu/png520)

Module Goal: To establish the basic framework for vapor-liquid equilibrium calculations.

Module Objective: To formulate the basic governing equations for performing flash calculations.

In a previous module we derived two different objective functions for the purpose of solving the flash equilibrium problem. Let us take a closer look at these equations:

$${F}_{y}\left({\alpha}_{g}\right)={\displaystyle \sum}_{i=1}^{n}\frac{{z}_{i}{K}_{i}}{1+{\alpha}_{g}\left({K}_{i}-1\right)}-1=0$$ (13.1a)

$${F}_{x}\left({\alpha}_{g}\right)={\displaystyle \sum}_{i=1}^{n}\frac{{z}_{i}}{1+{\alpha}_{g}\left({K}_{i}-1\right)}-1=0$$ (13.1b)

These equations arise from simple mole balances and the concept of equilibrium ratios. As we discussed before, if we are given {z_{i}; i = 1,2,…,n} and, if for some reason, let us say, we are able to obtain all the equilibrium ratios {K_{i}; i = 1,2,…,n}; the only unknown in these objective functions would be the vapor fraction ‘${\alpha}_{g}$’.

Once we are able to solve for this ${\alpha}_{g}$, we will have no problem applying any combination of equations (12.5), (12.7), and (12.11) to solve for all the vapor and liquid compositions at equilibrium {y_{i}, x_{i}}:

$${K}_{i}=\frac{{y}_{i}}{{x}_{i}}$$ (12.5)

$${y}_{i}=\frac{{z}_{i}{K}_{i}}{1+{\alpha}_{g}\left({K}_{i}-1\right)}$$ (12.7)

$${x}_{i}=\frac{{z}_{i}}{1+{\alpha}_{g}\left({K}_{i}-1\right)}$$ (12.11)

With all this information, the VLE problem would be completely solved. With the compositional information and the use of suitable equations of state and correlations, we can then find all other related properties, such as densities, viscosities, molecular weights, etc. This is why we call these the *objective functions* once we have solved them, we have achieved the objective of the VLE calculation.

Of course one question remains unanswered — how do you solve for that ‘$\alpha $ _{g}’, which is buried inside those expressions? And secondly, which of the two objective functions that are available to us shall we use? It turns out that the answers to these questions are tied to one another. In fact, the proper* objective function* that we shall use to solve for ‘$\alpha $_{g}’ is the same equation that simplifies the process of solving for that unknown.

To come up with these answers, the first thing that we have to notice is that both expressions are non-linear in $\alpha $_{g}. This means that we cannot express ‘$\alpha $_{g}’ *explicitly* as a function of the other variables. What do we use to solve equations that are non-linear in one variable? Doesn’t this ring a bell? We apply iterative techniques. And the classical iterative technique is the *Newton-Raphson procedure*. Now we can provide an answer to both of the questions that we just asked.

A distinctive characteristic of any Newton-Raphson procedure is that the success of the procedure depends greatly upon the choice of the initial guess for the variable considered. In fact, it is very commonly said that for Newton-Raphson to succeed, the initial guess must be as close as possible to the real solution. This ‘illness’ becomes worse when dealing with non-monotonic functions. In a monotonic function, derivatives at every point are always of the same sign — the function either increases or decreases monotonically. For Newton-Raphson, this means that there are neither valleys nor peaks that could lead the procedure to false solutions. If you apply Newton-Raphson to a monotonic and everywhere-continuous function, the success of the procedure is * not* dependent on the initial guess. In fact, if you apply Newton-Raphson to a monotonic function that is continuous at every single point of the domain as well, it does not matter at all where you start: you will always find

Why does this matter when dealing with equations (13.1)? The fact of the matter is that equations (13.1) are* not *monotonic, and this does not make things easier. If, as an exercise, you plot them as a function of $\alpha $_{g} or take the derivative, you will realize that both functions may change the sign of their first derivatives for different values of $\alpha $_{g}.

This poses a problem, obviously. You will not get a unique solution by applying Newton-Raphson, and you might end up with the wrong solution. Rachford and Rice (1952) recognized this problem and came up with a suggestion. They proposed a new *objective function,* based on equations (13.1), which simplifies the application of the Newton-Raphson procedure.

They combined equations (13.1) by subtraction to yield:

$$F\left({\alpha}_{g}\right)={F}_{y}\left({\alpha}_{g}\right)-{F}_{x}\left({\alpha}_{g}\right)$$ (13.2)

Hence, the new objective function becomes:

$$F\left({\alpha}_{g}\right)={\displaystyle \sum}_{i=1}^{n}\frac{{z}_{i}\left({K}_{i}-1\right)}{1+{\alpha}_{g}\left({K}_{i}-1\right)}=0$$ (13.3)

Equation (13.3) is known in the literature as the *Rachford-Rice Objective Function*. Rachford and Rice combined two very well known objective functions into a single objective function. Are there advantages to this “new” approach?

The wonderful news is that equation (13.3) is **monotonic**. The implication of this is that equation (13.3) is better suited for Newton-Raphson application than equations (13.1). How do you demonstrate the monotonic character of the Rachford and Rice objective function? To do this, we take the first derivative of the function:

$$F\text{'}\left({\alpha}_{g}\right)={\displaystyle \sum}_{i=1}^{n}\frac{{z}_{i}{\left({K}_{i}-1\right)}^{2}}{{\lfloor 1+{\alpha}_{g}\left({K}_{i}-1\right)\rfloor}^{2}}=0$$ (13.4)

Every item within the summation sign is positive — this is guaranteed by the squares in the numerator and the denominator and the fact that all compositions are positive. Hence, the derivative expression in (13.4) has no choice but to be* negative*, and the Rachford-Rice Objective function has been proven to be a *monotonically decreasing function*. With this approach, Rachford-Rice removed a major headache from the Vapor-Liquid Equilibria problem.

A remaining weakness of the Rachford-Rice objective function is that, although monotonic, it is not continuous at all points of the domain. By inspection, you can see that the function has ‘n’ singularities (as many singularities as components in the mixture), because it becomes singular at values of ‘${\alpha}_{g}$ ’ equal to:

$${\alpha}_{g}=\frac{1}{1-{K}_{i}}$$ (13.5)

Hence, you may still face convergence problems if the procedure crosses any of the singularities. If, by any means, one is able to keep the Newton-Raphson procedure within values of ${\alpha}_{g}$ where a physically meaningful solution is possible (within the two asymptotes where the values 0 < ${\alpha}_{g}$ < 1 are found), the monotonically decreasing feature of the Rachford-Rice equation would guarantee convergence.

$$F\left({\alpha}_{g}\right)={\displaystyle \sum}_{i=1}^{n}\frac{{z}_{i}\left({K}_{i}-1\right)}{1+{\alpha}_{g}\left({K}_{i}-1\right)}=0$$ (13.3)

Equation (13.3) is a non-linear equation in one variable, and the Newton Raphson procedure is usually implemented to solve it. In general, Newton Raphson is an iterative procedure with a fast rate of convergence. The method calculates a new estimate, ${\alpha}_{g}{}^{new}$ , which is closer to the real answer than the previous guess, ${\alpha}_{g}{}^{old}$ , as follows:

$${\alpha}_{g}^{new}={\alpha}_{g}^{old}-\frac{F\left({\alpha}_{g}^{old}\right)}{F\text{'}\left({\alpha}_{g}^{old}\right)}$$ (13.6)

Substituting (13.3) and (13.4) into (13.6),

$${\alpha}_{g}^{new}={\alpha}_{g}^{old}+\frac{{\displaystyle \sum}_{i=1}^{n}\frac{{z}_{i}\left({K}_{i}-1\right)}{1+{\alpha}_{g}^{old}\left({K}_{i}-1\right)}}{{\displaystyle \sum}_{i=1}^{n}\frac{{z}_{i}{\left({K}_{i}-1\right)}^{2}}{{\left[1+{\alpha}_{g}^{old}\left({K}_{i}-1\right)\right]}^{2}}}$$ (13.7)

In this iterative scheme, convergence is achieved when

$$\left|{\alpha}_{g}^{new}-{\alpha}_{g}^{old}\right|<\epsilon $$ (13.8)

where $\epsilon $ is a small number ($\epsilon \text{}=\text{}1.0\text{}x\text{}{10}^{\u2013\text{}9}$ is usually adequate). After solving for ${\alpha}_{g}$ , the liquid molar fraction and composition of each of the phases can be calculated as follows:

$$\text{LiquidMolarFraction:}{\alpha}_{l}=1-{\alpha}_{g}$$ (13.9a)

$$\text{PercentageofLiquid:}\%L=100*{\alpha}_{l}$$ (13.9b)

$$\text{PercentageofVapor:}\%V=100*{\alpha}_{g}$$ (13.9c)

$$\text{VaporPhaseComposition:}{y}_{i}=\frac{{z}_{i}{K}_{i}}{1+{\alpha}_{g}\left({K}_{i}-1\right)}$$ (12.7)

$$\text{LiquidPhaseComposition:}{x}_{i}=\frac{{z}_{i}}{1+{\alpha}_{g}\left({K}_{i}-1\right)}$$ (12.11)

In the previous development, we made one crucial assumption. We assumed that, somehow, we knew all the equilibrium ratios. The fact is, however, that we usually don’t. If we do not know all equilibrium ratios, then all of the previous discussion is meaningless. So far, the only conclusion we can draw is that if we *happen* to know K_{i}’s, the VLE problem is *solvable.*

The K_{i} value of each component in a real hydrocarbon mixture is a function of the pressure, temperature, and also of the composition of each of the phases. Since the compositions of the phases are not known beforehand, equilibrium constants are not known, either. If they were known, the VLE calculation would be performed in a straightforward manner. This is because once the equilibrium constants of each component of the mixture are known for the given pressure and temperature of the system, both gas and liquid molar fractions, ${\alpha}_{g}$ and al can be calculated by solving the Rachford-Rice Objective Function.

Nevertheless, the good news is that *sometimes* K_{i}’s are fairly independent of the phase’s composition. This is true at pressure and temperature conditions away from the critical point of the mixture. Therefore, numerous correlations have been developed throughout the years to estimate the values of K_{i} for each hydrocarbon component as a function of the pressure and temperature of the system.

To illustrate how K_{i} may be calculated as a function of pressure and temperature, let us take the case of an ideal mixture. For a mixture to behave ideally, it must be far removed from critical conditions. The fact of the matter is, in an ideal mixture, the partial pressure of a component in the vapor phase (p_{i}) is proportional to the vapor pressure (P_{sat}) of that component when in its pure form, at the given temperature. The constant of proportionality is the molar fraction of that component in the liquid (x_{i}). Then, we have:

$${p}_{i}={x}_{i}{P}^{sat}$$ (13.10)

Equation (13.10) is known as *Raoult’s law*. Additionally, if the vapor phase behaves ideally, Dalton’s law of partial pressures applies. Dalton’s law of partial pressures says that the total pressure in a vapor mixture is equal to the sum of the individual contributions (partial pressures) of each component. The partial pressure of each component is a function of the composition of that component in the vapor phase:

$${p}_{i}={y}_{i}P$$ (13.11)

Equalizing equation (13.10) and (13.11),

$${y}_{i}P={x}_{i}{P}^{sat}$$ (13.12)

We rearrange equation (13.12) to show:

$$\frac{{y}_{i}}{{x}_{i}}=\frac{{P}^{sat}}{P}$$ (13.13)

If we recall the definition of equilibrium ratios, $${K}_{i}={y}_{i}/{x}_{i}$$ , we readily see:

$${K}_{i}=\frac{{P}^{sat}}{P}$$ (13.14)

Since the vapor pressure of a pure substance (P_{sat}) is a function of temperature, we have just shown with equation (13.14) that the equilibrium ratios K_{i} are functions of pressure and temperature — and not of composition — when we are dealing with ideal substances. Vapor pressure can be calculated by a correlation, such as that of Lee and Kesler.

The estimation of equilibrium ratios (K_{i}) has been a very intensely researched subject in vapor-liquid equilibria. A number of methods have been proposed in the literature. In the early years, the most common way of estimating equilibrium ratios was with the aid of charts and graphs that provided K_{i} values as a function of pressure and temperature for various components. Charts provided better estimations of K_{i}’s than what came from the direct application of Raoult-Dalton’s derivation (equation 13.14). However, the compositional dependency could not be fully captured by the use of charts. Charts remained popular until the advent of computers, when the application of more rigorous thermodynamic models became possible.

Most of the K_{i}-charts available have been represented by empirical mathematical correlations to make them amenable for computer calculations. You may find a large number of correlations in the literature that would allow you to estimate K_{i}’s for a range of conditions. A very popular empirical correlation that is very often used in the petroleum and natural gas industry is *Wilson**’s empirical correlation.* This correlation gives the value of K_{i} as a function of reduced conditions (P_{ri}, T_{ri}: reduced pressure and temperature respectively) and Pitzer’s acentric factor and is written as:

$${K}_{i}=\frac{1}{{P}_{ri}}EXP\left[5.37\left(1+{\omega}_{i}\right)\left(1-\frac{1}{{T}_{ri}}\right)\right]$$ (13.15)

Wilson’s correlation is based on Raoult-Dalton’s derivation (equation 13.14). Therefore, it does not provide any compositional dependency for K_{i}’s and, as such, it is only applicable at low pressures (away from the critical conditions). We have included this correlation not because of its accuracy, but because it will become the initial guess that is needed to start the K_{i}-prediction procedure, which we will develop later.

In the following modules, we will develop a rigorous thermodynamic model for predicting equilibrium ratios. As you may anticipate, we will be taking advantage of what we have learned about equations of state in order to build up a phase behavior predictor model. This approach is known as the equation of state approach, or the “fugacity” approach. Even though the “fugacity” approach is the one that we will be covering in detail, you may at some point encounter the fact that equilibrium ratios can also be estimated by using the solution theory or “activity” approach. We will not be discussing the latter, as the former represents the most popular K_{i}-prediction method in the petroleum and natural gas business.

In order to guarantee a complete understanding of the thermodynamic model which we are about to apply, we need to review the basic concepts of classical thermodynamics. This is the goal for the next couple of modules (Modules 14 through 16). After all of the thermodynamic tools have been given to the reader, we will resume our discussion of Vapor-Liquid Equilibria (Module 17), along with discussing how equations of state come into the picture. Our final goal will be to tie the equilibrium ratio, K_{i}, to the thermodynamic concepts of chemical potential, fugacity, and equilibrium.

Answer the following problems, and submit your answers to the drop box in Canvas that has been created for this module.

*Please note:*

- Your answers must be submitted in the form of a Microsoft Word document.
- Include your Penn State Access Account user ID in the name of your file (for example, "module2_abc123.doc").
- The due date for this assignment will be sent to the class by e-mail in Canvas.
- Your grade for the assignment will appear in the drop box approximately one week after the due date.
- You can access the drop box for this module in Canvas by clicking on the Lessons tab, and then locating the drop box on the list that appears.

- In the older literature, “K
_{i}” were referred to as the “equilibrium constants.” Now, they are called the “equilibrium ratios.” Why should that make a difference? - In the analysis of the Rachford-Rice Objective Function, we assumed that there was only one unknown: the gas molar fraction (α
_{g}). It this strictly true? - Warren and Adewumi (1992) derived the following analytical expression for flash calculations for a binary mixture:

$${\alpha}_{g}=-\left[\frac{{z}_{1}}{{K}_{2}-1}+\frac{{z}_{2}}{{K}_{1}-1}\right]$$ Please demonstrate that this expression is true.