# Cubic EOS and Their Behavior (I)

### Module Goals

Module Goal: To demonstrate thermodynamic quantification using modern cubic EOS.

Module Objective: To introduce you to the basis premises of cubic EOS and their behavior.

# Introduction

If we multiply the vdW EOS (expression 7.11a in Module 7) by ${\stackrel{˜}{v}}^{2}$ and expand the factorized product by applying the distributive law, the result is the vdW EOS expressed in terms of molar volume, as follows:

${\stackrel{˜}{v}}^{3}-\left(b+\frac{RT}{P}\right){\stackrel{˜}{v}}^{2}+\left(\frac{a}{P}\right)\stackrel{˜}{v}-\frac{ab}{P}=0$ (9.1)

Note that equation (9.1) is a third order polynomial in $\stackrel{˜}{v}$ i.e., it is cubic in molar volume. Additionally, we can substitute the definition of compressibility factor Z,

$Z=\frac{P\stackrel{˜}{v}}{RT}$ (9.2)

into equation (9.1) and obtain a different cubic polynomial in Z, as shown:

${Z}^{3}-\left(1+\frac{bP}{RT}\right){Z}^{2}+\left(\frac{aP}{{R}^{2}{T}^{2}}\right)Z-\frac{ab{P}^{2}}{{\left(RT\right)}^{3}}=0$ (9.3)

As we see, vdW EOS is referred to as cubic because it is a polynomial of order 3 in molar volume (and hence in compressibility factor Z). In general, any equation of state that is cubic in volume (and Z) and explicit in pressure (equation 7.11b) is regarded as a cubic equation of state. vdW EOS is a cubic EOS, and all the transformations and modifications that it has undergone during the more than one hundred years since its publication are also cubic EOS; or better, they are in-the-van-der-Waals-spirit EOS or of-the-van-der-Waals-family EOS.

# The Cubic Behavior

Let us see what we have accomplished thus far. First, remember that we are not satisfied with the fact that the ideal gas law is not able to predict the discontinuity of Figure 6.3, which corresponds to the condensation that an all-gas pure substance undergoes during isothermal compression.The isotherm that we get from the ideal model was shown in Figure 6.2. Condensation, as we know, is to be expected at some point when you isothermally compress any pure gas while below its critical conditions (T < Tc). This lack of conformance to actual behavior led us to realize that the ideal model is quantitatively and qualitatively wrong at high pressures. To remove the weaknesses of the ideal model, we recall that vdW, in developing his equation of state (equation 7.11), introduced the concepts of co-volume and the attraction term.

Now, we wonder, what did he really accomplish? Are we now able to predict the condensation phenomena? Are these “new” cubic-type EOS capable of showing where such a discontinuity occurs? So far, we have not seen what a cubic isotherm looks like. Let us plot the cubic isotherm for conditions below critical (T < Tc), superimpose it on Figure 6.3, and see what we get.

Figure 9.1: P-v Isotherm for a cubic EOS

Figure 9.1 shows a typical cubic behavior. That is, around condensation conditions, the equation of state presents an S-shaped curve. This should not come as a surprise; because the equation is cubic in volume, it will provide three roots for volume, and hence the S-shaped cubic feature. For the sake of our discussion, we have shown the gas and liquid branches of the cubic equation lying directly on top of the experimental ones. This kind of matching is the “ultimate goal” of any cubic EOS, but in reality, the matching is not this good — especially for the case of the original vdW EOS, as we will discuss later.

The cubic behavior is bounded by the two extremes of real fluid behavior, given by the zero pressure and infinite pressure limits. On one hand, it is clear from equation (7.11b) that we have a singularity at $\stackrel{˜}{v}=b$ , where “b” represents the co-volume, or physical space that molecules themselves occupy. This singularity conveniently creates an asymptotic behavior (high pressure asymptote) of the cubic equation liquid branch, by which $\stackrel{˜}{v}\to b$ as $P\to \infty$ . Recall from our previous discussions that predictions for $\stackrel{˜}{v} are meaningless. Accordingly, we need an infinite amount of pressure if we are to compress a fluid to the extent that no free space is available among molecules $\left(v˜=b\right)$ . On the other hand, it is clear from vdW EOS (equation 7.11a) that as $v\to \infty \left(P\to 0\right)$ , the cubic EOS collapses to the ideal EOS (equation 7.6). At this low pressure limit, vdW
corrections to the ideal model become inconsequential .

Mathematically, this is the low-pressure asymptote of the cubic-equation gas branch by which $\stackrel{˜}{v}\to \infty$ as $P\to 0$ .

Let us see how “good” the cubic behavior is for the ideal-gas model. Refer to Figure 9.2, where we have superimposed the ideal gas isotherm of Figure 6.2 on the cubic EOS behavior.

Figure 9.2: P-v Isotherms for the ideal gas EOS and cubic EOS

A look at Figure 9.2 helps us confirm that the cubic equation of state collapses to the prediction of the ideal gas model at reduced pressure — i.e., they share the same low-pressure asymptote. This is to be expected, since the assumptions underlying the ideal model are satisfied. As we recall, these assumptions are that the attractive forces between molecules are very weak and that the physical volume of the molecules can be disregarded when compared to the total volume of the container. It is worthy to note that the high pressure asymptote is not the same for both models; as $P\to \infty$ , $\stackrel{˜}{v}\to b$ as for the cubic model while $\stackrel{˜}{v}\to 0$ for the ideal model. The latter is a direct consequence of neglecting the molecular volume in the ideal model.

# Implications of S-shaped curve (Sub-critical Conditions)

A major issue that has kept us banging our heads has been how to mathematically represent the discontinuity of the P-v isotherm during the vapor-liquid transition (Figure 6.3). Such a discontinuity shows up during the isothermal compression of any pure substance at sub-critical conditions (T < Tc). What we want here requires fitting a continuous mathematical function to a discontinuous, real-life event. Strictly speaking, it would be contradictory to find a single continuous mathematical function that can capture such a discontinuity in its full nature.

Can we really model the discontinuity? Not really, but we can get around it. van der Waals provided a possible solution in his dissertation on the “continuity of vapor and liquid.” Even though neither cubic equations nor any other continuous mathematical function is able to follow the discontinuity, what they can do is good enough for engineering purposes. The “cubic behavior” can reasonably match the liquid and vapor branches for the real, experimental isotherms.

Since van der Waals’ EOS, we have been able to consider the continuity between gas and liquid phases. Now, we need to learn how to deal with the S-shaped behavior, and to look at it as a minor, inconsequential price that we pay for the modeling of the vapor — liquid discontinuous transition with a continuous mathematical function. Let us zoom in on Figure 9.1, as shown in Figure 9.3.

Figure 9.3: S-shaped Feature of Cubic Equations

There are several features of the S-shaped behavior that should be noted.

1. The S-shaped transition represents the zone where gas and liquid coexist in equilibrium for a pure substance; hence, such a behavior will show up whenever a cubic equation of state is used for predictions at temperatures below critical conditions (sub-critical conditions).
2. Physically, changes in pressure and changes in volume in a fluid must have different signs in any isothermal process, such that:
This requirement is met by the liquid branch, gas branches and sections AA’ and BB’ of the cubic isotherm. Portions AA’ and BB’ have even been realized experimentally for metastable conditions, i.e., conditions of fragile or weak stability. However, changes of pressure and volume have the same sign in portion A’B’; consequently, this portion of the cubic isotherm is regarded as meaningless and unphysical.
3. Point A’ and B’ are the minimum and maximum values of the metastable behavior represented by the S-shaped curve. The saturation pressure (Psat) will naturally lie between these two extremes. Graphically, Psat can be attained as the pressure that make the areas AoA’ and BoB’ equal. This equal-area rule is known as the Maxwell principle. Similarly, we can also determine the saturation condition as the pressure where fugacity — a thermodynamic property we will be studying later — is equal both at the liquid and vapor branches.
4. It is not impossible for section AA’ of the cubic isotherm to reach negative pressures (i.e., PA’ < 0). This should not be of any concern to us because we are seldom interested in such metastable behavior. In fact, once Psat is determined, most practical applications call for cleaning up the cubic isotherm and suppressing areas AoA’ and BoB’. In this case, we are left with liquid branch, the discontinuity AB’ at Psat, and the gas branch; just as the experimental isotherm would look.
5. The most crucial implication of the S-shaped curve is that the cubic equation will certainly produce three distinct real roots for molar volume (or compressibility factor, if it is the case.) This will always be the case as long as you are making a prediction within the S-shaped curve (PA’ < P < PB’ ; T < Tc). At the vapor-liquid discontinuous transition — i.e., at the intercept of the cubic isotherm with the saturation pressure, Psat, at the given temperature — we end up with 3 mathematically possible real roots for volume. The extreme points of intersection represent liquid molar volume and gas molar volume respectively, represented as ${\stackrel{˜}{v}}_{liquid}$ and ${\stackrel{˜}{v}}_{gas}$ in Fig. 9.3. The third, middle root, given by point “o”, is always regarded as unphysical and is always discarded, because it belongs to the path A’B’.

# Action Item

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

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

#### Problem Set

1. Cubic EOS yield three roots for Z-factor for a pure substance at subcritical conditions. Does that mean that there is a “Z” factor for liquids? Isn’t the “Z” factor concept only applicable to gases? What does a “Z factor” connote for liquids?

# Cubic EOS and Their Behavior (II)

### Module Goals

Module Goal: To demonstrate thermodynamic quantification using modern cubic EOS.

Module Objective: To highlight the most often used cubic EOS.

# Multiple Roots and Cubic Behavior

It comes as no surprise that cubic equations of state yield three different roots for volume and compressibility factor. This is simply because they are algebraic equations, and any nth order algebraic equation will always yield “n” roots. However, those “n” roots are not required to be distinct, and that is not all: they are not required be real numbers, either. A quadratic expression (n = 2) may have zero real roots (e.g., x2 + 1 = 0); this is because those roots are complex numbers. In the case of cubic expressions (n = 3), we will either have one or three real roots; this is because complex roots always show up in pairs (i.e., once you have a complex root, its conjugate must also be a solution.) In our case, and because we are dealing with physical quantities (densities, volumes, compressibility factors), only real roots are of interest. More specifically, we look for real, positive roots such that $\stackrel{˜}{v}>b$ in the case of molar volume and $Z>Pb/RT$ in the case of compressibility factor.

In a cubic equation of state, the possibility of three real roots is restricted to the case of sub-critical conditions (T < Tc), because the S-shaped behavior, which represents the vapor-liquid transition, takes place only at temperatures below critical. This restriction is mathematically imposed by the criticality conditions. Anywhere else, beyond the S-shaped curve, we will only get one real root of the type $\stackrel{˜}{v}>b$ . Figure 10.1 illustrates this point.

Figure 10.1: Multiple Roots in Cubic EOS

Let us examine the three cases presented in Figure 10.1:

1. Supercritical isotherms (T > Tc): At temperatures beyond critical, the cubic equation will have only one real root (the other two are imaginary complex conjugates). In this case, there is no ambiguity in the assignment of the volume root since we have single-phase conditions. The occurrence of a unique real root remains valid at any pressure: any horizontal (isobaric) line cuts the supercritical isotherm just once in Figure 10.1.
2. Critical isotherm (T = Tc): At the critical point (P = Pc), vapor and liquid properties are the same. Consequently, the cubic equation predicts three real and equal roots at this special and particular point. However, for any other pressure along the critical isotherm (P < Pc or P > Pc,) the cubic equation gives a unique real root with two complex conjugates.
3. Subcritical isotherm (T < Tc): Predictions for pressures within the pressure range for metastability (PA’  < P < PB’) or for the saturation condition (P = Psat) will always yield three real, different roots. In fact, this is the only region in Figure 10.1 where an isobar cuts the same isotherm more than once. The smallest root is taken as the specific volume of the liquid phase; the largest is the specific volume of the vapor phase; the intermediate root is not computed as it is physically meaningless. However, do not get carried away. Subcritical conditions will not always yield three real roots of the type $\stackrel{˜}{v}>b$ . If the pressure is higher than the maximum of the S-shaped curve, PB, we will only have one (liquid) real root that satisfies $\stackrel{˜}{v}>b$ . By the same token, pressures between 0 < P < PA’ yield only one (vapor) root. In the case of PA’ being a negative number, three real roots are to be found even for very low pressures when the ideal gas law applies. This can be seen in Figure 10.1 as well. The largest root is always the correct choice for the gas phase molar volume of pure components.

Most of these considerations apply to the cubic equation of state in Z (compressibility factor). The most common graphical representation of compressibility factor is the well-known chart of Standing and Katz, where Z is plotted against pressure. Standing and Katz presented their chart for the compressibility factor (Z) of sweet natural gases in 1942. This chart was based on experimental data. Graphical determination of properties was widespread until the advent of computers, and thus the Standing and Katz Z-chart became very popular in the natural gas industry. Typical Standing and Katz charts are given for high temperature conditions (T > Tc or Tr > 1). Figure 10.2, using a cubic equation of state for a pure gas, presents the qualitative behavior of the solution of Z versus pressure. Isotherms (T > Tc) show the typical qualitative behavior we are accustomed to seeing in the Standing and Katz chart. Cases T < Tc (Tr < 1) are not as familiar to us, as they were not considered by Standing and Katz. For such isotherms, it is clear that you come up with two values of Z (liquid, gas) at saturation conditions.

Figure 10.2: Compressibility Factor versus Pressure

# Modern Cubic EOS

At this point, we have a couple of comments on the cubic behavior that the pioneer work of vdW introduced to the field of equations of state. First, we can say that the vdW cubic behavior is qualitatively reasonable; and second, we can say that it is capable of describing the continuity between liquid and vapor. Nevertheless, vdW cubic EOS has been proven not to be quantitatively suitable for most engineering purposes. Certainly, it yields unacceptable errors for the quantitative prediction of densities and any other related thermodynamic property. However, all of the development in the field of phase behavior that has been achieved today is due to the work of van der Waals. Although his own equation is seldom used because of its lack of accuracy, his principles are still the foundations of the current developments. vdW concepts were so far reaching that he won the Nobel Prize for his equation.

The truth is that van der Waals’ accomplishment in 1873 triggered a tremendous effort among scientists to make modifications to his EOS which would remove from it large disagreements with experimental data. This effort has not yet ceased today and is not likely to stop in the near future. Much of this endeavor has focused on how to better model the attractive parameter “a” and the repulsive term “b”, with the hope that we can get better quantitative predictions. Naturally, the qualitative cubic-nature of vdW’s original EOS is always preserved, and hence all subsequent refinements belong to the same family of modified-van-der-Waals equations of state. We refer to vdW EOS and all its descendents as cubic equations of state, because, as we have said, they take a cubic form when expressed in terms of volume or compressibility factor and are explicit in pressure.

It is fair to claim that modern cubic EOS started to make a difference when a temperature dependency was introduced to the attractive parameter “a”. Interestingly enough, van der Waals was convinced that the parameters “a” (and even “b”) of his equation of state were not necessarily constants and suggested that, indeed, some dependency on temperature could be found. A very interesting discussion on this, from van der Waals himself, is found in the lecture speech that he offered during his acceptance of the Nobel Prize in Physics, in 1910, for his work on the continuity of vapor and liquid. This speech and the biography of this great physicist, Johannes Diderik van der Waals (1837-1923), can be found in the web resources of the Nobel Prize organization.

The most popular cubic EOS, which time has proven to be most reliable, are:

• Redlich-Kwong EOS,
• Soave-Redlich-Kwong EOS (very popular among chemical engineers),
• Peng-Robinson EOS (very popular among petroleum and natural gas engineers).

Keep in mind that, once you have an EOS, you can derive virtually any property of the fluid.

# Redlich-Kwong EOS (1949)

vdW cubic equation of state had to wait almost 100 years before a real, successful improvement was introduced to it. As we stated before, this progress occurred once researchers committed themselves to finding the empirical temperature dependency of the attraction parameter “a” proposed by van der Waals. In contrast, very little attention has been paid to modifying the parameter “b” for co-volume. It makes a lot of sense that “b” would not be modified by temperature, because it represents the volume of the molecules, which should not be affected by their kinetic energy (measured in terms of temperature).

The very first noteworthy successful modification to the attraction parameter came with the publication of the equation of state of Redlich-Kwong in 1949. Redlich and Kwong revised van der Waals EOS and proposed the following expression:

$\left(P+\frac{a}{{T}^{0.5}v\left(v+b\right)}\right)\left(v-b\right)=RT$ (10.1)

Notice that the fundamental change they introduced was to the functional form of (equation 7.8, Module 7). Additionally, they introduced the co-volume “b” into the denominator of this functional form.

The important concept here is that the attraction parameter “a” of van der Waals needed to be made a function of temperature before any cubic EOS was able to do a better job of quantitatively matching experimental data. This was a realization that vdW himself had suggested, but no actual functional dependency had been introduced until the Redlich-Kwong EOS.

We know what follows at this point. To come up with an expression for “a” and “b” of equation (10.1), we apply the criticality conditions to this EOS. As we recall, imposing the criticality conditions allows us to relate the coefficients “a” and “b” to the critical properties (Pc, Tc) of the substance. Once we have done that, we obtain the definition of “a” and “b” for the Redlich-Kwong EOS,

$a=0.427480\frac{{R}^{2}{T}_{c}^{2.5}}{{P}_{c}}$ (10.2a)

$b=0.086640\frac{R{T}_{c}}{{P}_{c}}$ (10.2b)

This EOS radically improved, in a quantitative sense, the predictions of vdW EOS. We now recall that vdW-type equations are cubic because they are cubic polynomials in molar volume and compressibility factor. It comes as no surprise then, that we can transform equation (10.1) into:

${\stackrel{˜}{v}}^{3}-\left(\frac{RT}{P}\right){\stackrel{˜}{v}}^{2}+\frac{1}{P}\left(\frac{a}{{T}^{0.5}}-bRT-p{b}^{2}\right)\stackrel{˜}{v}-\frac{ab}{P{T}^{0.5}}=0$ (10.3)

and, by defining the following parameters,

$A=\frac{aP}{{R}^{2}{T}^{2.5}}$ (10.3a)

$B=\frac{bP}{RT}$ (10.3b)

and introducing the compressibility factor definition $\left(Z=\frac{P\stackrel{˜}{v}}{RT}\right)$ , we get:

${Z}^{3}-{Z}^{2}+\left(A-B-{B}^{2}\right)Z-AB=0$ (10.4)

We may also verify the two-parameter corresponding state theory by introducing equations (10.3) and (10.2) into (10.4),

${Z}^{3}-{Z}^{2}+\frac{{P}_{r}}{{T}_{r}}\left(\frac{0.42748}{{T}_{r}^{1.5}}-0.08664-0.007506\frac{{P}_{r}}{{T}_{r}}\right)Z-0.03704\frac{{P}_{r}^{2}}{{T}_{r}^{3.5}}=0$ (10.5)

In equation (10.5) we can observe the same thing that we saw with vdW EOS: gases at corresponding states have the same properties. Equation (10.5) is particularly clear about it: any two different gases at the same Pr, Tr condition have the same compressibility factor.

Just as any other cubic equation of state, equations (10.1) through (10.5), as they stand, are to be applied to pure substances. For mixtures, however, we apply the same equation, but we impose certain mixing rules to obtain “a” and “b”, which are functions of the properties of the pure components. Strictly speaking, we create a new “pseudo” pure substance that has the average properties of the mixture. Redlich-Kwong preserved the same mixing rules that vdW proposed for his EOS:

(10.6a)

${b}_{m}=\sum _{i}{y}_{i}{b}_{i}$ (10.6b)

Naturally, Redlich and Kwong did not have the last word on possible improvements to the vdW EOS. The Redlich-Kwong EOS, as shown here, is no longer used in practical applications. Research continued and brought with it new attempts to improve the RK EOS. After more than two decades, a modified RK EOS with very good potential was developed. The Soave-RK EOS was born.

# Soave-Redlich-Kwong EOS (1972)

In 1972, Soave proposed an important modification to the RK EOS — or shall we say, a modification to vdW EOS. Between the time of vdW EOS and Redlich-Kwong’s, a new concept for fluid characterization was being discussed. Pitzer had introduced the concept of acentric factor in 1955.

All modifications to the vdW EOS had focused on the temperature dependency of the attractive parameter. Soave expanded this by proposing a two-variable dependency for “a”:

$a=a\left(T,\omega \right)$ (10.7)

It was the first time that “a” was expressed not only as a function of temperature, but also as a function of the shape (sphericity) of the molecules (through w, Pitzer’s acentric factor). As we recall, Pitzer’s acentric factor is a measure of the configuration and sphericity of the molecule. It can also be seen as a measure of the deformity of the molecule.

The Soave-Redlich-Kwong EOS is given by the expression:

$\left(P+\frac{\alpha a}{\stackrel{˜}{v}\left(\stackrel{˜}{v}+b\right)}\right)\left(\stackrel{˜}{v}-b\right)=RT$ (10.8a)

Like all cubic equations of state, the SRK EOS is also explicit in pressure. Notice, for example, how the SRK EOS readily becomes:

$P=\frac{RT}{\stackrel{˜}{v}-b}-\frac{\alpha a}{\stackrel{˜}{v}\left(\stackrel{˜}{v}+b\right)}$ (10.8b)

where,

$\alpha ={⌊1+\left(0.48508+1.55171\omega -0.15613{\omega }^{2}\right)\left(1-\sqrt{Tr}\right)⌋}^{2}$ (10.8c)

The influence of acentric factor and temperature on the attractive term is introduced now through “a”. What do we do next? We apply the criticality conditions to equation (10.8b). Notice that expression (10.8c) becomes unity at Tr=1, throughout the critical isotherm. We obtain:

$a=0.427480\frac{{R}^{2}{T}_{c}^{2}}{{P}_{c}}$ (10.9a)

$b=0.086640\frac{R{T}_{c}}{{P}_{c}}$ (10.9b)

Now we show the cubic form (in compressibility factor) of Soave-Redlich-Kwong EOS. Defining,

$A=\frac{\alpha aP}{{R}^{2}{T}^{2}}$ (10.10a)

$B=\frac{bP}{RT}$ (10.10b)

we are able to obtain:

${Z}^{3}-{Z}^{2}+\left(A-B-{B}^{2}\right)Z-AB=0$ (10.11)

For mixtures, Soave proposed a “little” modification to the mixing rules with which we have dealt with so far by introducing the use of “binary interaction parameters” (kij):

(10.12a)

${b}_{m}=\sum _{i}{y}_{i}{b}_{i}$ (10.12b)

The use of binary interaction parameters (kij) generated a lot of resistance upon their first introduction. This is because there is no analytical, science-based derivation that justifies their existence. Nowadays, they are regarded just as they are, empirical factors used to tune equations of state and make them match experimental data for mixtures. This has become the heuristic justification for their existence: with them, EOS can do a better job of matching experimental data. Heuristically speaking, they are a measure of interaction between a pair of dislike molecules. Based on this “definition,” their value is zero for pairs of molecules that are alike. Actually, this is no more than a mathematical requirement in order for equation (10.12a) to give ${\left(\alpha a\right)}_{ij}={\left(\alpha a\right)}_{i}$ when j=i. The determination of kij is based on experimental data from binary systems; “kij” results from the value that allows the given equation of state (through the expression in 10.12a) to yield the closest match. These values are assumed to be constant (and so are used) when the same two components are part of a more complex multi-component mixture.

# Action Item

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

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

#### Problem Set

1. An engineer is seeking your advice on which EOS he should use to model various hydrocarbon mixtures he is dealing with. He has access to a program that provides him with several options for different cubic EOS. He presents you with the following list:
• A sweet natural gas from Nigeria
• A sour natural gas (high CO2, H2S contents)
• Flue gas
Which EOS would you recommend him to use for each of these cases?

# Cubic EOS and Their Behavior (III)

### Module Goals

Module Goal: To demonstrate thermodynamic quantification using modern cubic EOS.

Module Objective: To assess the relative merit of applying the most common EOS.

# Peng-Robinson EOS (1976)

The Peng-Robinson EOS has become the most popular equation of state for natural gas systems in the petroleum industry. During the decade of the 1970’s, D. Peng was a PhD student of Prof. D.B. Robinson at the University of University of Alberta (Edmonton, Canada). The Canadian Energy Board sponsored them to develop an EOS specifically focused on natural gas systems. When you compare the performance of the PR EOS and the SRK EOS, they are pretty close to a tie; they are “neck to neck,” except for a slightly better behavior by the PR EOS at the critical point. A slightly better performance around critical conditions makes the PR EOS somewhat better suited to gas/condensate systems.

Peng and Robinson introduced the following modified vdW EOS:

$\left(P+\frac{\alpha a}{{\stackrel{˜}{v}}^{2}+2b\stackrel{˜}{v}-{b}^{2}}\right)\left(\stackrel{˜}{v}-b\right)=RT$
(11.1a)

or, explicitly in pressure,

$P=\frac{RT}{\stackrel{˜}{v}-b}-\frac{\alpha a}{{\stackrel{˜}{v}}^{2}+2b\stackrel{˜}{v}-{b}^{2}}$
(11.1b)

where:

Peng and Robinson conserved the temperature dependency of the attractive term and the acentric factor introduced by Soave. However, they presented different fitting parameters to describe this dependency (equation 4.11c), and further manipulated the denominator of the pressure correction (attractive) term. As we have seen before, coefficients “a” and “b” are made functions of the critical properties by imposing the criticality conditions. This yields:

$\alpha ={⌊1+\left(0.37464+1.54226\omega -0.26992{\omega }^{2}\right)\left(1-\sqrt{{T}_{r}}\right)⌋}^{2}$
(11.1c)

$a=0.45724\frac{{R}^{2}{T}_{c}^{2}}{{P}_{c}}$
(11.2a)

The PR cubic expression in Z becomes:

$b=0.07780\frac{R{T}_{c}}{{P}_{c}}$
(11.2b)

where:

${Z}^{3}-\left(1-B\right){Z}^{2}+\left(A-2B-3{B}^{2}\right)Z-\left(AB-{B}^{2}-{B}^{3}\right)=0$
(11.3a)

$A=\frac{\alpha aP}{{R}^{2}{T}^{2}}$
(11.3b)

$B=\frac{bP}{RT}$
(11.3c)

Similar to SRK, the PR mixing rules are:

(11.4a)

${b}_{m}=\sum _{i}{y}_{i}{b}_{i}$
(11.4b)

where binary interaction parameters (kij) again play the important empirical role of helping to better fit experimental data. Due to the empirical character of these interaction parameters, kij’s calculated for PR EOS are unlikely to be the same as the kij’s calculated for SRK EOS for the same pair of molecules.

# Comparative Assessment of RK, SRK, and PR EOS

Over the years, these EOS have been tested, and some comparisons can be made. As an engineer, you have to be able to decide which EOS best fits your purposes.

Redlich Kwong EOS:

• Generally good for gas phase properties.
• Poor for liquid phase properties.
• Better when used in conjunction with a correlation for liquid phase behavior.
• Satisfactory for gas phase fugacity calculation @ Pr < Tr/3.
• Satisfactory for enthalpy departure and entropy departure calculations.

Soave Redlich Kwong & Peng Robinson EOS

• Serve similar functions as the Redlich Kwong EOS but require more parameters.
• PR obtains better liquid densities than SRK.
• Overall, PR does a better job (slightly) for gas and condensate systems than SRK. However, for polar systems, SRK always makes a better prediction, but in the petroleum engineering business we do not usually deal with those.

# Zc as a Measure of Goodness of an EOS

### Some values for critical compressibility factors are shown below:

#### Experimental Zc for Various Substances

CO2 = 0.2744
CH4 = 0.2862
C2H6 = 0.2793
nC5 = 0.2693
nC6 = 0.2659

The values of critical compressibility factors shown here are relatively close to each other, but, in actuality, they are different. They are, in fact, substance-dependent. This is a striking finding if we recall our highly-praised Principle of Corresponding States. Didn’t we say that at the same reduced conditions, all substances “must” have, at the very least, the same compressibility factor, Z? Here is a case where we have different substances at the same corresponding states (Pr = Tr = 1, right at the critical point) but different “Z” values. After all the Corresponding States Principle is not infallible (as was stated by Pitzer). As we recall, he proposed the introduction of a third parameter (acentric factor) into the corresponding state definition to alleviate these kinds of “problems.”

At the very least, we may say that the values of Zc (compressibility factor at the critical point) of different substances are “close enough” among themselves. That is, they are not “grossly different,” so as to say that the application of the two-parameter corresponding states principle would be outrageous at the critical point. The fact of the matter is, as a consequence of the Corresponding States Principle, all cubic EOS predict a “unique” and “universal” value of Z at the critical point, regardless of the substance. The list below tells us how they perform.

#### “Universal” Zc predicted by different EOS

Ideal EOS = 1.000
vdW EOS = 0.375
RK EOS = 0.333
SRK EOS = 0.333
PR EOS = 0.301

From the list above, we would have been expecting some kind of “average” Zc of 0.27 or so. But none of the equations of state we have studied is capable of predicting a value that low. The “best” job is done by PR EOS, which provides the “closest” match to the real values observed for most substances. This illustrates why the PR EOS performs somewhat better near critical conditions.

# Advantages of Using Cubic Equations of State

All cubic equations of state have their foundation in vdW EOS. The use of cubic equations of state has become widespread because of their advantages:

• Simplicity of application
• Only a few parameters need to be determined
• Low computational overhead is required to implement them. This was a critical issue, particularly in the early stages of computers; it is not really anymore. Nevertheless, this feature is still a “plus.”

The engineer using cubic equations of state must also be aware of the disadvantages that they all share. The most important ones are the limited accuracy that they can provide, particularly for complex systems. In these cases, the procurement of empirical adjustments through the use of the binary interaction parameters (kij) is essential.

# Solution Techniques for Cubic Expressions and Root Finding

Hopefully, we have convinced you that the use of cubic equations of state can represent a very meaningful and advantageous way of modeling the PVT behavior of petroleum fluids. What we need now are the tools that will allow us to get the information that we want out of them. Even though cubic equations of state are explicit in pressure, pressure is not the common unknown to be calculated in the typical problem. In the most common problem, pressure and temperature are known and we want either molar volume (or its reciprocal, molar density) or compressibility factor (the most likely case). Therefore, we are faced very often with the need to solve for the roots of a cubic expression. Here we present a number of approaches that may be followed.

#### Analytical scheme

Given the cubic polynomial with real coefficients: x3 + ax2 + bx + c = 0, the first step is to calculate the parameters:

(11.5)

Now let M = R2 – Q3 be the discriminant. We then consider the following cases:

1. If M < 0 (R2 < Q3), the polynomial has three real roots. For this case, compute $\theta =\mathrm{arccos}\left(\frac{R}{\sqrt{{Q}^{3}}}\right)$ and calculate the three distinct real roots as:

${x}_{1}=-\left(2\sqrt{Q}\mathrm{cos}\frac{\theta }{3}\right)-\frac{a}{3}$ (11.6a)

${x}_{2}=-\left(2\sqrt{Q}\mathrm{cos}\frac{\theta +2\pi }{3}\right)-\frac{a}{3}$ (11.6b)

${x}_{3}=-\left(2\sqrt{Q}\mathrm{cos}\frac{\theta -2\pi }{3}\right)-\frac{a}{3}$ (11.6c)

Note that x1, x2, x3 are not given in any special order, and that $\theta$ has to be calculated in radians.

2. If M > 0 (R2 > Q3), the polynomial has only one real root. Compute:

$S=\sqrt[3]{-R+\sqrt{M}}$ (11.7a)

$T=\sqrt[3]{-R-\sqrt{M}}$ (11.7b)

and calculate the real root as follows:

${x}_{1}=S+T-\frac{a}{3}$ (11.7c)

Two complex roots (complex conjugates) may be found as well. However, they are of no interest for our purposes, and thus no formulas are provided. Such formulas can be found in the following suggested readings:

• W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran, 2nd Edition, Cambridge Univ. Press, (1992), p. 179.
• Spiegel, M., Liu, J., Mathematical handbook of formulas and tables, 2nd Edition, Schaum’s Outline Series, McGraw Hill, p.10.

Sometimes, the equations for S and T listed above cause problems while programming. This usually happens whenever the computer/calculator performs the cubic root of a negative quantity. If you want to avoid such a situation, you may compute S’ and T’ instead:

$S\text{'}=-sign\left(R\right)\sqrt[3]{abs\left(R\right)+\sqrt{M}}$ $T\text{'}=Q/S\text{'}$ (making T’=0 when S’=0)

where:

abs(R) = Absolute value of R

and

sign(R) = (+1) or (– 1) if R is positive or negative respectively.

It may be defined as:

sign(R) = R/Abs(R)

and then the real root is:

${x}_{1}=S\text{'}+T\text{'}-\frac{a}{3}$

Keep in mind the following useful relationships among the roots of any cubic expression:

${x}_{1}+{x}_{2}+{x}_{3}=-a$ (11.8a)

${x}_{1}{x}_{2}+{x}_{2}{x}_{3}+{x}_{3}{x}_{1}=+b$ (11.8b)

${x}_{1}{x}_{2}{x}_{3}=-c$ (11.8c)

#### Self-Check

Test your understanding of cubic root calculations by analytical means by solving the following examples.

With general cubic expressions,

$\begin{array}{l}{x}^{3}-6{x}^{2}+11x-6=0⇒\hfill \\ {x}_{1}=1\hfill \\ {x}_{2}=2\hfill \\ {x}_{3}=3\hfill \end{array}$
$\begin{array}{l}{x}^{3}+7{x}^{2}+49x+343=0\hfill \\ {x}_{1}=-7\hfill \\ {x}_{2}=0+7i\hfill \\ {x}_{3}=0-7i\hfill \end{array}$
$\begin{array}{l}{x}^{3}+2{x}^{2}+3x+4=0⇒\hfill \\ {x}_{1}=-1.65063\hfill \\ {x}_{2}=-0.174685+1.54687i\hfill \\ {x}_{3}=-0.174685-1.54687i\hfill \end{array}$

with EOS in terms of volume (v), for a pure component,

${v}^{3}-7.8693{v}^{2}+13.3771v-6.5354=0$
=> One real root (one phase),
${v}_{1}=5.7357{v}_{3}-15.6368{v}_{2}+30.315v-14.8104=0$
=> Two possible phases (three real roots),
v1 = 0.807582 (liquid phase)
v2 = 1.36174  (rejected)
v3 = 13.4675(gas phase)

with EOS in terms of the compressibility factor (z), for a pure component,

${z}^{3}-1.0595{z}^{2}+0.2215z-0.01317=0$
=> One phase,
$\begin{array}{l}z=0.8045\hfill \\ {z}^{3}-{z}^{2}+0.089z-0.0013=0\hfill \end{array}$ => Two possible phases,

#### Numerical Scheme

The Newton-Raphson method provides a useful scheme for solving for a non-explicit variable from any form of equation (not only cubic ones). Newton-Raphson is an iterative procedure with a fast convergence, although it is not always capable of providing an answer — because a first guess close enough to the actual answer must be provided.

In solving for “x” in any equation of the type f(x)=0, the method provides a new estimate (new guess) closer to the actual answer, based on the previous estimate (old guess). This is written as follows:

${x}_{new}={x}_{old}-\frac{f\left({x}_{old}\right)}{f\text{'}\left({x}_{old}\right)}$ (11.9)

Considering a cubic equation $f\left(x\right)={x}^{3}+a{x}^{2}+bx+c=0$ , the previous equation takes the form:

${x}_{new}={x}_{old}-\frac{{x}_{old}^{3}+a{x}_{old}^{2}+b{x}_{old}+c}{3{x}_{old}^{2}+2a{x}_{old}+b}$ (11.10)

The iterations continue until no significant improvement for “xnew” is achieved, i.e., | xnew – xold | < tolerance. An educated guess must be provided as the starting value for the iterations. If you are solving a cubic equation in Z (compressibility factor), it is usually recommend to take Z = bP/RT as the starting guess for the compressibility of the liquid phase and Z = 1 for the vapor root.

Semi-analytical Scheme

If you used the previous numerical approach to calculate one of the roots of the cubic expression, the semi-analytical scheme can give you the other two real roots (when they exist). By using the relationships given before, with the value ‘x1’ as the root already known, the other two roots are calculated by solving the system of equations:

${x}_{2}+{x}_{3}=-a-{x}_{1}$ (11.11a)

${x}_{2}{x}_{3}=-c/{x}_{1}$ (11.11b)

This procedure can be reduced to the following steps:

1. Let ${x}^{3}+a{x}^{2}+bx+c=0$ be the original cubic polynomial and “E” the root which is already known (x1 = E). Then, we may factorize such a cubic expression as:

$\left(x-E\right)\left({x}^{2}+Fx+G\right)=0$ (11.12a)

where:

$F=a+E$ (11.12b)

$G=-c/E$ (11.12c)

2. Solve for x2, x3 by using the quadratic expression formulae,

${x}_{1}=E$ (11.13a)

${x}_{2}=\frac{-F+\sqrt{{F}^{2}-4G}}{2}$ (11.13b)

${x}_{3}=\frac{-F-\sqrt{{F}^{2}-4G}}{2}$ (11.13c)

# Action Item

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

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

#### Problem Set

1. What is the meaning of Zc? (critical compressibility factor). Why do you think it is used as a measure of the performance of an EOS?
2. What do you think should be the value of “Z” at the critical point? Is it “one” (1)? Why? Is this to be expected?

# Elementary Vapor-Liquid Equilibrium (I)

### Module Goals

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

Module Objective: To establish the connection between real production processes and vapor liquid equilibrium calculation.

# Introduction

In the next series of modules, we are going to look at how we can apply what we have learned so far to perform vapor/liquid equilibria calculations.

1. First of all, we are going to review engineering systems and how the phenomenon of VLE is related to them,
2. Then we will proceed by looking at how to describe the problem itself,
3. We then will discuss the formulation the problem,
4. And finally, we will discuss solution strategies.

These are the four main topics that we will look at in this module. As far as VLE is concerned, we can list a number of systems that are at the heart of petroleum fluid production that involve this phenomenon:

• Separators
• Reservoir
• Pipelines
• Wellbore
• LNG Processing
• NGL Processing
• Storage
• Oil and LNG Tankers.

Vapor/liquid equilibrium pertains to all aspects of petroleum production with which we are concerned. It is no wonder, then, that we devote a new module to the subject itself.

Consider the case of a typical transmission pipeline. As gas is injected at the inlet, the pressure will drop continuously along the length of the pipeline, due to friction. Even though we usually think of liquid as forming with increasing pressure — i.e., upon compression — we have to recall that the phenomena of retrograde condensation (discussed in Module 4) takes place in hydrocarbon mixtures. Therefore, contrary to expectations, most single-phase natural gasses yield liquid upon expansion. Therefore, as pressure drops in the pipeline, liquid may drop out as the thermodynamic path crosses the dew point line and enters the phase envelope. In this case, what started as a single phase flow became two-phase flow within the system.

We also encounter this phenomenon in gas condensate reservoirs. Your initial reservoir conditions may be outside the phase envelope, but as you deplete the reservoir, your production path may take the system inside the two-phase region. In these previous two examples, the single-most important property with which you are concerned is the dew point. You want to know dew point, since you would like to know at which point of the pipeline or at which stage of production liquid may start to form.

We also may have an oil reservoir, where the initial pressure and temperature conditions place us in the single-phase liquid region. As you produce, you deplete the reservoir and enter the two-phase region by crossing the bubble point curve. At this point, we would like to know the bubble point of the system so that we may anticipate the appearance of a gas phase within an originally-all-liquid reservoir.

In all these cases, by taking a sample of the fluid to the lab, we may be able to find the composition of the fluid. Hence, in these kinds of problems, composition is usually known, and so are temperature and pressure. Your unknowns are the dew point or bubble point condition.

Suppose we are not interested in what is happening in the reservoir, but rather in what is happening at the surface. You would then like to know how much liquid or gas you will have in your separators. In this case, you are no longer interested in bubble points or dew points, but rather the extent of the phases: how much liquid and how much gas the reservoir will be able to deliver at the surface. In this case, the composition of what is coming to your separator may be known, and the pressure and temperature of operation of each separation stage may be specified. We would want to know the quality and the quantity of what is coming out; that  is, we need the compositions of the gas and oil that leave the separators and the flow rates of gas (MSCF/D) and oil (STB/D).

# Types of VLE Problems

In a typical problem of liquid and vapor coexistence, we are usually required to know one or more of the following:

• The phase boundaries,
• The extent of each phase,
• The quality of each phase.

The main emphasis is on the quantitative prediction of the above. These three represent the three basic types of VLE problems. A more detailed description of each of them is given below.

1. Phase Boundary Determination Problem

These types of problems are either a bubble-point or a dew-point calculation. They are mathematically stated as follows:

• Bubble-point T calculation: Given liquid composition (xi) and pressure (P), determine the equilibrium temperature (T),
• Bubble-point P calculation: Given liquid composition (xi) and temperature (T), determine the equilibrium pressure (P),
• Dew-point T calculation: Given vapor composition (yi) and pressure (P), determine the equilibrium temperature (T),
• Dew-point P calculation: Given vapor composition (yi) and temperature (T), determine the equilibrium pressure (P).
2. Relative Phase Quantity Determination

In this type of problem, overall composition (zi), pressure (P), and temperature (T) are given, and the extent of the phases (molar fractions of gas and liquid) are required.

3. Phase Quality Determination

In this type of problem, overall composition (zi), pressure (P), and temperature (T) are given, and the composition of the liquid and vapor phases is required.

Problems of types 2 and 3 are collectively referred to as flash calculation problems. All three are problems that we encounter in production operations as petroleum engineers. Our focus now is on solving these sorts of problems. We want to use a predictive approach to do so. This is, we want to use mathematical models — the most economical and convenient approach — to accomplish the task.

One of the assumptions that we are making here is that of equilibrium. We assume that at all times, vapor and liquid that are coexisting together are in equilibrium. Are they really in equilibrium? No! Nevertheless, the state of current knolwedge requires us to assume equilibrium so as to be able to proceed.

# Formulation of the VLE Problem

We assume that the system is at steady state and at a state of equilibrium. Adopting these two assumptions is essential to the developing of the equations that we use to solve these problems. These assumptions are convenient for modeling and have proven useful in representing the real phenomena. In conclusion, we claim that the systems maintain a state that resembles equilibrium and does not depart from it greatly. And what does steady state means? Simply stated, we say that a system is at steady state when whatever comes into the system goes out. This is, no accumulation takes place within the system.

Let us consider the equilibrium cell shown in Figure 12.1. “F” moles of a feed enter our equilibrium cell with a composition “zi” and “nc” is the number of components that we have in the mixture. A flash vaporization takes place at a given pressure and temperature, and two streams come out: “V” moles of a vapor of composition “yi” and “L” moles of a liquid of composition “xi”. In steady state, a simple overall balance yields:

Figure 12.1: The Flash Problem

$F=L+V$ (12.1a)

Now we define the fractions of gas and liquid to be, respectively:

${\alpha }_{x}=\frac{V}{F}$ (12.1b)

${\alpha }_{l}=\frac{L}{F}$ (12.1c)

Therefore, if we divide equation (12.1a) by “F”, we get:

$1={\alpha }_{l}+{\alpha }_{g}$ (12.2)

The same steady state assumption applies for the mass of each component separately. Here we revisit a concept we used in Module 5, when we studied the lever rule. At that point, we said that the number of moles of a component “i” per mole of mixture in the liquid phase is given by the product “xiaL”, while the number of moles of “i” per mole of mixture in the gas is given by “yiaG”. Since there are “zi” moles of component “i” per mole of mixture coming into the system, the conservation of each component in the system imposes:

(12.3)

Equation (12.3) is true for each of the components in the system. Equation (12.2) can be introduced into equation (12.3) to yield:

${y}_{i}{\alpha }_{g}+{x}_{i}\left(1-{\alpha }_{g}\right)={z}_{i}$ (12.4)

One of the concepts that we normally use in vapor-liquid equilibria is that of the equilibrium ratio, Ki. In fact, most of the computations of phase behavior of natural gas mixtures are carried out through the concept of the equilibrium ratio. By definition, the equilibrium ratio of a component “i” in a vapor-liquid mixture is defined as the ratio of the molar composition of that component in the vapor phase to that in the liquid phase,

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

In earlier literature, this concept was referred to as the equilibrium constant. In actuality, Ki is not constant but a function of the pressure, temperature, and composition of the system. However, equilibrium ratios can be fairly independent of composition when the pressure and temperature conditions are far from critical.

Therefore, today we refer to it as the vapor-liquid equilibrium ratio, Ki. We can introduce this concept into the balance in (12.4), as shown:

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

Now, solving for yi,

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

A constraint that mole fractions must satisfy is that they must add up to unity. Since we solved for yi, we can impose that the summation of all molar vapor fractions must be equal to one, i.e.,

$\underset{{n}_{c}}{\overset{i=1}{{\sum }^{\text{​}}}}{y}_{i}=1$ (12.8)

If we now substitute (12.7) into (12.8), we get:

$\underset{{n}_{c}}{\overset{i=1}{{\sum }^{\text{​}}}}\frac{{z}_{i}{K}_{i}}{1+{\alpha }_{g}\left({K}_{i}-1\right)}=1$ (12.9)

This equation is important for us; we call it an objective function because we can use it as the starting point for solving the vapor-liquid equilibrium problems we have posed.

However, as you may be thinking right now, this is not the only choice that we have for an objective function. In fact, we may obtain another objective function if we repeat the previous steps, while solving instead for xi. In this case, we may introduce the concept of equilibrium ratio in (12.5) into (12.6) as follows:

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

We now solve for xi,

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

If we apply the constraint that all mole fractions must add up to one,

$\sum _{i=1}^{{n}_{c}}\frac{{z}_{i}}{1+{\alpha }_{g}\left({K}_{i}-1\right)}=1$ (12.12)

Both (12.9) and (12.12) are plausible objective functions. Either of them allows us to solve the flash problem that we are dealing with. The variables that make up both equations are:

nc = Number of components,
zi = Overall composition, or composition of the feed,
Ki = Equilibrium ratios of each of the components of the mixture,
ag = Vapor fraction in the system.

What is it that we are looking for? Go back and look at the types of VLE problems that we would like to solve, as we presented them in the previous section. If we are interested in solving the flash problem, we want to know how much liquid and gas we will have inside the flash equilibrium cell. This is, given a liquid-vapor mixture of composition zi, and nc number of components, what percent of the total number of moles is liquid, and what percent is vapor? How do we split it? In this case, we would like to come up with a value for αl and αg respectively.

Equations (12.9) and (12.12) tell us that if we are able to come up with the proper values for the equilibrium ratios, Ki, which are functions of the pressure, temperature, and composition of the system, the only unknown left to solve for would be αg — exactly what we want!

Well, do not rush. We would have to come up with a way of calculating Ki’s first, and this may not be a trivial task. For the time being, let us say we “know” Ki’s. Two questions remain unanswered:

1. First, is it “better” to solve the problem using equation (12.9) or (12.12)? (Recall, either of them would lead us to the answer!).
2. Second, how do we solve for αg? For a complex mixture of many components, “αg” cannot be calculated explicitly.

We will address both of these questions in the next module. As for now, let us give you a hint: we will not use either equation (12.9) or (12.12) to solve the flash problem!

# Action Item

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

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

#### Problem Set

1. You are given two different phase envelopes (VLE region) that represent two different gas reservoir fluids. The first of them has a larger PT envelope covering larger P,T ranges. The other covers a fairly narrow P,T range. What do you conclude? Speculate on the type of gas reservoir for each case. Justify your answer.
2. In order to design the appropriate production scheme for a gas reservoir recently discovered, what type of VLE calculation do you think is more valuable? Explain in detail.

# Elementary Vapor-Liquid Equilibrium (II)

### Module Goals

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

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

# Analysis of Objective Functions

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)=\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)=\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 {zi; i = 1,2,…,n} and, if for some reason, let us say, we are able to obtain all the equilibrium ratios {Ki; 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 {yi, xi}:

${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$gexplicitly 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 the solution. It might take time, but you will be able to converge to a unique solution.

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)=\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)=\sum _{i=1}^{n}\frac{{z}_{i}{\left({K}_{i}-1\right)}^{2}}{{⌊1+{\alpha }_{g}\left({K}_{i}-1\right)⌋}^{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.

# Objective Function and Newton-Raphson Procedure

### We have seen that from a molar material balance applied to a two-phase system in equilibrium, and the definition of Ki, we can derive the Rachford and Rice Objective Function:

$F\left({\alpha }_{g}\right)=\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, , 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{\sum _{i=1}^{n}\frac{{z}_{i}\left({K}_{i}-1\right)}{1+{\alpha }_{g}^{old}\left({K}_{i}-1\right)}}{\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

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

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

(13.9a)

(13.9b)

(13.9c)

(12.7)

(12.11)

# Do We Really Know Ki?

### The Need For More Advanced Tools

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 Ki’s, the VLE problem is solvable.

The Ki 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 Ki’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 Ki for each hydrocarbon component as a function of the pressure and temperature of the system.

To illustrate how Ki 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 (pi) is proportional to the vapor pressure (Psat) 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 (xi). 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 (Psat) is a function of temperature, we have just shown with equation (13.14) that the equilibrium ratios Ki 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 (Ki) 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 Ki values as a function of pressure and temperature for various components. Charts provided better estimations of Ki’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 Ki-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 Ki’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 Ki as a function of reduced conditions (Pri, Tri: 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 Ki’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 Ki-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 Ki-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, Ki, to the thermodynamic concepts of chemical potential, fugacity, and equilibrium.

# Action Item

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

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

#### Problem Set

1. In the older literature, “Ki” were referred to as the “equilibrium constants.” Now, they are called the “equilibrium ratios.” Why should that make a difference?
2. 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?
3. 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.

# Thermodynamic Tools (I)

### Module Goals

Module Goal: To establish the mathematical framework for thermodynamics of phase equilibrium.

Module Objective: To highlight the importance of thermodynamic functions as functions of state.

# Introduction

Thermodynamics plays an important role in science and engineering. Most physical processes of engineering interest operate on thermodynamic principles. Thermodynamics is the science dealing exclusively with the principles of energy conversion. This message is conveyed by the word itself: “thermo” means “heat” — a manifestation of energy — and “dynamics” deals with “motion.”

Often times, we do not use thermodynamics as much as we should. This is either because we do not know how useful a tool it can be, or, simply because we do not know how to use it.

The basic ingredients for the utilization of thermodynamics are:

• The ability to identify and define the system that best characterizes a given process
• The availability of relevant information
• Sound engineering judgment

The goal of this module, combined with all the background that we already have built on fluid behavior, is to establish a foundation for us to cultivate a culture of the appropriate utilization of thermodynamics.

What kinds of problems can be solved with thermodynamics? There are three general classes of thermodynamic problems.

1. System property variation determination. In this case, we are given a process that takes place under certain known constraints. We are required to determine how the system properties vary.
2. Interactions that cause change. In these kinds of problems, the changes desired in the system properties are prescribed. We are required to determine the amount of external interactions that are needed to cause these changes — i.e., how to vary constraints to obtain the required final system state.
3. The best-path problem. Here we are given a system with constraints and a desired variation. We are required to determine the best way to accomplish the change — i.e., the best way to reach a goal.

In science and engineering, mathematical rigor is of essence. The next sections will take us through the mathematics of thermodynamics.

# Basic Definitions

In thermodynamics, a system is a spatial domain bounded for the purpose of describing a problem; while the surroundings are the entire spatial domain outside of the system. The communication between them is established through the of the system. The system and the surroundings make up the universe.

The system and the surroundings interact with each other. As we discussed above, one type of thermodynamic problem is that of predicting changes in a system due to interactions with its surroundings. A system is open if it can exchange mass with the surroundings, and closed if it does not exchange mass with the surroundings. A system is adiabatic if it does not exchange heat energy with the surroundings. We called a system isolated if there is neither heat nor mass crossing its boundaries.

Thermodynamic properties can be divided into two general classes: intensive and extensive properties. An intensive property is one whose value is independent of the size, extent, or mass of the system, and includes pressure, temperature, and density. By contrast, the value of an extensive property changes directly with the mass. Mass and volume are examples of extensive properties. Extensive properties per unit mass, such as specific volume, are intensive properties.

A system is homogeneous if it is has uniform properties throughout, i.e. a property such as density has the same value from point to point in a macroscopic sense. A phase is defined as a quantity of matter that is homogeneous throughout. Hence, a homogeneous system is, essentially, a single-phase system. A heterogeneous system is one with non-uniform properties, and hence, is made up of phases which can be distinguished from one another by the presence of interfaces.

Thestate is the thermodynamic coordinate of the system, specified by a number of intensive variables. The degree of freedom is the number of intensive variables needed to define the state of the system. State functions are those whose changes depend on their end states only and are independent of the path between them.

A process is the series of successive, intermediate states that the system goes through in order to go from an initial to a final state. A process is isothermal or isobaric if the temperatures or pressures of all successive steps are the same, respectively. A reversible process is one for which the exchange of energy between the system and its surroundings takes place under vanishing gradients or driving forces; that is, the properties of the system and surroundings are balanced. In a reversible process, each step of the process can be “reversed” and the original states of the system and surroundings can be restored. Any process that does not take place under infinitesimal gradients is irreversible. Strictly speaking, a reversible process is an idealization.

# Thermodynamics of Systems of Constant Composition (Closed Systems)

Thermodynamics cannot tell about the rate (kinetics) of a process, but it can tell whether or not it is possible for a process to occur. For this, we make use of the first and second laws of thermodynamics.

From our basic courses in thermodynamics, we recall that the first law of thermodynamics for a closed system is written as follows:

$\partial W=PdV$ (14.1)

where,

U = internal energy,
Q = heat added to or extracted from the system,
W = work done by or to the system.

If we want to improve the definition of the first law of thermodynamics stated in (14.1), we need to come up with an expression for the amount of work and heat. The amount of work need to accomplish a reversible process is given by the expression:

$\partial Q=TdS$ (14.2)

where:

P = pressure,
V = volume.

Additionally, we can compute the heat required to accomplish a reversible process by virtue of the second law of thermodynamics:

$\partial Q=TdS$ (14.3)

where:

T = temperature,
S = entropy.

We substitute (14.2) and (14.3) into (14.1) to get:

$dU=TdS-PdV$ (14.4)

which is the fundamental thermodynamic relationship used to compute changes in Internal Energy (U) for a closed system. This is only a restatement of the first law of thermodynamics.

Although equations (14.2) and (14.3) are applicable strictly to reversible processes, equation (14.4) is quite general and does not have such a constraint. Equation (14.4) applies to reversible and irreversible processes. Are you surprised by this? Okay, it is about time that we start to get a feeling for the implications of the fact that most thermodynamic properties are state functions. Internal energy (U) is a state function, and as such, its changes (dU) do not depend on the path that was taken (reversible/irreversible), but on the end points of process. Hence, if taking any path to compute “dU” is okay, why not take the reversible path? We do not have nice, explicit equations to describe work and heat for irreversible processes. Why should we bother? Generally speaking, “reversibility” is a thermodynamic trick that helps us to get away with many thermodynamic manipulations by taking advantage of the properties of state functions.

The formal, fundamental thermodynamic definition of Enthalpy (H) is the following:

H = U + P $\stackrel{˜}{v}$ (14.5)

Enthalpy is a state function as well. If we are interested in computing changes in enthalpy, we write:

$H=U+P\stackrel{˜}{v}$ (14.6)

Since we already have a way to compute dU (equation 14.4), we can now write the fundamental thermodynamic relationship used to compute changes of enthalpy in a closed system:

$dH=dU=PdV=VdP$

By the same token, the formal, fundamental definition of Gibbs Free Energy (G) is the following:

$G=H=TS$ (14.8)

from which we see that changes in this property may be calculated as:

$dG=dH=TdS=SdT$ (14.9a)

We now substitute (14.7) to get the fundamental thermodynamic relationship for changes in Gibbs free energy in a closed system:

$dG=VdP=SdT$ (14.9b)

We proceed with Helmholtz Free Energy (A) accordingly. Its definition:

$A=U=TS$ (14.10)

The expression of change,

$dA=dU=TdS=SdT$ (14.11)

We substitute (6.4) into it and get:

$dA=-PdV-SdT$ (14.12)

We have just derived the following fundamental thermodynamic relationships for fluids of constant composition:

$dU=TdS-PdV$ (14.4)

$dH=TdS+VdP$ (14.7)

$dG=VdP-SdT$ (14.9)

$dA=-PdV-SdT$ (14.12)

All these properties (U, H, G and A) are state functions and extensive properties of the system. Notice that the only assumption that we took throughout their development is that the system was closed (for the derivation of 14.4). Hence, these equations strictly apply to systems of constant composition.

# Functions of State or State Functions

A function of state is one in which the differential change is determined only by the end states and not by intervening states. Most thermodynamic variables are state functions and hence property changes are determined by the end states and not by the process path. Notable exceptions are work and heat.

The most common thermodynamic state functions with which we shall deal include

• Internal Energy (U),
• Entropy (S),
• Enthalpy (H),
• Helmotz Energy (A), and
• Gibbs Energy (G).

# Mechanics of Manipulating a Function of State

Given that f(x,y,z) is any state function that characterizes the system and (x,y,z) is a set of independent variable properties of that system, we know that any change $\Delta$ f will be only a function of the value of “f” at the final and initial states,

$\Delta f={f}_{2}-{f}_{1}=f\left({x}_{2},{y}_{2},{z}_{2}\right)-f\left({x}_{1},{y}_{1},{z}_{1}\right)$ (14.13)

Since f=f(x,y,z), we can mathematically relate the total differential change (df) to the partial derivatives $\frac{\partial f}{\partial x}$ ,$\frac{\partial f}{\partial y}$ , and $\frac{\partial f}{\partial z}$ of the function, as follows:

$df={\left(\frac{\partial f}{\partial x}\right)}_{y,z}dx+{\left(\frac{\partial f}{\partial y}\right)}_{x,z}dy+{\left(\frac{\partial f}{\partial z}\right)}_{x,y}dz$ (14.14)

where, in general:

${\left(\frac{\partial f}{\partial x}\right)}_{y,z}=$ the change of f with respect to x, while y and z are unchanged.

If we want to come up with the total change, $\Delta f$ , of a property (we want to go from 14.14. to 14.13), we integrate the expression in (14.14) to get:

$\Delta f={f}_{2}-{f}_{1}=\underset{{x}_{1}}{\overset{{x}_{2}}{\int }}{\left(\frac{\partial f}{\partial x}\right)}_{y,z}dx+\underset{{y}_{{}_{1}}}{\overset{{y}_{2}}{\int }}{\left(\frac{\partial f}{\partial y}\right)}_{x,z}dy+\underset{{z}_{1}}{\overset{{z}_{2}}{\int }}{\left(\frac{\partial f}{\partial x}\right)}_{x,y}dz$ (14.15)

Let us visualize this with an example. For a system of constant composition, its thermodynamic state is completely defined when two properties of the system are fixed. Let us say we have a pure component at a fixed pressure (P) and temperature (T). Hence, all other thermodynamic properties, for example, enthalpy (H), are fixed as well. Since H is only a function of P and T, we write:

$H=H\left(P,T\right)$ (14.16)

and hence, applying 6.2, any differential change in enthalpy can be computed as:

$dH={\left(\frac{\partial H}{\partial P}\right)}_{T}dP+{\left(\frac{\partial H}{\partial T}\right)}_{P}dT$ (14.17)

The total change in enthalpy of the pure-component system becomes:

$\Delta H={H}_{2}-{H}_{1}=\underset{{P}_{1}}{\overset{{P}_{2}}{\int }}{\left(\frac{\partial H}{\partial P}\right)}_{T}dP+\underset{{T}_{1}}{\overset{{T}_{2}}{\int }}{\left(\frac{\partial H}{\partial T}\right)}_{P}dT$ (14.18)

Now we are ready to spell out the exactness condition, which is the mathematical condition for a function to be a state function. The fact of the matter is, that for a function to be a state function — i.e., its integrated path shown in (14.15) is only a function of the end states, as shown in (14.13) — its total differential must be exact. In other words, if the total differential shown in (14.14) is exact, then f(x,y,z) is a state function. How do we know if a total differential is exact or not?

Given a function $\Psi$ (x,y,z),

$d\Psi =M\left(x,y,z\right)dx+N\left(x,y,z\right)dy+Q\left(x,y,z\right)dz$ (14.19a)

where:

$M\left(x,y,z\right)={\left(\frac{\partial \Psi }{\partial x}\right)}_{y,z}$ (14.19b)

$N\left(x,y,z\right)={\left(\frac{\partial \Psi }{\partial y}\right)}_{x,z}$ (14.19c)

$Q\left(x,y,z\right)={\left(\frac{\partial \Psi }{\partial z}\right)}_{x,y}$ (14.19d)

we say that $d\Psi$ is an exact differential and consequently $\psi \left(x,y,z\right)$ a state function if all the following conditions are satisfied:

${\left(\frac{\partial M}{\partial y}\right)}_{x,z}={\left(\frac{\partial N}{\partial x}\right)}_{y,z}$ (14.20a)

${\left(\frac{\partial N}{\partial z}\right)}_{x,y}={\left(\frac{\partial Q}{\partial y}\right)}_{x,z}$ (14.20b)

${\left(\frac{\partial M}{\partial z}\right)}_{x,y}={\left(\frac{\partial Q}{\partial x}\right)}_{y,z}$ (14.20c)

Equations (14.20) are called the exactness condition.

# Computing Property Changes in Closed Systems

Now we recall the fundamental thermodynamic relationships for closed systems, which we derived above:

(14.4)

(14.7)

(14.9)

(14.12)

where U, H, G, and A are state functions. There is something remarkable about the above expressions: they allow for the direct calculation of the change in a state function as a function of changes in other two. An important lesson to be learned here: when dealing with thermodynamic state properties, we are most interested in changes in the value of state properties rather than their actual values.

As we have said, the above relationships allow us to visualize that the changes of each thermodynamic property are dependent on the change of two others in a closed system.

(14.21a)

(14.21b)

(14.21c)

(14.21d)

Therefore, recalling what we did in equation (14.14), we express the total differential of each of these properties as:

$dU={\left(\frac{\partial U}{\partial S}\right)}_{v}dS+{\left(\frac{\partial U}{\partial V}\right)}_{s}dV$
(14.22a)

$dH={\left(\frac{\partial H}{\partial S}\right)}_{p}dS+{\left(\frac{\partial H}{\partial P}\right)}_{s}dP$
(14.22b)

$dG={\left(\frac{\partial G}{\partial P}\right)}_{r}dP+{\left(\frac{\partial G}{\partial T}\right)}_{p}dT$
(14.22c)

$dA={\left(\frac{\partial A}{\partial V}\right)}_{r}dV+{\left(\frac{\partial A}{\partial T}\right)}_{v}dT$
(14.22d)

Comparing these equations term to term, the following realizations can be made:

$T={\left(\frac{\partial U}{\partial S}\right)}_{v};P=-{\left(\frac{\partial U}{\partial V}\right)}_{s}$
(14.23a)

$T={\left(\frac{\partial H}{\partial S}\right)}_{p};V=-{\left(\frac{\partial H}{\partial P}\right)}_{s}$
(14.23b)

$V={\left(\frac{\partial G}{\partial P}\right)}_{r};S=-{\left(\frac{\partial G}{\partial T}\right)}_{p}$
(14.23c)

$P={\left(\frac{\partial A}{\partial V}\right)}_{r};S=-{\left(\frac{\partial A}{\partial T}\right)}_{v}$
(14.23d)

Additionally, we can go one step further. U, H, G, and A are state functions, and as such, their total differentials (equations 14.4, 14.7, 14.9, and 14.12) must be exact. Recall that for a total differential df=Mdx+Ndy to be an exact differential, it must satisfy the equation:

${\left(\frac{\partial M}{\partial y}\right)}_{x}={\left(\frac{\partial N}{\partial x}\right)}_{y}$
(14.24)

Equation (14.24) is the exactness criteria for a function of two independent variables. It was previously stated in (14.20) for a state function of three independent variables. Its application yields:

${\left(\frac{\partial T}{\partial V}\right)}_{s}={\left(\frac{\partial P}{\partial S}\right)}_{v}$
(14.25a)

${\left(\frac{\partial T}{\partial P}\right)}_{s}={\left(\frac{\partial V}{\partial S}\right)}_{p}$
(14.25b)

${\left(\frac{\partial V}{\partial T}\right)}_{p}={\left(\frac{\partial S}{\partial P}\right)}_{r}$
(14.25c)

${\left(\frac{\partial P}{\partial T}\right)}_{v}={\left(\frac{\partial S}{\partial V}\right)}_{r}$
(14.25d)

Equations (14.25) are known as Maxwell’s relationships. Maxwell’s relationships are very useful for manipulating thermodynamic equations. For instance, it is always desirable for practical purposes to express thermodynamic properties such as enthalpy (H) and entropy (S) as functions of measurable properties such as pressure (P) and temperature (T):

(14.26a)

(14.26a)

Starting with the total differential of H and S as a function of P and T, we can prove that the relationships between the parameters H and S and the parameters P and T are given by the expressions:

$dH={C}_{p}dT+\left[V-T{\left(\frac{\partial V}{\partial T}\right)}_{p}\right]dP$
(14.27a)

$dS={C}_{p}\frac{dT}{T}-{\left(\frac{\partial V}{\partial T}\right)}_{p}dP$
(14.27b)

Additionally, since these expressions for dH and dS must also be exact differentials, you could prove that the heat capacity at constant pressure (CP) for an ideal gas (PV = RT) does not depend on pressure. The thermodynamic definition of CP and Cv are:

(14.28a)

(14.28b)

The heat capacity at constant volume (CV) for an ideal gas does not depend on pressure, either. You could prove this by proving first that Cp=CV+R (R=universal gas constant) for ideal gases, using above tools.

# Action Item

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

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

#### Problem Set

1. In you own words, state how you would explain the implications of the 1st law and 2nd law of thermodynamics. Assume that you need to explain this to a 5 year old child that has no idea about thermodynamic jargon.
2. Suppose that most thermodynamic properties were not state functions. What impact would that have in process design?

# Thermodynamic Tools (II)

#### Module Goal:

To establish the mathematical framework for thermodynamics of phase equilibrium.

#### Module Objective:

To establish that there is an unique relationship between partial molar quantities in any mixture.

# Homogeneous Functions, Euler's Theorem and Partial Molar Quantities

Any function f(x) that possesses the characteristic mapping:

$\begin{array}{l}x\to \lambda x\\ f\left(x\right)\to \lambda f\left(x\right)\end{array}$
(15.1)

is said to be homogeneous, with respect to x, to degree 1. By the same token, if f(x) obeys the mapping:

$\begin{array}{l}x\to \lambda x\\ f\left(x\right)\to {\lambda }^{k}f\left(x\right)\end{array}$
(15.2)

then f(x) is homogeneous to degree “k”. In general, a multivariable function f(x1,x2,x3,…) is said to be homogeneous of degree “k” in variables xi(i=1,2,3,…) if for any value of $\lambda$ ,

$f\left(\lambda {x}_{1},\lambda {x}_{2},...\right)={\lambda }^{k}f\left({x}_{1},{x}_{2},...\right)$
(15.3)

For example, let us consider the function:

$f\left(x,y\right)=\frac{x}{{x}^{2}+{y}^{2}}$
(15.4)

How do we find out if this particular function is homogeneous, and if it is, to what degree? We evaluate this function at to obtain:

$f\left(\lambda x,\lambda y\right)=\frac{\lambda x}{{\lambda }^{2}{x}^{2}+{\lambda }^{2}{y}^{2}}={\lambda }^{-1}\frac{x}{{x}^{2}+{y}^{2}}{\lambda }^{-1}f\left(x,y\right)$
(15.5)

hence, the function f(x,y) in (15.4) is homogeneous to degree -1.

In regard to thermodynamics, extensive variables are homogeneous with degree “1” with respect to the number of moles of each component. They are, in fact, proportional to the mass of the system to the power of one (k=1 in equation 15.2 or 15.3). This is, if we triple the amount of mass in the system, the value of any given extensive property will be tripled as well. Notice that this is not the case for intensive properties of the system (such as temperature or pressure), simply because they are independent of mass. Hence, intensive thermodynamic properties are homogeneous functions with degree “0” — in such a case, k=0 in equation (15.2) or (15.3).

From the previous section, it is clear that we are not only interested in looking at thermodynamic functions alone, but that it is also very important to compute how thermodynamic functions change and how that change is mathematically related to their partial derivatives . Hence, to complete the discussion on homogeneous functions, it is useful to study the mathematical theorem that establishes a relationship between a homogeneous function and its partial derivatives. This is Euler’s theorem.

Euler’s theorem states that if a function f(ai, i = 1,2,…) is homogeneous to degree “k”, then such a function can be written in terms of its partial derivatives, as follows:

$k{\lambda }^{k-1}f\left({a}_{i}\right)=\sum _{i}^{}{a}_{i}{\left(\frac{\partial f\left({a}_{i}\right)}{\partial \left(\lambda {a}_{i}\right)}\right)|}_{\lambda x}$
15.6a

Since (15.6a) is true for all values of $\lambda$ , it must be true for $\lambda -1$ . In this case, (15.6a) takes a special form:

$kf\left({a}_{i}\right)=\sum _{i}^{}{a}_{i}{\left(\frac{\partial f\left({a}_{i}\right)}{\partial \left({a}_{i}\right)}\right)|}_{x}$
15.6b

So far, so good. But…what is the application of all this? Well, first of all, we have to know something more about extensive thermodynamic properties. A very neat thing about them is that they can be written as a function of a sufficient number of independent variables to completely define the thermodynamic state of the system. Such a set is said to be a complete set. As it turns out, any thermodynamic system is completely defined when both the masses of all the substances within it are defined and two additional independent variables are fixed. This is Duhem’s theorem. From a real-life perspective, it is natural to choose pressure and temperature as those “independent variables” — physical quantities that we have a “feel” for and we think we can control — rather than specific volume or entropy. As we will see later, they are also convenient variables of choice because they are homogeneous of degree zero in mass.

Let “$\Im$ ” be a given extensive property of a multi-component system. From the previous section, we know that the value of “$\Im$ ” must be fixed and uniquely determined once we fix the pressure, temperature, and number of moles of each component in the system. This is,

$\Im =\Im \left(P,T,{n}_{1},{n}_{2},...,{n}_{N}\right)$
15.7a

Additionally, we recall that extensive properties are homogeneous of degree one with respect to number of moles and homogeneous of degree zero with respect to pressure and temperature. Thus, expression (15.6b) is readily applicable:

$\Im =\sum _{i}{n}_{i}{\left(\frac{\partial \Im }{\partial {n}_{i}}\right)}_{P,T,{n}_{i=1}}=\sum _{i}{n}_{i}{\Im }_{i}$
15.7b

where we have just defined:

${\Im }_{i}{\left(\frac{\partial \Im }{\partial {n}_{i}}\right)}_{P,T,{n}_{i=1}}$
15.7c

Equation (15.7c) is a very important definition. It defines the concept of a partial molar quantity. A partial molar quantity $\overline{{\Im }_{i}}$ represents the change in the total quantity ($\Im$ ) due to the addition of an infinitesimal amount of species “i” to the system at constant pressure and temperature. Euler’s theorem gave birth to the concept of partial molar quantity and provides the functional link between it (calculated for each component) and the total quantity. The selection of pressure and temperature in (15.7c) was not trivial. First, they are convenient variables to work with because we can measure them in the lab. But most important, they are intensive variables, homogeneous functions of degree zero in number of moles (and mass). This allowed us to use Euler’s theorem and jump to (15.7b), where only a summation with respect to number of moles survived. The definition of the partial molar quantity followed.

The conventional notation we are going to follow throughout the following section is:

$\Im$ = Total quantity (e.g., total volume, total internal energy, etc),

$\stackrel{˜}{\Im }$ = Molar quantity, i.e., total quantity per unit mole:

15.8a

$\overline{\Im }$ = Partial molar quantity,

$\stackrel{^}{\Im }$ = Mass or specific quantity, i.e., total quantity per unit mass:

15.8b

We can rewrite equation (15.7b) in terms of molar quantity using the definition in (15.8a),

$\stackrel{˜}{\Im }=\frac{1}{n}\sum _{i}{n}_{i}{\overline{\Im }}_{i}=\sum _{i}{x}_{i}{\overline{\Im }}_{i}$
(15.9)

where:

Any molar quantity in thermodynamics can be written in terms of the partial molar quantity of its constituents. If we set $\Im =G$ , we end up with,

$\stackrel{˜}{G}=\sum _{i}{x}_{i}{\overline{G}}_{i}$
(15.10a)

Equivalently, if we set $\Im =V$ (total volume),

$\stackrel{˜}{v}=\sum _{i}{x}_{i}{\overline{v}}_{i}$
(15.10b)

Notice that for single component systems (xi=1), partial molar properties are just equal to the molar property:

E
(15.11a)

(15.11b)

This is also a consequence of the definition in (15.7c),

${\overline{\Im }}_{i}={\left(\frac{\partial \Im }{\partial {n}_{i}}\right)}_{P,T,{n}_{i=1}}$
(15.7c)

For a pure component, and:

$\overline{\Im }={\left(\frac{\partial n\stackrel{˜}{\Im }}{\partial n}\right)}_{PI}=\stackrel{˜}{\Im }$
(15.12)

The reason for the introduction of the concept of a partial molar quantity is that often times we deal with mixtures rather than pure-component systems. The way to characterize the state of the mixtures is via partial molar properties. This concept provides the bridge between the thermodynamics of systems of constant composition, which we have studied so far, and the thermodynamics of systems of variable composition, which we will deal with in the next section. Basically, the definition in (15.7c):

${\overline{\Im }}_{i}={\left(\frac{\partial \Im }{\partial {n}_{i}}\right)}_{P,T,{n}_{i=1}}$
(15.7c)

allows us to quantify how the total, extensive property $\Im$ changes with additions of ni at constant pressure and temperature. If you look at (15.7b) and (15.9), you will also realize that (15.7c) is just an allocation formula that allows assigning to each species “i” a share of the total mixture property, such that:

$\Im =\sum _{i}{n}_{i}{\overline{\Im }}_{i}$
(15.7b)

$\stackrel{˜}{\Im }=\sum _{i}{x}_{i}{\overline{\Im }}_{i}$
(15.9)

We can play with “$\Im$ ” a little more. Let us say that we are now interested in looking at the differential changes of $\Im$. Since we also know that $\Im$ is a state function, and given the functional relationship in (15.7a), the total differential for $\Im$ is written:

$d\Im ={\left(\frac{\partial \Im }{\partial P}\right)}_{T{n}_{i}}dP+{\left(\frac{\partial \Im }{\partial T}\right)}_{P,{n}_{i}}dT+\sum _{i}{\left(\frac{\partial \Im }{\partial {n}_{i}}\right)}_{P,T,n=1}d{n}_{i}$
(15.13a)

or

$d\Im ={\left(\frac{\partial \Im }{\partial P}\right)}_{T{n}_{i}}dP+{\left(\frac{\partial \Im }{\partial T}\right)}_{P,{n}_{i}}dT+\sum _{i}{\overline{\Im }}_{i}d{n}_{i}$
(15.13b)

Basically, equations (15.13) tell us that any change in P, T, or ni will cause a corresponding change in the total property, $\Im$ . This is a reinforcement of what is explicitly declared in (15.7a). If we recall (15.7b), an alternate expression for the total differential in (15.13) is written:

$d\Im =d\left(\sum _{i}{n}_{i}{\overline{\Im }}_{i}\right)=\sum _{i}{\overline{\Im }}_{i}d{n}_{i}+\sum _{i}{n}_{i}d{\overline{\Im }}_{i}$
(15.14)

If we subtract (15.14) from (15.13b), we get:

$d\Im -d\Im ={\left(\frac{\partial \Im }{\partial P}\right)}_{T,{n}_{i}}dP+{\left(\frac{\partial \Im }{\partial T}\right)}_{P,{n}_{i}}dT+\sum _{i}{\overline{\Im }}_{i}d{n}_{i}-\sum _{i}{\overline{\Im }}_{i}d{n}_{i}-\sum _{i}{n}_{i}d{\overline{\Im }}_{i}$
(15.15)

Therefore,

${\left(\frac{\partial \Im }{\partial P}\right)}_{T,{n}_{i}}dP+{\left(\frac{\partial \Im }{\partial T}\right)}_{P,{n}_{i}}dT-\sum _{i}{n}_{i}d{\overline{\Im }}_{i}=0$
(15.16)

Equation (15.16) is the well-known Gibbs-Duhem equation. It can be applied to any extensive thermodynamic property: U, S, H, G, A, and it must hold true. It represents a thermodynamic constraint between the intensive variables P, T and $\stackrel{˜}{\Im }$ . Pressure, temperature and partial molar properties cannot vary in just any fashion; any change taking place among them must satisfy (15.16). The change in any one of them can be calculated as a function of the change in the other two by means of the Gibbs-Duhem equation. This equation is the basis for thermodynamic consistency checks of experimental data.

# Thermodynamics of Systems of Variable Composition (Open Multicomponent Systems)

To extend all of the previous concepts to systems of variable mass, we must now consider at least one new variable: number of moles, n. To account for this effect, you will see below that we will have to introduce a “new” thermodynamic property. Let us take the case of the internal energy. For a constant composition system, we wrote:

$dU=TdS-PdV$
(14.4)

and,

$dU={\left(\frac{\partial U}{\partial S}\right)}_{V}dS+{\left(\frac{\partial U}{\partial V}\right)}_{S}dV$
(14.22a)

If this system has a single component, and now we allow its mass (and hence, number of moles, n) to change (an open system), the change in U (dU) is no longer just a function of dS and dV. We now have to account for changes in ‘n’, thus:

$dU={\left(\frac{\partial U}{\partial S}\right)}_{V,n}dS+{\left(\frac{\partial U}{\partial V}\right)}_{S,n}dV+{\left(\frac{\partial U}{\partial n}\right)}_{S,V}dn$
(15.17)

If we had a binary system, we would have two new variables, ‘n1’ and ‘n2’ (n=n1+n2) and we would have to expand (15.17) accordingly,

$dU={\left(\frac{\partial U}{\partial S}\right)}_{V,n}dS+{\left(\frac{\partial U}{\partial V}\right)}_{S,n}dV+{\left(\frac{\partial U}{\partial {n}_{1}}\right)}_{S,V,{n}_{2}}d{n}_{1}+{\left(\frac{\partial U}{\partial {n}_{2}}\right)}_{S,V,{n}_{1}}d{n}_{2}$
(15.18)

Hence, for a multi-component system, we just keep on adding terms ($n=\sum _{i}{n}_{i}$ ):

$dU={\left(\frac{dU}{dS}\right)}_{V,n}dS+{\left(\frac{dU}{dV}\right)}_{S,n}dV+\sum _{i}{\left(\frac{\partial \Im }{\partial {n}_{i}}\right)}_{S,V,{n}_{i}\ne {n}_{1}}d{n}_{i}$
(15.19)

Thermodynamics defines this ‘coefficient’ which multiplies the change in the number of moles of each component (dni) as the ‘chemical potential’ of that component(μi). See how the chemical potential is a thermodynamic property that must be defined for the proper description of a system of variable composition — that is, an open system.

Then, we write:

${\mu }_{i}={\left(\frac{dU}{d{n}_{i}}\right)}_{{}_{S,V,{n}_{i}\ne {n}_{1}}}$
(15.20)

and finally,

$dU=TdS-PdV+\sum _{i}{\mu }_{i}d{n}_{i}$
(15.21)

The same can be done with the thermodynamic definitions for dH, dG, and dA (the rest of equations 14.22). In fact, the chemical potential may be defined in at least four different and equivalent ways. You can now show that for equations (14.22) to account for systems of variable composition, we would have to expand them into:

$dH={\left(\frac{\partial H}{\partial S}\right)}_{P,n}dS={\left(\frac{\partial H}{\partial P}\right)}_{S,n}dP+{\sum _{i}\left(\frac{\partial H}{\partial {n}_{i}}\right)}_{S,P,n}{}_{{}_{i}\ne {n}_{1}}d{n}_{i}$
(15.22a)

$dG={\left(\frac{\partial G}{\partial P}\right)}_{T,n}dP={\left(\frac{\partial G}{\partial T}\right)}_{P,n}dT+{\sum _{i}\left(\frac{\partial G}{\partial {n}_{i}}\right)}_{P,T,n}{}_{{}_{i}\ne {n}_{1}}d{n}_{i}$
(15.22b)

$dA={\left(\frac{\partial A}{\partial V}\right)}_{T,n}dS={\left(\frac{\partial A}{\partial T}\right)}_{V,n}dT+{\sum _{i}\left(\frac{\partial A}{\partial {n}_{i}}\right)}_{T,V,n}{}_{{}_{i}\ne {n}_{1}}d{n}_{i}$
(15.22c)

Because thermodynamics defines the ‘coefficient’ which multiplies the change in the number of moles of each component (dni) as the ‘chemical potential’ of that component(${\mu }_{i}$ ), we have three new ways to express it:

${\mu }_{i}={\left(\frac{\partial H}{\partial {n}_{i}}\right)}_{S,P,n}{}_{{}_{i}\ne {n}_{j}}$
(15.23a)

${\mu }_{i}={\left(\frac{\partial G}{\partial {n}_{i}}\right)}_{P,T,n}{}_{{}_{i}\ne {n}_{j}}$
(15.23b)

${\mu }_{i}={\left(\frac{\partial A}{\partial {n}_{i}}\right)}_{T,V,n}{}_{{}_{i}\ne {n}_{j}}$
(15.23c)

Don’t any of these three equations ring a bell? Don’t forgot the definition of (15.7c). Of all the four available definitions for chemical potential, there is one that suits our definition of a partial molar quantity perfectly. Compare equation (15.23b) to (15.7c). What we can say is that the chemical potential of a component “i” is equal to the partial molar Gibbs energy of such component:

${\mu }_{i}={\left(\frac{\partial G}{\partial {n}_{i}}\right)}_{P,T,n}{}_{{}_{i}\ne {n}_{j}}={\overline{G}}_{l}$
(15.24)

Notice that, for a pure component, the chemical potential is equal to the molar Gibbs energy of the substance (see equation 15.9),

E$\mu =\stackrel{˜}{G}$
(15.25)

When we introduce the definition of chemical potential into each of equations (15.22), the fundamental thermodynamic expressions that apply to systems of variable composition become:

$dU=TdS-PdV+\sum _{i}{\mu }_{i}d{n}_{i}$
(15.26a)

$dH=TdS+VdP+\sum _{i}{\mu }_{i}d{n}_{i}$
(15.26b)

$dG=VdP-SdT+\sum _{i}{\mu }_{i}d{n}_{i}$
(15.26c)

$dA=-PdV-SdT+\sum _{i}{\mu }_{i}d{n}_{i}$
(15.26d)

It is clear that equations (14.23) still hold — evaluated at constant “n” (total number of moles) — as shown:

$T={\left(\frac{\partial U}{\partial S}\right)}_{V,n}P=-{\left(\frac{\partial U}{\partial V}\right)}_{S,n}$
(14.23a)

$T={\left(\frac{\partial H}{\partial S}\right)}_{P,n}V={\left(\frac{\partial H}{\partial P}\right)}_{S,n}$
(14.23b)

$V={\left(\frac{\partial G}{\partial P}\right)}_{T,n}S=-{\left(\frac{\partial G}{\partial T}\right)}_{P,n}$
(14.23c)

$P=-{\left(\frac{\partial A}{\partial V}\right)}_{T,n}S=-{\left(\frac{\partial A}{\partial T}\right)}_{V,n}$
(14.23d)

The same reasoning applies to Maxwell’s relationships. From the above expressions, equation (15.22b) is the only one that matches equation (15.13). Notice that the following additional identities can be also identified (see equations 14.20):

${\left(\frac{\partial {\mu }_{i}}{\partial S}\right)}_{V,n}={\left(\frac{\partial T}{\partial {n}_{i}}\right)}_{S,V,{n}_{i\ne 1}}{\left(\frac{\partial {\mu }_{i}}{\partial V}\right)}_{S,n}={\left(\frac{\partial P}{\partial {n}_{i}}\right)}_{S,V,{n}_{i\ne 1}}$
(15.27a)

${\left(\frac{\partial {\mu }_{i}}{\partial S}\right)}_{P,n}={\left(\frac{\partial T}{\partial {n}_{i}}\right)}_{S,P,{n}_{i\ne 1}}{\left(\frac{\partial {\mu }_{i}}{\partial P}\right)}_{S,n}={\left(\frac{\partial V}{\partial {n}_{i}}\right)}_{S,P,{n}_{i\ne 1}}$
(15.27b)

${\left(\frac{\partial {\mu }_{i}}{\partial P}\right)}_{T,n}={\left(\frac{\partial V}{\partial {n}_{i}}\right)}_{P,T,{n}_{i\ne 1}}={\overline{V}}_{i}{\left(\frac{\partial {\mu }_{i}}{\partial T}\right)}_{P,n}={\left(\frac{\partial S}{\partial {n}_{i}}\right)}_{P,T,{n}_{i\ne 1}}={\overline{S}}_{i}$
(15.27c)

${\left(\frac{\partial {\mu }_{i}}{\partial V}\right)}_{T,n}={\left(\frac{\partial P}{\partial {n}_{i}}\right)}_{T,V,{n}_{i\ne 1}}{\left(\frac{\partial {\mu }_{i}}{\partial T}\right)}_{V,n}={\left(\frac{\partial S}{\partial {n}_{i}}\right)}_{T,V,{n}_{i\ne 1}}$
(15.27d)

We can see that chemical potential can be calculated by solving any of these differential expressions. To do this, experimental information is needed on how the other properties (T, V, S, P) change with additions of the given species (ni) under certain restraining conditions.

In some instances, we may have an open system of constant composition. This is the case of a system of a pure component exchanging mass with its surroundings. For such a system, nc = 1 and equations (15.26) are written as:

$dU=TdS-PdV+\mu dn$
(15.28a)

$dH=TdS+VdP+\mu dn$
(15.28b)

$dG=VdP-SdT+\mu dn$
(15.28c)

$dA=-PdV-SdT+\mu dn$
(15.28d)

Equations (14.23) and (15.27) will still hold with .

# Action Item

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

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

#### Problem Set

1. For a design project involving a binary C1/C2 system, an engineer needs to know the partial molar volumes of C1 and C2 at the same conditions. Due to budget constraints, he was able to measure only one of them. Speculate on how to go about obtaining the second one.

# Thermodynamic Tools (III)

### Module Goals

Module Goal: To establish the mathematical framework for thermodynamics of phase equilibrium.

Module Objective: To establish that there is an unique relationship between partial molar quantities in any mixture.

# The Chemical Potentials

We have just seen that the chemical potential is a thermodynamic property which is related to all thermodynamic properties with units of energy. Its most useful definition is given in terms of constant pressure and temperature:

${\mu }_{i}={\left(\frac{\partial G}{\partial {n}_{i}}\right)}_{P,T,{n}_{i}\ne {n}_{1}}={\overline{G}}_{i}$
(15.24)

This constitutes the working definition of chemical potential, the one that relates it to the partial molar quantity concept studied before.

Although the mathematical definition of chemical potential can be stated clearly, its “physical” meaning is not as easy to grasp. It is much easier to understand the physical implications of “pressure,” “temperature,” and “internal energy” than it is to undestand the physical interpretation of “chemical potential.” Given equation (15.24), and recognizing Gibbs free energy (G) as the capacity of a system to do work, we may write the following formal definition of chemical potential:

The chemical potential of a component in a given phase is the rate of increase of the capacity of the phase to do work per unit addition of the substance to the phase, at constant temperature and pressure.

We may also quote the definition that J.W. Gibbs provided for it:

If to any homogeneous mass, we suppose an infinitesimal quantity of any substance to be added, the mass remaining homogeneous and its entropy and volume remaining unchanged, the increase of energy of the mass divided by the quantity of the substance added is the potential for that substance in the mass considered.

This definition is closely related to the mathematical definition given in (15.20).

${\mu }_{i}={\left(\frac{\partial U}{\partial {n}_{i}}\right)}_{S,V,{n}_{i}\ne {n}_{1}}$
(15.20)

To understand the physical implications of the chemical potential of a species, we have to recall that for any thermodynamic process to be carried out, a driving force must be causing it. For instance, a pressure gradient is the driving force that causes the bulk movement of fluids from one point to the other, and a temperature gradient provides the potential difference needed for heat to flow. We also know that if we have a higher concentration of solute in a homogeneous system, it will diffuse to the zones of lower concentration. Here, the chemical potential is responsible for the diffusion of species within two points in space, or even its exchange between two different phases, without the presence of either pressure or temperature gradients. The chemical potential is the potential describing the ability of species to move from one phase to another.

# The Thermodynamic Concept of Equilibrium

Intuitively, the concept of equilibrium conveys the message that something “balances out.” Equilibrium describes a state of vanishing driving forces or gradients, where everything remains as it is. If a system is in equilibrium, it retains its current state because there are no driving forces causing anything to change.

If two materials have the same temperature, we say that they are in thermal equilibrium. No exchange of heat takes places because there are no thermal gradients. For instance, a liquid and a vapor phase are in thermal equilibrium when:

${T}_{l}=Tv$
(16.1)

We achieve mechanical equilibrium if two substances are found at the same pressure. No bulk movement of fluids takes place because there are no pressure gradients. A liquid and a vapor phase are in mechanical equilibrium when:

${P}_{l}=Pv$
(16.2)

For a thermodynamic system to be in equilibrium, all intensive (temperature, pressure) and extensive thermodynamic properties (U, G, A, H, S, etc) must be constants. Hence, the total change in any of those properties ($d\Im$ ) must be zero at equilibrium.

Now we would like to have a concept of thermodynamic equilibrium for a vapor-liquid equilibrium. Let us consider a closed, heterogeneous vapor-liquid system. Any changes in a total property of the system will be the result of the changes of that property in the liquid phase plus the changes of that property in the vapor phase.

$d{\Im }^{\left(total\right)}=d{\Im }^{\left(liquid\right)}+d{\Im }^{\left(vapor\right)}$
(16.3)

In this case, liquid and vapor by themselves are not closed systems; they can exchange matter between themselves but not with the surroundings. To elaborate more upon the concept of equilibrium, let’s look at equation (15.26c). Because it is written in terms of changes in pressure and temperature, two measurable laboratory quantities, it is the “friendliest” of all fundamental equations. We write it for both of the phases:

$d{\Im }^{\left(liquid\right)}={\left(VdP-SdT\right)}^{l}+\sum _{i}{\mu }_{i}^{l}d{n}_{i}^{l}$
(16.4a)

$d{G}^{\left(vapor\right)}={\left(VdP-SdT\right)}^{v}+\sum _{i}{\mu }_{i}^{v}d{n}_{i}^{v}$
(16.4b)

In (16.3), for $\Im =G$ , we get:

$d{G}^{\left(total\right)}=d{G}^{\left(liquid\right)}+d{G}^{\left(vapor\right)}$
(16.5)

Hence,

$d{G}^{\left(total\right)}={\left(VdP-SdT\right)}^{l}+{\left(VdP-SdT\right)}^{v}+\sum _{i}{\mu }_{i}^{l}d{n}_{i}^{l}+\sum _{i}{\mu }_{i}^{v}d{n}_{i}^{v}$
(16.16)

Since at equilibrium all extensive properties, such as G, must remain constant, dG(total) must be zero. For this to hold true, and by inspection of equation (16.6), the conditions for thermodynamic equilibrium are:

(16.7)

E
(16.8)

(16.9)

It can be also proven that, at equilibrium, the total free energy of the system (G(total)) must take a minimum value; this reinforces the fact that dG(total)=0 at equilibrium. The minimum Gibbs energy criterion for equilibrium is a restatement of the second law of thermodynamics, from which we know that the entropy of a system in equilibrium must be at its maximum, considering all of the possible states for equilibrium.

It is somehow reasonable that for a true equilibrium condition there should be neither pressure nor temperature gradients (equations 16.7 and 16.8). This is because equilibrium is, at the very least, a state of lack of gradients. But what is equation (16.9) trying to tell us? To demystify equation (16.9), we recall that we are dealing with a closed system, hence, the total amount of moles per species:

${n}_{i}^{\left(total\right)}={n}_{i}^{\left(l\right)}+{n}_{i}^{\left(v\right)}$
(16.10)

must be constant (we do not allow for chemical reactions within the system). Thus we write:

$d{n}_{i}^{\left(total\right)}=d{n}_{i}^{\left(l\right)}+d{n}_{i}^{\left(v\right)}=0$
(16.11)

Therefore,

$d{n}_{i}^{\left(v\right)}=-d{n}_{i}^{\left(l\right)}$
(16.12)

(16.12) into (16.11) yields:

$\sum _{i}\left({\mu }_{i}^{\left(l\right)}-{\mu }_{i}^{\left(v\right)}\right)d{n}_{i}^{\left(v\right)}=0$
(16.13)

For equation (16.13) to hold true,

(16.14)

We have then arrived at the criteria for vapor-liquid equilibria for a system at constant pressure and temperature: the chemical potential of every species must be the same in both phases. We may generalize this finding to any number of phases, for which the chemical potential of every species must be the same in all phases. The chemical potential being the driving force which moves a species from one phase to the other, equation (16.14) is physically reasonable. If the chemical potential of a species in one phase is the same as that in the other, there is zero driving force and thus a zero net transfer of species at equilibrium.

# Fugacity

We have seen that, for a closed system, the Gibbs energy is related to pressure and temperature as follows:

$dG=VdP-SdT$
(16.15)

For a constant temperature process,

(16.16)

For an ideal gas,

$dG=\frac{RT}{P}dP$
(16.17)
(16.18)

This expression by itself is strictly applicable to ideal gases. However, Lewis, in 1905, suggested extending the applicability of this expression to all substances by defining a new thermodynamic property called fugacity, f, such that:

(16.19)

This definition implies that for ideal gases, ‘f’ must be equal to ‘P’. For mixtures, this expression is written as:

(16.20)

where are the partial molar Gibbs energy and fugacity of the i-th component, respectively. Fugacity can be readily related to chemical potential because of the one-to-one relationship of Gibbs energy to chemical potential, which we have discussed previously. Therefore, the definition of fugacity in terms of chemical potential becomes:

For a pure substance,

(16.21a)
(16.21b)

For a component in a mixture,

(16.22c)
(16.22d)

The fugacity coefficient (${\varphi }_{i}$ ) is defined as the ratio of fugacity to its value at the ideal state. Hence, for pure substances:

$\varphi =\frac{f}{P}$
(16.23a)

and for a component in a mixture,

${\varphi }_{i}=\frac{{f}_{i}}{{y}_{i}P}$
(16.23b)

The fugacity coefficient takes a value of unity when the substance behaves like an ideal gas. Therefore, the fugacity coefficient is also regarded as a measure of non-ideality; the closer the value of the fugacity coefficient is to unity, the closer we are to the ideal state.

Fugacity turns out to be an auxiliary function to chemical potential. Even though the concept of thermodynamic equilibrium which we discussed in the previous section is given in terms of chemical potentials, above definitions allow us to restate the same principle in terms of fugacity. To do this, previous expressions can be integrated for the change of state from liquid to vapor at saturation conditions to obtain:

(16.24a)
(16.24b)

For equilibrium, ${\mu }_{i}^{\left(l\right)}={\mu }_{i}^{\left(v\right)}$ ,hence,

$In\left(\frac{{f}_{i}^{\left(v\right)}}{{f}_{i}^{\left(l\right)}}\right)=0$
(16.24c)

Therefore:

(16.25)

For equilibrium, fugacities must be the same as well! This is, for a system to be in equilibrium, both the fugacity and the chemical potential of each component in each of the phases must be equal. Conditions (16.14) and (16.25) are equivalent. Once one of them is satisfied, the other is satisfied immediately. Using to describe equilibrium is a matter of choice, but generally the fugacity approach is preferred.

# Expressions for Fugacity Calculation

It is clear that, if we want to take advantage of the fugacity criteria to perform equilibrium calculations, we need to have a means of calculating it. Let us develop a general expression for fugacity calculations. Let us begin with the definition of fugacity in terms of chemical potential for a pure component shown in (16.21a):

(16.26)

The Maxwell’s Relationships presented in equation (15.27c) is written for a pure component system as:

${\left(\frac{\partial \mu }{\partial P}\right)}_{r}=\overline{V}=\stackrel{˜}{v}$
(16.27)

Consequently,

(16.28)

Substituting (16.28) into (16.26),

(16.29)

Introducing the concept of fugacity coefficient given in equation (16.23a),

$\varphi =\frac{f}{P}$
(16.23a)

(16.30)

We end up with:

(16.31a)

or equivalently,

(16.31b)

Integrating expression (16.31b),

(16.32)

It is convenient to define the lower limit of integration as the ideal state, for which the values of fugacity coefficient, volume, and compressibility factor are known.

At the ideal state, in the limit $P->0$ ,

${\varphi }^{*}->1\therefore \mathrm{ln}{\varphi }^{*}->0$
(16.33)

Substituting into (16.32),

E
(16.34)

Equation (16.34) is the expression of fugacity coefficient as a function of pressure, temperature, and volume. Notice that this expression can be readily rewritten in terms of compressibility factor:

(16.35)

Let us also derive the expression for the fugacity coefficient for a component in a multicomponent mixture. Following a pattern similar to that which we have presented, beginning with the definition of fugacity for a component in terms of chemical potential:

(16.36)

This time, it is more convenient to use the Maxwell’s Relationships presented in equation (15.27d):

${\left(\frac{\partial {\mu }_{i}}{\partial V}\right)}_{T,n}=-{\left(\frac{\partial P}{\partial {n}_{i}}\right)}_{T,V,{n}_{i\ne 1}}$
(16.37)

After you introduce the definitions of fugacity coefficient and compressibility factor:

${\varphi }_{i}=\frac{{f}_{i}}{{y}_{i}P}$
(16.38a)

$P=\frac{ZnRT}{V}$
(16.38b)

and recalling that our lower limit of integration is the ideal state, for which, at the limit $P->0$ :

${V}^{*}->\infty$
(16.39c)

(16.39a)

(16.39b)

it can be proven that:

(16.40)

The multi-component mixture counterpart of equation (16.35) becomes:

(16.41a)

where:

${\overline{Z}}_{i}={\left(\frac{\partial Z}{\partial {n}_{i}}\right)}_{P,T,{n}_{i\ne 1}}=\frac{P}{RT}{\left(\frac{\partial V}{\partial {n}_{i}}\right)}_{P,T,{n}_{i\ne 1}}=\frac{P\overline{{V}_{i}}}{RT}$
(16.41b)

Equations (16.34), (16.35), (16.40), and (16.41) are very important for us. Basically, they show that fugacity, or the fugacity coefficient, is a function of pressure, temperature and volume:

$f=f\left(P,V,T\right)$

This tells us that if we are able to come up with a PVT relationship for the volumetric behavior of a substance, we can calculate its fugacity by solving such expressions. It is becoming clear why we have studied equations of state — they are just what we need right now: PVT relationships for various substances. Once we have chosen the equation of state that we want to work with, we can calculate the fugacity of each component in the mixture by applying the above expression. Now that we know how to calculate fugacity, we are ready to apply the criteria for equilibrium that we just studied! That is the goal of the next module.

# Cubic EOS Fugacity Expressions

Expressions (16.34) and (16.40) are particularly suitable for the calculation of fugacity with P-explicit equations of state, which cubic equations of state are. One can take every cubic EOS we have presented and proceed with the integration, coming up with the expression for fugacity for that particular equation of state. We will spare the reader these derivations. The fugacity expressions for the cubic EOS of most interest to us (SRK and PR EOS) are presented below:

(16.42)

(16.43)

### PR EOS

(16.44)

#### Mixtures

(16.45)

The expressions for are the same as given before in the previous modules. (AA)i and (BB)i are calculated as:

${\left(AA\right)}_{i}=\frac{2}{{\left(a\alpha \right)}_{m}}\left[\sum _{j}^{{n}_{c}}{\left(a\alpha \right)}_{ij}\right]$
(16.46a)

${\left(BB\right)}_{i}=\frac{{b}_{i}}{{b}_{m}}$
(16.47b)

A and B parameters and the Z-factor of each phase are needed in order to calculate the corresponding fugacity coefficients. Now that we know how to calculate fugacity via EOS, and how this concept can be applied for equilibrium calculations (Section 16.2), we are ready to resume our discussion on Vapor-Liquid Equilibrium as we left it in Module 13. In our next module, we will concentrate on Vapor-Liquid Equilibrium via EOS.

# Action Item

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

• Your answer 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.

#### Problem Set

1. Provide a short lay person definition of Chemical Potential and Fugacity. Explain how these variables can be used to describe phase equilibrium.

# Vapor-Liquid Equilibrium via EOS

### Module Goals

Module Goal: To use cubic equations of state for the description of phase equilibrium.

Module Objective: To integrate the knowledge of qualitative phase behavior and the framework for quantification to actually build a phase equilibrium model.

# Introduction

Gas and liquid co-existence is common to all petroleum and natural gas applications, as we have mentioned already. Additionally, in our discussions presented in Modules 12 and 13, we tried to mathematically model the problem of vapor-liquid co-existence in equilibrium. As we recall, at that time we concluded that we were lacking tools necessary to complete the definition of the problem. Then we went on and concentrated on building those thermodynamic tools which we needed. Now we are ready to couple the concepts reviewed in Modules 14 through 16 (“Thermodynamic Tools”) with those reviewed in Modules 12 and 13. In order for us to do that, we will need to apply our knowledge of equations of state (studied in Modules 6 through 11).

Once the mathematics of vapor-liquid equilibrium is understood, we will be able to solve for a variety of applications. They will be discussed in the last modules of this course. One of those applications will be to generate any sort of phase diagrams, like the ones that we studied in our earlier discussions.

As you see, everything that we have studied in this course is completely interlinked. That’s one of the main ideas that I want you to glean from these discussions.

# Equilibrium and Equilibrium Ratios (Ki)

Consider a liquid-vapor in equilibrium. As we have discussed previously, a condition for equilibrium is that the chemical potential of each component in both phases are equal, thus:

${\mu }_{i}^{v}={\mu }_{i}^{L}$
(17.1)

We showed that this is equivalent to:

${f}_{i}^{V}={f}_{i}^{L}$
(17.2)

This is, for a system to be in equilibrium, the fugacity of each component in each of the phases must be equal as well. The fugacity of a component in a mixture can be expressed in terms of the fugacity coefficient. Therefore, the fugacity of a component in either phase can be written as:

${f}_{i}^{V}={y}_{i}{\varphi }_{i}^{V}P$
(17.3a)

${f}_{i}^{L}={y}_{i}{\varphi }_{i}^{L}P$
(17.3b)

Introducing (17.3) into (17.2),

${y}_{i}{\varphi }_{i}^{V}P={x}_{i}{\varphi }_{i}^{L}P$
(17.4)

This equilibrium condition can be written in terms of the equilibrium ratio ${K}_{i}={y}_{i}l{x}_{i}$ , to get:

${K}_{i}=\frac{{y}_{i}}{{x}_{i}}=\frac{{\varphi }_{i}^{L}}{{\varphi }_{i}^{V}}$
(17.5)

Do you recall the problem at the end of Module 13? At that point we needed a more reliable way to calculate the equilibrium ratios that showed up in the Rachford-Rice objective function. We demonstrated that once we know all values of Ki’s, the problem of vapor-liquid equilibrium is reduced to solving the Rachford-Rice objective function, using the Newton-Raphson Procedure.

We can now calculate equilibrium ratios, using (17.5), in terms of fugacity coefficients. We also know that we have an analytic expression for the calculation of fugacity coefficients via EOS — this was shown in the last section of the previous module. This is why we call this module “Vapor Liquid Equilibrium via EOS.”

Is this the end to our problems? Not quite. Take a look at the expression for fugacity coefficients in mixtures both for SRK EOS and PR EOS. It is clear that they are functions of the pressure, temperature, and composition of the phases:

${\varphi }_{i}^{L}={\varphi }_{i}^{L}\left(P,T,{x}_{i}\right)$
(17.6a)

${\varphi }_{i}^{V}={\varphi }_{i}^{V}\left(P,T,{y}_{i}\right)$
(17.6b)

Do we know the composition of the phases “xi”, “yi” in advance? In a typical flash problem, we are given pressure, temperature and overall composition (zi). What do we want to know? How much gas, how much liquid, and the compositions of the phases: . So, we do not know those compositions in advance; therefore, as it stands, we cannot calculate (17.6) or (17.5). Thus far, it seems that the flash problem is unsolvable.

If we are bold enough, we could try to overcome this problem by “guessing” those compositions, and proceed by solving (17.6) and (17.5). With this “rough” estimate for Ki’s, we could solve for “${\alpha }_{g}$ ” with the procedure outlined in Module 13 (“Objective Function and Newton-Raphson Procedure”). Once “${\alpha }_{g}$” is known, we could back calculate the compositions of the phases using equations (12.7) and (12.11). If we were correct, those compositions would match each other (the “guessed” ones with respect to the “back-calculated”). More than likely, this would not happen, and we would have to make a new “guess.” This is, fundamentally, an iterative procedure. Although this is not what we do, it does illustrate that this problem is solvable by implementing the appropriate iterative scheme.

In equations (17.4) and (17.5), the fugacity of the liquid and vapor phases were computed in terms of the fugacity coefficient. Hence, this method of expressing the equilibrium criteria is known as the dual-fugacity coefficient method. For the sake of completeness, it is necessary to indicate that the fugacity of a component in a mixture can also be expressed in terms of a thermodynamic concept called the activity coefficient. While the fugacity coefficient is seen as a measure of the deviation of behavior with respect to the ideal gas model, the activity coefficient measures the deviation of behavior with respect to the ideal liquid model. This approach is called the dual-activity coefficient method, in which both liquid and vapor phase fugacities are expressed in terms of the activity coefficient and substituted into the equilibrium criteria in (17.2). A mixed activity coefficient-fugacity coefficient method can be also devised by expressing the liquid phase fugacities in terms of activity coefficients and the vapor phase fugacities in terms of fugacity coefficients. Each of the aforementioned methods for the calculation of phase equilibria has its advantages and disadvantages. The dual-fugacity-coefficient method is simpler both conceptually and computationally, but if the equation of state does not predict liquid and vapor densities well, the results may be inaccurate. The activity coefficient method can be more accurate, but it is more complicated to implement. For the rest of the discussion, the dual-fugacity coefficient approach will be used.

# Solution Algorithms for VLE Problems

We have essentially reduced the typical VLE problem to that of solving a system of non-linear algebraic equations. The objective function is the nucleus of the VLE calculation. Rather than solve the material balance equations in their elemental form, the use of an objective function will allow for a more stable iteration scheme. We have seen that three objective functions are available for this calculation:

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

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

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

Expression (17.7c) is the Rachford-Rice Objective Function; we have proven that its numerical solution is cleaner than that of equations (17.7a) and (17.7b). We also discussed that the key problem in solving for “${\alpha }_{g}$” is the need for values for the equilibrium ratios (Ki’s). With the aid of equations (17.5) and (17.6) these problems can be solved simultaneously by iteration techniques.

In the problem we have posed, we have (3n+1) unknowns:

And, ${\alpha }_{g}$

And we also have (3n+1) equations:

${K}_{i}=\frac{{\varphi }_{i}^{L}\left(P,T,{x}_{i}\right)}{{\varphi }_{i}^{V}\left(P,T,{y}_{i}\right)}i=1,2,...n$
(17.8a)

(17.8b)

(17.8c)

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

where P, T, and zi are the values that we know. See how this is only a restatement of the original set of equations that contain the equilibrium conditions, material balances, and molar fraction constraints:

(17.9a)

(17.9b)

(17.9c)

(17.9d)

Since we have the same number of equations as unknowns, this system is determined and has a unique solution. To solve this system of simultaneous non-linear equations, there are two categories of solution techniques:

• Newton-type methods,
• Substitution-type methods.

Newton-type methods for more than one unknown require finding (n x n) elements of a Jacobian matrix at every iteration step. This could be computationally expensive. Additionally, the Newton Raphson method requires a very good initial guess (i.e., close to the actual solution) to guarantee convergence to the true values. This is not always possible, especially at the start of the procedure.

The most popular method, and the easiest to implement, is the substitution type. However, substitution type methods can be quite slow for some conditions of interest. In these sort of cases, we either switch to a Newton-Raphson solver (that has a much better rate of convergence) or implement an acceleration procedure to the substitution method itself.

# Successive Substitution Method (SSM)

In a substitution-type method, we start with initial guesses for all of the unknowns and loop around the equations to obtain “better” approximations for each of them. We test the goodness of the solution at every time step by comparing the new, better approximation to the previous guess. If the correction is small under certain convergence criteria, the procedure is stopped and we use the results from the last iteration as the final answer.

As we discussed in the previous section, reliable values for the equilibrium constants (Ki’s) must be obtained before we can solve the Rachford-Rice Objective Function. Generally, first estimates for equilibrium constants are calculated by using Wilson’s empirical equation (equation 13.15). However, Wilson’s correlation yields only approximate values for equilibrium ratios. Here, we apply thermodynamic equilibrium considerations in order to obtain more reliable predictions for Ki’s.

Let us recall that for a system to be in equilibrium, any net transfer (of heat, momentum, or mass) must be zero. For this to be true, all potentials (temperature, pressure, and chemical potential) must be the same in all of the phases. Provided that the temperature and pressure of both phases are the same, a zero net transfer for all components in the mixture results when all chemical potentials are the same. A restatement of this: all fugacities for all components in each phase are equal. Since fugacity is a measure of the potential for transfer of a component between two phases, equal fugacities of a component in both phases results in zero net transfer. SSM (Successive Substitution Method) takes advantage of this fact. Equation (17.5) showed that the fugacity criterion for equilibrium allows writing the equilibrium ratios Ki as a function of fugacity coefficients as follows:

${K}_{i}=\frac{{\varphi }_{i}^{L}}{{\varphi }_{i}^{V}}$
(17.10)

Equation (17.10) presupposes equality of fugacity (${f}_{i}^{V}={f}_{i}^{L}$), but during an iterative procedure, fugacities may not be equal (${f}_{i}^{V}\ne {f}_{i}^{L}$) until convergence conditions have been attained. Therefore, if fugacities are not equal (convergence has not been achieved), equation (17.10) becomes:

${K}_{i}=\frac{{\varphi }_{i}^{L}}{{\varphi }_{i}^{V}}=\frac{{f}_{i}^{L}/\left({x}_{i}P\right)}{{f}_{i}^{V}/\left({y}_{i}P\right)}=\frac{{y}_{i}}{{x}_{i}}\left(\frac{{f}_{i}^{L}}{{f}_{i}^{V}}\right)$
(17.11)

Using the above equation, a correction step can be formulated to improve the current values of Ki. The SSM-updating step is written as:

${K}_{i}^{n+l}={\left(\frac{{y}_{i}}{{x}_{i}}\right)}^{n}{\left(\frac{{f}_{li}}{{f}_{gi}}\right)}^{n}$
${K}_{i}^{n+l}={K}_{i}^{n}{\left(\frac{{f}_{li}}{{f}_{gi}}\right)}^{n}$
(17.12a)

SSM updates all previous equilibrium ratios (Ki) using the fugacities predicted by the equation of state. This iteration method requires an initial estimation of Ki-values — for which Wilson’s correlation is used. It can easily be concluded that the convergence criteria will be satisfied whenever the fugacity ratios of all the components in the system are close to unity. Such condition is achieved when the following inequality is satisfied:

${\sum _{i}^{n}\left(\frac{{f}_{li}}{{f}_{gi}}-1\right)}^{2}<{10}^{-14}$
(17.12b)

# The Accelerated Successive Substitution Method (ASSM)

When the system is close to the critical point and fugacities are strongly composition-dependent, a slowing-down of the convergence rate of the SSM (Successive Substitution Method) is to be expected. In an attempt to avoid slow convergence problems, some methods have been proposed. Among the most popular are the Minimum Variable Newton Raphson (MVNR) Method and the Accelerated and Stabilized Successive Substitution Method (ASSM).

The ASSM is basically an accelerated version of the SSM procedure, and thus follows a similar theory. Such procedure is implemented to accelerate the calculation of Ki-values, especially in the region close to critical point where the use of the SSM alone will not be efficient. The ASSM technique was presented by Rinses et al. (1981) and consists of the following steps:

1. Use the SSM technique to initiate the updating of the Ki-values the first time.
2. Check all following criteria at every step during iterations using the SSM:
$\frac{\sum _{i}^{nc}{\left(R{r}_{i}^{new}-1\right)}^{2}}{\sum _{i}^{nc}{\left(R{r}_{i}^{old}-1\right)}^{2}}>0.8$
(17.13a)

(17.13b)

${10}^{-5}<\sum _{l}^{nc}\left(R{r}_{i}^{new}-1{\right)}^{2}<{10}^{-3}$
(17.13c)

$0
(17.13d)

These criteria show that you have sufficient proximity to the conditions to ensure the efficiency of the method. Rri is the ratio of liquid fugacity to gas fugacity of the i-th component and ‘${\alpha }_{g}$’ is molar gas fraction of the two-phase system.

$R{r}_{i}=\frac{{f}_{i}^{L}}{{f}_{i}^{V}}$
(17.14)

3. If the system satisfies ALL above criteria, the iteration technique is then switched from the SSM to the ASSM. Otherwise, SSM is used for the update of the Ki-values. The following expressions are used to update Ki-values in ASSM:

${K}_{i}^{new}={K}_{i}^{old}R{r}_{i}^{{\lambda }_{i}}$

where ${\lambda }_{i}=\left[\left(R{r}_{i}^{old}-1\right)/\left(R{r}_{i}^{old}-R{r}_{i}^{new}\right)\right]$

In some cases, using a constant acceleration value of ${\lambda }_{i}=2$ is good enough.

4. Once all the criteria in step (2) are satisfied, skip step (2) for the subsequent iterations and use the ASSM technique to update Ki-values until convergence is attained, unless it does not give acceptable new estimates (as stated next).
5. When ASSM is used, it must always be tested to show that it leads to an improved solution (i.e., that it brings fugacity ratios closer to unity). If not, it must be rejected and switched back to SSM.

Even though we are outlining Risnes et al.’s version of the accelerated Successive Substitution Method, there are several other published algorithms whose main purpose have also been to accelerate the successive substitution method. Fussel and Yanosik (1978), Michelsen (1982), and Mehra et al. (1983) are examples of such attempts. Risnes et al. version is the easiest and most straightforward to implement, but it is subjected to limitations.

# The Stability Criteria

Interestingly enough, one of the most difficult aspects of making VLE calculations may not be the two-phase splitting calculation itself, but knowing whether or not a mixture will actually split into two (or even more) phases for a pressure and temperature condition.

A single-phase detection routine has to be simultaneously introduced at this stage to detect whether the system is in a true single-phase condition at the given pressure and temperature or whether it will actually split into two-phases. Several approaches may be used here: the Bring-Back technique outlined by Risnes et al. (1981), and Phase Stability Criteria introduced by Michelsen (1982), among others. Here we describe Michelsen’s stability test.

Michelsen (1982) suggested creating a second-phase inside any given mixture to verify whether such a system is stable or not. It is the same idea behind the Bring-Back procedure (Risnes et al., 1981), but this test additionally provides straightforward interpretation for the cases where trivial solutions are found
(Ki’s —> 1). The test must be performed in two parts, considering two possibilities: the second phase can be either vapor-like or liquid-like. The outline of the method is described below, following the approach presented by Whitson and Brule (2000).

1. Calculate the mixture fugacity (fzi) using overall composition zi.
2. Create a vapor-like second phase,
1. Use Wilson’s correlation to obtain initial Ki-values.
2. Calculate second-phase mole numbers, Yi:
${Y}_{i}={z}_{i}{K}_{i}$
(17.15)

3. Obtain the sum of the mole numbers,
${S}_{v}=\sum _{i}^{n}{Y}_{i}$
(17.16)

4. Normalize the second-phase mole numbers to get mole fractions:
${y}_{i}=\frac{{Y}_{i}}{{S}_{v}}$
(17.17)

5. Calculate the second-phase fugacity (fyi) using the corresponding EOS and the previous composition.
6. Calculate corrections for the K-values:
${R}_{i}=\frac{{f}_{xi}}{{f}_{yi}}\frac{1}{{S}_{v}}$
(17.18)

${K}_{i}^{\left(n+1\right)}={K}_{i}^{\left(n\right)}{R}_{i}$
(17.19)

7. Check if:
1. Convergence is achieved:
$\sum _{i}^{n}{\left({R}_{i}-1\right)}^{2}<1\cdot {10}^{-10}$
(17.20)

2. A trivial solution is approached:
(17.21)

If a trivial solution is approached, stop the procedure.

If convergence has not been attained, use the new K-values and go back to step (b).

3. Create a liquid-like second phase,

Follow the previous steps by replacing equations (17.15), (17.16), (17.17), and (17.18) by (17.22), (17.23), (17.24), and (17.25) respectively.

${Y}_{i}={z}_{i}/{K}_{i}$
(17.22)

${S}_{L}=\sum _{i}^{n}{Y}_{i}$
(17.23)

${x}_{i}=\frac{{Y}_{i}}{{S}_{L}}$
(17.24)

${R}_{i}=\frac{{f}_{xi}}{{f}_{zi}}{S}_{L}$
(17.25)

The interpretation of the results of this method follows:

• The mixture is stable (single-phase condition prevails) if:
• Both tests yield S < 1 (SL < 1 and SV < 1),
• Or both tests converge to trivial solution,
• Or one test converges to trivial solution and the other gives S < 1.
• Only one test indicating S > 1 is sufficient to determine that the mixture is unstable and that the two-phase condition prevails. The same conclusion is made if both tests give S > 1, or if one of the tests converges to the trivial solution and the other gives S > 1.

# Action Item

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

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

#### Problem Set

1. Provide a comprehensive flowchart of how to build a flash calculation model for a multi-component mixture. Assume that the data you have available is pressure, temperature, fluid composition, and all the properties of the pure components making up the mixture. This flowchart should be clear enough so that anyone could code it and generate results!

# Properties of Natural Gas and Condensates (I)

### Module Goals

Module Goal: To highlight the important properties used to characterize natural gas and condensate systems.

Module Objective: To present the most popular models for estimating properties of natural gas and condensate systems.

# Introduction

Fluid properties are used to characterize the condition of a fluid at a given state. A reliable estimation and description of the properties of hydrocarbon mixtures is fundamental in petroleum and natural engineering analysis and design. Fluid properties are not independent, just as pressure, temperature, and volume are not independent of each other. Equations of State provide the means for the estimation of the P-V-T relationship, and from them many other thermodynamic properties can be derived. Compositions are usually required for the calculation of the properties of each phase. For a VLE system, using the tools we have discussed in the previous lectures, we are ready to predict some important properties of both the liquid (condensate) and vapor (natural gas) phases — by means of the known values of composition of both phases. Some of the most relevant are discussed next.

# Molecular Weight

The molecular weight (MW) of each of the phases in a VLE system is calculated as a function of the molecular weight of the individual components (MWi), provided that both the composition of the gas (yi) and the liquid (xi) are known:

$M{W}_{g}=\sum _{i=1}^{n}{y}_{i}M{W}_{i}$
(18.1a)

$M{W}_{l}=\sum _{i=1}^{n}{x}_{i}M{W}_{i}$
(18.1b)

# Density

Density is customarily defined as the amount of mass contained in a unit volume of fluid. Density is the single-most important property of a fluid, once we realize that most other properties can be obtained or related to density. Both specific volume and density — which are inversely proportionally related to each other — tell us the story of how far apart the molecules in a fluid are from each other. For liquids, density is high — which translates to a very high molecular concentration and short intermolecular distances. For gases, density is low — which translates to low molecular concentrations large intermolecular distances.

The question then is: Given this, how can we obtain this all-important property called density? This takes us back to Equations of State (EOS). Since very early times, there have been correlations for the estimation of density of the liquids (oil, condensates) and gases/vapors (dry gases, wet gases). In modern times, equations of state (EOS) are a natural way of obtaining densities. The density of the fluid ‘f’ is calculated using its compressibility factor (Zf) as predicted by an appropriate equation of state. From the real gas law, the density can be expressed as:

${\rho }_{f}=\frac{P}{RT}\left(\frac{M{W}_{f}}{{Z}_{f}}\right)$
(18.2)

where: MWf is the molecular weight of fluid ‘f’. Expression (18.2) is used for both the gas and liquid density. In either case, the proper value for MWf (either MWg or MWl) and Zf (either Zg or Zo) has to be used. This takes us back to the discussion of equations of state. From Equation (18.2) it is clear that all that we need is the Z-factor.

The all-important parameter to calculate density is the Z-factor, both for the liquid and vapor phases. The relation between liquid behavior and Z-factor is not obvious, because Z-factor has been traditionally defined for gases. However, we can get “Z” for liquids. “Z” is, indeed, a measure of departure from the ideal gas behavior. Fair enough, for defining “Z” for liquids, we still measure the departure of liquid behavior from ideal gas behavior. A “liquid state” is a tremendous departure from ideal-gas conditions, and as such, “Z” for a liquid is always very far from unity. Typical values of “Z” for liquids are small.

Equations of State have proven very reliable for the estimation of vapor densities, but they do not do as good a job for liquid densities. There is actually a debate among different authors about the reliability of Z-factor estimations for liquids using EOS. In fact, people still believe the EOS are not reliable for liquid density predictions and that we should use correlations instead. However, Peng-Robinson EOS provides fair estimates for vapor and liquid densities as long as we are dealing with natural gas and condensate systems.

Empirical correlations for Z-factor for natural gases were developed before the advent of digital computers. Although their use is in decline, they can still be used for fast estimates of the Z-factor. The most popular of such correlations include those of Hall-Yarboroug and Dranchuk-Abou-Kassem.

Chart look-up is another means of determining Z-factor of natural gas mixtures. These methods are invariably based on some type of corresponding states development. According to the theory of corresponding states, substances at corresponding states will exhibit the same behavior (and hence the same Z-factor). The chart of Standing and Katz is the most commonly used Z-factor chart for natural gas mixtures.

Methods of direct calculation using corresponding states have also been developed, ranging from correlations of chart values to sophisticated equation sets based on theoretical developments.

However, the use of equations of state to determine Z-factors has grown in popularity as computing capabilities have improved. Equations of state represent the most complex method of calculating Z-factor, but also the most accurate. A variety of equations of state have been developed to describe gas mixtures, ranging from the ideal EOS (which yields only one root for the vapor and poor predictions at high pressures and low temperatures), cubic EOS (which yields up to three roots, including one for the liquid phase), and more advanced EOS such as BWR and AGA8.

#### References:

Hall K., and Yarborough, L. (1973), “A New Equation of State for Z-factor Calculations”, Oil and Gas Journal, June 1973, pp. 82-92.

Dranchuk, P. and Abou-Kassem, J. (1975), “Calculation of Z-factors for Natural Gases Using Equations-of-State”, JCPT, July-September 1975, p. 34-36.

Standing, M. and Katz, D. (1942), “Density of Natural Gases”, Trans. AIME, v. 146, pp. 140-149.

# Specific Gravity

Specific gravity is defined as the ratio of fluid density to the density of a reference substance, both defined at the same pressure and temperature. These densities are usually defined at standard conditions (14.7 psia and 60°F). For a condensate, oil or a liquid, the reference substance is water:

${\gamma }_{o}=\frac{{\left({\rho }_{0}\right)}_{sc}}{{\left({\rho }_{w}\right)}_{sc}}$
(18.3)

The value of water density at standard conditions is 62.4 lbm/ft3 approximately. For a natural gas, or any other gas for this matter, the reference substance is air:

${\gamma }_{g}=\frac{{\left({\rho }_{g}\right)}_{sc}}{{\left({\rho }_{air}\right)}_{sc}}$
(18.3a)

Or, equivalently, substituting Equation (18.2) evaluated at standard conditions (${Z}_{sc}\approx 1$ for most gases),

${\gamma }_{g}=\frac{M{W}_{g}}{M{W}_{air}}$
(18.3b)

where the value of the molecular weight for air is MWair = 28.96 lbm/lbmol. Specific gravity is nondimensional because both numerator and denominator have the same units.

# API

Petroleum and Natural Gas Engineers also use another gravity term which is called API gravity. It is used for liquids (e.g., condensates) and is defined as:

$°API=\frac{141.5}{{\gamma }_{o}|{}_{sc}}-131.5$
(18.4)

By definition (see Equation 18.3), the specific gravity of water is unity. Therefore, water has an API gravity of 10. The API gravity of 10 is associated with very heavy, almost asphaltic, oils. Light crude oils have an API greater than or equal to 45°. Condensate gravities range between 50° and 70° API. Liquid condensates are normally light in color.

# Volumetric Factors (Bo and Bg)

Due to the dramatically different conditions prevailing at the reservoir when compared to the conditions at the surface, we do not expect that 1 barrel of fluid at reservoir conditions could contain the same amount of matter as 1 barrel of fluid at surface conditions. Volumetric factors were introduced in petroleum and natural gas calculations in order to readily relate the volume of fluids that are obtained at the surface (stock tank) to the volume that the fluid actually occupied when it was compressed in the reservoir.

For example, the volume that a live oil occupies at the reservoir is more than the volume of oil that leaves the stock tank at the surface. This may be counter-intuitive. However, this is a result of the evolution of gas from oil as pressure decreases from reservoir pressure to surface pressure. If an oil had no gas in solution (i.e., a dead oil), the volume that it would occupy at reservoir conditions is less than the volume that it occupies at the surface. In this case, only liquid compressibility plays a role in the change of volume.

The formation volume factor of a natural gas (Bg) relates the volume of 1 lbmol of gas at reservoir conditions to the volume of the same lbmol of gas at standard conditions, as follows:

(18.5)

Those volumes are, evidently, the specific molar volumes of the gas at the given conditions. The reciprocal of the specific molar volume is the molar density, and thus, Equation (18.5) could be written:

${B}_{g}=\frac{{\stackrel{˜}{v}}_{g}/{}_{res}}{{\stackrel{˜}{v}}_{g}/{}_{sc}}=\frac{{\overline{\rho }}_{g}/{}_{sc}}{{\overline{\rho }}_{g}/{}_{res}}=\frac{\left({\rho }_{g}/M{W}_{g}\right){|}_{sc}}{\left({\rho }_{g}/M{W}_{g}\right){|}_{res}}$
(18.6)

Introducing the definition for densities in terms of compressibility factor,

${B}_{g}=\frac{\frac{{P}_{sc}}{R{T}_{sc}{Z}_{sc}}}{\frac{P}{RTZ}}$
(18.7)

Therefore, recalling that ${Z}_{sc}\approx 1$ ,

${B}_{g}=\frac{{P}_{sc}}{{T}_{sc}}\frac{ZT}{P}=0.005035{\frac{ZT}{P}}_{\left[RCF/SCF\right]}$
(18.8)

Gas formation volume factors can be also expressed in terms of [RB/SCF]. In such a case, 1 RB = 5.615 RCF and we write:

${B}_{g}=0.005035{\frac{ZT}{P}}_{\left[RCF/SCF\right]}$
(18.9)

The formation volume factor of an oil or condensate (Bo) relates the volume of 1 lbmol of liquid at reservoir conditions to the volume of that liquid once it has gone through the surface separation facility.

(18.10)

The total volume occupied by 1 lbmol of liquid at reservoir conditions (Vo)res can be calculated through the compressibility factor of that liquid, as follows:

(18.11)

Upon separation, some gas is going to be taken out of the liquid stream feeding the surface facility. Let us call “nst” the moles of liquid leaving the stock tank per mole of feed entering the separation facility. The volume that 1 lbmol of reservoir liquid is going to occupy after going through the separation facility is given by:

(18.12)

Here we assume that the last stage of separation, the stock tank, operates at standard conditions. Introducing Equations (18.12) and (18.11) into (18.10), we end up with:

${B}_{o}=\frac{{\left(\frac{n{Z}_{o}RT}{P}\right)}_{res}}{{\left(\frac{{n}_{st}{Z}_{0}RT}{P}\right)}_{sc}}$
(18.13)

or,

${B}_{o}=\frac{1}{{n}_{st}}\frac{{\left({Z}_{o}\right)}_{res}}{{\left({Z}_{o}\right)}_{sc}}\frac{T}{P}{\frac{{P}_{sc}}{{T}_{sc}}}_{\left[RB/STB\right]}$
(18.14)

Please notice that (Zo)sc — unlike Zsc for a gas — is never equal to one. Oil formation volume factor can be also seen as the volume of reservoir fluid required to produce one barrel of oil in the stock tank.

# Isothermal Compressibilities

The isothermal compressibility of a fluid is defined as follows:

${c}_{f}=-\frac{1}{V}{\left(\frac{\partial V}{\partial \rho }\right)}_{T}$
(18.15)

This expression can be also given in term of fluid density, as follows:

${c}_{f}=-\frac{1}{\rho }{\left(\frac{\partial \rho }{\partial P}\right)}_{T}$
(18.16)

For liquids, the value of isothermal compressibility is very small because a unitary change in pressure causes a very small change in volume for a liquid. In fact, for slightly compressible liquid, the value of compressibility (co) is usually assumed independent of pressure. Therefore, for small ranges of pressure across which co is nearly constant, Equation (18.16) can be integrated to get:

${c}_{o}\left(p-{p}_{b}\right)=\mathrm{ln}\left(\frac{{\rho }_{o}}{{\rho }_{ob}}\right)$
(18.17)

In such a case, the following expression can be derived to relate two different liquid densities at two different pressures (p, pb):

${\rho }_{o}={\rho }_{ob}\left[1+{c}_{o}\left(p-{p}_{b}\right)\right]$
(18.18)

The Vasquez-Beggs correlation is the most commonly used relationship for co.

For natural gases, isothermal compressibility varies significantly with pressure. By introducing the real gas law into Equation (18.16), it is easy to prove that, for gases:

${c}_{g}=\frac{1}{P}-\frac{1}{Z}{\left(\frac{\partial Z}{\partial P}\right)}_{r}$
(18.19)

Note that for an ideal gas, cg is just the reciprocal of the pressure. “cg” can be readily calculated by graphical means (chart of Z versus P) or by introducing an equation of state into Equation (18.19).

# Surface Tension

Surface tension is a measure of the surface free energy of liquids, i.e., the extent of energy stored at the surface of liquids. Although it is also known as interfacial force or interfacial tension, the name surface tension is usually used in systems where the liquid is in contact with gas.

Qualitatively, it is described as the force acting on the surface of a liquid that tends to minimize the area of its surface, resulting in liquids forming droplets with spherical shape, for instance. Quantitatively, since its dimension is of force over length (lbf/ft in English units), it is expressed as the force (in lbf) required to break a film of 1 ft of length. Equivalently, it may be restated as being the amount of surface energy (in lbf-ft) per square feet.

Katz et al. (1959) presented the Macleod-Sudgen equation for surface tension ($\sigma$) calculations in dynes/cm for hydrocarbon mixtures:

${\sigma }^{1/4}=\sum _{i=1}^{n}Pc{h}_{i}\left[\frac{{\rho }_{i}}{62.4\left(M{W}_{l}\right)}{x}_{i}-\frac{{\rho }_{g}}{62.4\left(M{W}_{g}\right)}{y}_{i}\right]$
(18.20)

where:

Pchi = Parachor of component “i”,

xi = mole fraction of component “i” in the liquid phase,

yi = mole fraction of component “i” in the gas phase.

In order to express surface tension in field units (lbf/ft), multiply the surface tension in (dynes/cm) by 6.852177x10-3. The parachor is a temperature independent parameter that is calculated experimentally. Parachors for pure substances have been presented by Weinaug and Katz (1943) and are listed in Table 18.1.

Table 18.1. Parachors for pure substances (Weinaug and Katz, 1943)
Component Parachor
CO2 78.0
N2 41.0
C1 77.0
C2 108.0
C3 150.3
iC4 181.5
nC4 189.9
iC5 225.0
nC5 231.5
nC6 271.0
nC7 312.5
nC8 351.5

Weinaug and Katz (1943) also presented the following empirical relationship for the parachor of hydrocarbons in terms of their molecular weight:

$Pc{h}_{i}=-4.6148734+2.558855M{W}_{i}+3.404065\cdot {10}^{-4}M{W}_{i}^{2}+\frac{3.767396\cdot {10}^{3}}{M{W}_{i}}$
(18.21)

• This correlation may be used for pseudo-components or for those hydrocarbons not shown in Table 18.1.

# Action Item

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

• Your answer 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.

#### Problem Set

1. Which of all the properties we have studied so far would you use to compare and distinguish among wet natural gas, dry natural gas, and gas-condensate systems? Explain why.

# Properties of Natural Gas and Condensates (II)

### Module Goals

Module Goal: To highlight the important properties used to characterize natural gas and condensate systems.

Module Objective: To present the most popular models for estimating properties of natural gas and condensate systems.

# Heat Capacities

The constant volume heat capacity is defined by:

${C}_{v}{\left(\frac{\partial U}{\partial T}\right)}_{V}$
(19.1)

To see the physical significance of the constant volume heat capacity, let us consider a 1 lbmol of gas within a rigid-wall (constant volume) container. Heat is added to the system through the walls of the container and the gas temperature rises. It is evident that the temperature rise ($\Delta T$ ) is proportional to the amount of heat added,

$Q\propto \Delta T$
(19.2)

Introducing a constant of proportionality “cv”,

$Q={c}_{v}\Delta T$
(19.3)

In our experiment, no work was done because the boundaries (walls) of the system remained unchanged. Applying the first law of thermodynamics to this closed system, we have:

$\Delta U={c}_{v}\Delta T$
(19.4)

Therefore, for infinitesimal changes,

${C}_{v}={\left(\frac{\partial U}{\partial T}\right)}_{V}$
(19.5)

As we have seen, constant volume heat capacity is the amount of heat required to raise the temperature of a gas by one degree while retaining its volume.

Let us now consider the same 1 lbmol of gas confined in a piston-cylinder equipment (i.e., a system with non-rigid walls or boundaries). When heat is added to the system, the gas temperature rises and the gas expands so that the pressure in the system remains the same at any time. The piston displaces a volume $\Delta$ V and the gas increases its temperature in $\Delta$T degrees. Again, the temperature rise ($\Delta$T) is proportional to the amount of heat added, and the new constant of proportionality we use here is “cp”,

$Q={c}_{P}\Delta T$
(19.6)

This time, some work was done because the boundaries (walls) of the system changed from their original position. Applying the first law of thermodynamics to this closed system, we have that:

$\Delta U=Q-W$
(19.7)

If the pressure remained the same both inside and outside the container, the system made some work against the surroundings in the amount of $W={P}_{\Delta }V$. Introducing (19.7) into (19.6),

$\Delta U+P\Delta V={c}_{P}\Delta T$
(19.8)

The left hand side of this equation represents the definition of enthalpy change ($\Delta H$ ) for a constant-pressure process. Therefore:

$\Delta H={c}_{P}\Delta T$
(19.9)

Finally, for infinitesimal changes,

${c}_{P}={\left(\frac{\partial H}{\partial T}\right)}_{V}$
(19.10)

The function “cp” is called the constant pressure heat capacity. The constant pressure heat capacity is the amount of heat required to raise the temperature of a gas by one degree while retaining its pressure.

The units of both heat capacities are (Btu/lbmol-°F) and (cal/gr-°C). Their values are never equal to each other, not even for ideal gases. In fact, the ratio “cp/cv” of a gas is known as “k” — the heat capacity ratio — and it is never equal to unity. This ratio is frequently used in gas-dynamics studies.

$k=\frac{{c}_{p}}{{c}_{v}}$
(19.11)

Heat capacities can be calculated using equations of state. For instance, Peng and Robinson (1976) presented an expression for the departure enthalpy of a fluid mixture, shown below:

$\stackrel{¨}{H}=H-{H}^{·}=RT\left(Z-1\right)+\frac{T\frac{d{\left(a\alpha \right)}_{m}}{dT}-{\left(a\alpha \right)}_{m}}{2\sqrt{2{b}_{m}}}\mathrm{ln}\left(\frac{Z+\left(\sqrt{2}+1\right)B}{Z-\left(\sqrt{2}-1\right)B}\right)$
(19.12)

The value of the enthalpy of the fluid (H) is obtained by adding up this enthalpy of departure (shown above) to the ideal gas enthalpy (H*). Ideal enthalpies are sole functions of temperature. For hydrocarbons, Passut and Danner (1972) developed correlations for ideal gas properties such as enthalpy, heat capacity and entropy as a function of temperature. Therefore, an analytical relationship for “cp” can be derived taking the derivative of (19.12), as shown below:

$\begin{array}{l}{C}_{P}={C}_{P}^{·}+R\left(T{\left(\frac{\partial Z}{\partial T}\right)}_{P}+Z-1\right)+\frac{T\frac{d{\left(a\alpha \right)}_{m}}{dT}-{\left(a\alpha \right)}_{m}}{2\sqrt{2{b}_{m}}}\\ \left[\frac{{\left(\frac{\partial Z}{\partial T}\right)}_{P}+2.414{\left(\frac{\partial B}{\partial T}\right)}_{P}}{Z+2.414B}-\frac{{\left(\frac{\partial Z}{\partial T}\right)}_{P}-0.414{\left(\frac{\partial B}{\partial T}\right)}_{P}}{Z-0.414B}\right]+\frac{T\frac{{d}^{2}{\left(a\alpha \right)}_{m}}{d{T}^{2}}}{2\sqrt{2{b}_{m}}}\mathrm{ln}\left(\frac{Z+2.414B}{Z-414B}\right)\end{array}$
(19.13)

where:

${C}_{P}^{·}={\left(\frac{\partial {H}^{·}}{\partial T}\right)}_{P}$ = ideal gas CP,

also found in the work of Passut and Danner (1972).

The second derivative of ${\left(a\alpha \right)}_{m}$ with respect to temperature can be calculated through the expression:

$\begin{array}{l}\frac{{d}^{2}}{d{T}^{2}}{\left(a\alpha \right)}_{m}=-\frac{0.45724{R}^{2}}{2\sqrt{T}}\underset{i}{{\sum }^{\text{​}}}\sum _{j}{c}_{i}{c}_{j}\left(1-{k}_{ij}\right)\\ \left[f\left({w}_{j}\right)\left(\frac{{\alpha }_{i}^{0.5}{T}_{ci}}{{P}_{ci}^{0.5}}\right){\left(\frac{{T}_{cj}}{{P}_{cj}}\right)}^{0.5}{\psi }_{i}+f\left({w}_{i}\right)\left(\frac{{\alpha }_{j}^{0.5}{T}_{cj}}{{P}_{cj}^{0.5}}\right){\left(\frac{{T}_{ci}}{{P}_{ci}}\right)}^{0.5}{\psi }_{j}\right]\end{array}$
(19.14a)

where,

${\psi }_{i}=-\frac{f\left({w}_{i}\right)}{2\sqrt{{T}_{ci}T{\alpha }_{i}}}-\frac{1}{2T}$
(19.14b)

For the evaluation of expression (19.13), the derivative of the compressibility factor with respect to temperature is also required. Using the cubic version of Peng-Robinson EOS, this derivative can be written as:

${\left(\frac{\partial Z}{\partial T}\right)}_{P}=-\left(\frac{{\left(\frac{\partial {\Omega }_{2}}{\partial T}\right)}_{P}{Z}^{2}+{\left(\frac{\partial {\Omega }_{3}}{\partial T}\right)}_{P}Z+{\left(\frac{\partial {\Omega }_{4}}{\partial T}\right)}_{P}}{3{Z}^{2}+2{\Omega }_{2}Z+{\Omega }_{3}}\right)$
(19.15)

where,

$\begin{array}{l}{\left(\frac{\partial {\Omega }_{2}}{\partial T}\right)}_{P}={\left(\frac{\partial B}{\partial T}\right)}_{P}\\ {\left(\frac{\partial {\Omega }_{3}}{\partial T}\right)}_{P}={\left(\frac{\partial A}{\partial T}\right)}_{P}-6B{\left(\frac{\partial B}{\partial T}\right)}_{P}-2{\left(\frac{\partial B}{\partial T}\right)}_{P}\\ {\left(\frac{\partial {\Omega }_{4}}{\partial T}\right)}_{P}=-\left[A{\left(\frac{\partial B}{\partial T}\right)}_{P}+B{\left(\frac{\partial A}{\partial T}\right)}_{P}-2B{\left(\frac{\partial B}{\partial T}\right)}_{P}-3{B}^{2}{\left(\frac{\partial B}{\partial T}\right)}_{P}\right]\\ {\left(\frac{\partial A}{\partial T}\right)}_{P}=\frac{A}{{\left(a\alpha \right)}_{m}}\frac{d{\left(a\alpha \right)}_{m}}{dT}-2\frac{A}{T}\\ {\left(\frac{\partial B}{\partial T}\right)}_{P}=-\frac{B}{T}\end{array}$

“cp” and “cv” values are thermodynamically related. It can be proven that this relationship is controlled by the P-V-T behavior of the substances through the relationship:

${c}_{p}-{c}_{v}=T{\left(\frac{\partial V}{\partial T}\right)}_{P}{\left(\frac{\partial P}{\partial T}\right)}_{V}$
(19.16)

For ideal gases, $PV=nRT$ and Equation (18.28) collapses to:

${c}_{p}^{·}-{c}_{v}^{·}=R$
(19.17)

# Joule-Thomson Coefficient

One remarkable difference between flow of condensate (or liquid) and natural gases through a pipeline is that of the effect of pressure drop on temperature changes along the pipeline. This is especially true when heat losses to the environment do not control these temperature variations. Natural gas pipelines usually cool with distance (effect commonly called ‘Joule–Thomson cooling’), while oil lines heat. The reason for such dissimilarity pertains to the different effect that pressure drop has on the entropy of a natural gas than on the entropy of an oil mixture. Katz (1972) and Katz and Lee (1990) presented a very enlightening discussion on this regard.

Whether or not a gas cools upon expansion or compression — that is, when subjected to pressure changes — depends on the value of its Joule–Thomson coefficient. This is not only important for natural gas pipeline flow, but also for the recovery of condensate from wet natural gases. In the cryogenic industry, turboexpanders are used to subject a wet gas to a sudden expansion (sharp pressure drop) in order to cool the gas stream beyond its dew point and recover the liquid dropout.

Thermodynamically, the Joule–Thomson coefficient is defined as the isenthalpic change in temperature in a fluid caused by a unitary pressure drop, as shown:

$\eta ={\left(\frac{\partial T}{\partial P}\right)}_{H}$
(19.18)

Using thermodynamic relationships, alternative expressions can be written. For example, using the cycling rule we may write:

${\left(\frac{\partial H}{\partial P}\right)}_{T}{\left(\frac{\partial P}{\partial T}\right)}_{H}{\left(\frac{\partial T}{\partial H}\right)}_{P}=-1$
(19.19)

or,

${\left(\frac{\partial H}{\partial P}\right)}_{T}=-{\left(\frac{\partial H}{\partial T}\right)}_{P}{\left(\frac{\partial T}{\partial P}\right)}_{H}$
(19.20)
${\left(\frac{\partial H}{\partial P}\right)}_{T}=-{c}_{P}\eta$
(19.21)

We have also seen that we can express enthalpy changes in terms of pressure, temperature and volume changes:

${\left(\frac{\partial H}{\partial P}\right)}_{T}=\left[\stackrel{˜}{v}-T{\left(\frac{\partial \stackrel{˜}{v}}{\partial T}\right)}_{P}\right]$
(19.22)

Additionally, the following identity can be derived:

$\eta =\frac{R{T}^{2}}{P{c}_{P}}{\left(\frac{\partial Z}{\partial T}\right)}_{P}$
(19.23)

All together, we have several ways of calculating the Joule–Thompson coefficient for a fluid, as shown next:

$\eta ={\left(\frac{\partial T}{\partial P}\right)}_{H}=\frac{1}{{c}_{P}}\left[T{\left(\frac{\partial \stackrel{˜}{v}}{\partial T}\right)}_{P}-\stackrel{˜}{v}\right]=-\frac{1}{{c}_{P}}{\left(\frac{\partial H}{\partial P}\right)}_{T}=\frac{R{T}^{2}}{P{c}_{P}}{\left(\frac{\partial Z}{\partial T}\right)}_{P}$
(19.24)

Once the constant pressure specific heat “cp” is calculated as discussed in the previous lecture, all the entries in the previous expression are known and the Joule–Thomson coefficient can be analytically calculated. An interesting observation from all above expressions for “$\eta$ ” is that the Joule–Thompson coefficient of an ideal gas is identically equal to zero. However, real fluids take positive or negative Joule–Thompson values.

# Viscosity

What other properties are we interested in? We are interested in flow properties. Whether you are interested in flow in pipes or in porous media, one of the most important transport properties is viscosity. Fluid viscosity is a measure of its internal resistance to flow. The most commonly used unit of viscosity is the centi-poise, which is related to other units as follows:

1 cp = 0.01 poise = 0.000672 lbm/ft-s = 0.001 Pa-s

Natural gas viscosity is usually expected to increase both with pressure and temperature. A number of methods have been developed to calculate gas viscosity. The method of Lee, Gonzalez and Eakin is a simple relation which gives quite accurate results for typical natural gas mixtures with low non-hydrocarbon content. Lee, Gonzalez and Eakin (1966) presented the following correlation for the calculation of the viscosity of a natural gas:

${\mu }_{g}=1\cdot {10}^{-4}{k}_{v}EXP\left({x}_{v}{\left(\frac{{\rho }_{g}}{62.4}\right)}^{{y}_{v}}\right)$
(19.25a)

where:

${k}_{v}=\frac{\left(9.4+0.02M{W}_{g}\right){T}^{1.5}}{209+19M{W}_{g}+T}$
(19.25b)

${y}_{v}=2.4-0.2{x}_{v}$
(19.25c)

${x}_{v}$
(19.25d)

In this expression, temperature is given in (°R), the density of the fluid (${\rho }_{g}$ ) in lbm/ft3 (calculated at the pressure and temperature of the system), and the resulting viscosity is expressed in centipoises (cp).

The most commonly used oil viscosity correlations are those of Beggs-Robinson and Vasquez-Beggs. Corrections must be applied for under-saturated systems and for systems where dissolved gas is present in the oil. However, in compositional simulation, where both gas and condensate compositions are known at every point of the reservoir, it is customary to calculate condensate viscosity using Lohrenz, Bray & Clark correlation. It this type of simulation, it is usual to calculate gas viscosities based on Lohrenz, Bray & Clark correlation as well. This serves the purpose of guaranteeing that the gas phase and condensate phase converge to the same value of viscosity as they approach near-critical conditions.

Lohrenz, Bray and Clark (1964) proposed an empirical correlation for the prediction of the viscosity of a liquid hydrocarbon mixture from its composition. Such expression, originally proposed by Jossi, Stiel and Thodos (1962) for the prediction of the viscosity of dense-gas mixtures, is given below:

$\mu ={\mu }^{·}+{\xi }_{m}^{-1}{\left(0.1023+0.023364{\rho }_{r}+0.058533{\rho }_{r}^{2}-0.040758{\rho }_{r}^{3}+0.0093724{\rho }_{r}^{4}\right)}^{4}-1\cdot {10}^{-4}$
(19.26)

where:
$\mu$ = fluid viscosity (cp),
${\mu }^{\ast }$ = viscosity at atmospheric pressure (cp),
${\xi }_{m}$ = mixture viscosity parameter (cp-1),
${\rho }_{r}$ = reduced liquid density (unitless),

Lohrenz et al. original paper presents a typographical error in Equation (19.26). Here it is written as originally proposed by Jossi, Stiel and Thodos (1962). All four parameters listed above have to be calculated as a function of critical properties in order to apply Equation (19.26). Lohrenz et al. original paper uses scientific units, here we present the equivalent equations in field (English) units.

For the viscosity of the mixture at atmospheric pressure (${\mu }^{\ast }$ ), Lohrenz et al. suggested using the following Herning & Zipperer equation:

${\mu }^{·}=\frac{\sum _{i}{z}_{i}{\mu }_{i}^{·}\sqrt{M{W}_{i}}}{\sum _{i}{z}_{i}\sqrt{M{W}_{i}}}$
(19.27)

where:
zj = mole composition of the i-th component in the mixture,
MWi = molecular weight of the i-th component (lbm/lbmol)
${\mu }_{i}^{·}$ = viscosity of the i-th component at low pressure (cp):
${\mu }_{i}^{·}=\frac{34\cdot {10}^{-5}{T}_{ri}^{0.94}}{{\xi }_{i}}$ [ if Tri ≤ 1.5 ]
${\mu }_{i}^{·}=\frac{17.78\cdot {10}^{-5}{\left(4.5{T}_{ri}-1.67\right)}^{0.625}}{{\xi }_{i}}$ [ if Tri > 1.5
where:
Tri = reduced temperature for the i-th component (T/Tci),
MWi = viscosity parameter of the i-th component, given by: ${\xi }_{i}=\frac{5.4402{T}_{ci}^{1/6}}{\sqrt{M{W}_{i}}{P}_{ci}^{2/3}}$
For the mixture viscosity parameter ($\xi m$ ), Lohrenz et al. applied an equivalent expression to that shown above but using pseudo-properties for the mixture:

$\xi m=\frac{5.4402{T}_{pc}^{1/6}}{\sqrt{M{W}_{l}}{P}_{pc}^{2/3}}$
(19.28)

where:

Tpc = pseudocritical temperature (oR),
Ppc = pseudocritical pressure (psia),
MWl = liquid mixture molecular weight (lbm/lbmol).

The reduced density of the liquid mixture (${\rho }_{r}$ ) is calculated as:

${\rho }_{r}=\frac{{\rho }_{l}}{{\rho }_{pc}}=\left(\frac{{\rho }_{l}}{M{W}_{l}}\right){V}_{pc}$
(19.29)

where:
${\rho }_{pc}$ = mixture pseudocritical density (lbm/ft3),
Vpc = mixture pseudocritical molar volume (ft3/lbmol),

All mixture pseudocritical properties are calculated using Kay’s mixing rule, as shown:

${T}_{pc}=\sum {z}_{i}{T}_{ci}$
(19.30a)

${P}_{pc}=\sum {z}_{i}{P}_{ci}$
(19.30b)

${V}_{pc}=\sum {z}_{i}{V}_{ci}$
(19.30c)

“zi” pertains to the fluid molar composition, Tci is given in oR, Pci in psia, and Vci in ft3/lbmol. When the critical volumes are known in a mass basis (ft3/lbm), each of them is to be multiplied by the corresponding molecular weight. In the case of lumped C7+ heavy fractions, Lorentz et al. (1969) presented a correlation for the estimation C7+ critical volumes.

#### References:

Lee, A., Gonzalez, M., Eakin, B. (1966), “The Viscosity of Natural Gases”, SPE Paper 1340, Journal of Petroleum Technology, vol. 18, p. 997-1000.

Lohrenz, J., Bray, B.G., Clark, C.R. (1964), “Calculating Viscosities of Reservoir Fluids from their compositions”, SPE Paper 915, Journal of Petroleum Technology, p. 1171-1176.

# Action Item

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

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

#### Problem Set

1. Gas metering is a most important activity in the natural gas business. Among the properties we have studied, which one would you emphasize in terms of accuracy for most common gas meters?

# Engineering Applications (I)

### Module Goals

Module Goal: To highlight some of the important applications of phase behavior in production operations.

Module Objective: To describe the use of flash calculations in separator optimization and gas-condensate reservoir description.

# Introduction: Phase Behavior as the Quintessential Tool

Optimal design and the safe and efficient operation of hydrocarbon production handling and processing systems strongly depends on accurate knowledge of fluid phase behavior. In fact, in contrast with other disciplines, the practice of petroleum and natural gas engineering centers around understanding the interaction between fluids and various environments, including the reservoir, pipeline, separator, pumps and compressors, etc. Another distinguishing characteristic is the complexity of the fluid involved here — petroleum. Whereas, for instance, mechanical engineers deal mainly with water (single-component system) and air (considered ideal for most applications), here we are dealing with complex hydrocarbon mixtures where the composition dependence of thermophysical properties is very strong.

Another complication is the wide range of pressures and temperatures associated with the processes of interest, from ultralow temperatures (LNG) to as much as 210 °F and pressures ranging from atmospheric to several thousand psia. Within these ranges, the fluid can transcend the three principal phases, namely gas, liquid, and solid, and worse yet, any combination of these.

The combination of the complex mixtures involved, the wide compositional variability from reservoir to reservoir, from system to system, and the wide range of pressures and temperatures to which systems are often subjected (e.g. a pipeline) make the phase behavior of these systems a very challenging undertaking. Unless one has a good descriptive and predictive understanding of the fluid’s phase behavior, their interactions and responses cannot be successfully described.

In these last two modules of the course, we will examine some of the applications of our current knowledge of phase behavior and thermodynamics in Petroleum and Natural Gas Engineering. The message we would like to provide is very simple: the phase behavior of the hydrocarbon system must be fully grasped in order to fully understand the responses of condensate and natural gas systems and optimize their performance. For example, maximization of condensate yield is virtually impossible without the tools for accurate prediction of just how much liquid will exist under given conditions of pressure, temperature and composition. Therefore, having advanced predictive tools for the characterization of hydrocarbon phase behavior with the highest accuracy possible is the key to mastering the economics of hydrocarbon systems. In the next sections, we will explore some specific areas where the mastering of phase behavior concepts is a must.

# Design and Optimization of Separators

Once oil and gas are brought to the surface, our main goal becomes that of transportation of the oil and gas from the wellhead to the refinery (for final processing). All equipment and processes required to accomplish this are found at the surface production facility. Hence, all surface production starts right at the wellhead. Starting at the wellhead, the complex mixture of produced fluids makes its way from the production tubing into the flow line. Normally, many wells are drilled to effectively produce the hydrocarbons contained in the field. From each of these wells emerge one or more flow lines depending on how many layers are being produced simultaneously. Depending on the physical terrain of the area and several other environmental factors, each of the flow lines may be allowed to continue from the wellhead to a central processing facility commonly referred as a production platform or a flow station, which then carries the fluids to the production platform. The combination of the wellhead, the flow lines, bulk headers, valves and fittings needed to collect and transport the produced fluid to the production platform is referred to as the gathering system.

The gathered fluids must be processed to enhance their value. First of all, fluids must be separated into their main phasial components - namely, oil, water, and natural gas. The separation system performs this function. For this, the system is usually made up of a free water knock-out (FWKO), flow line heater, and oil-gas (two-phase) separators. We will be looking at the design of this last component.

The physical separation of these three phases is carried out in several steps. Water is separated first from the hydrocarbon mixture (by means of the FWKO), and then the hydrocarbon mixture is separated into two hydrocarbon phases (gas and oil/condensate). A successful hydrocarbon separation maximizes production of condensate or oil, and enhances its properties. In field applications, this is accomplished by means of stage separation. Stage separation of oil and gas is carried out with a series of separators operating at consecutively reduced pressures. Liquid is discharged from a higher-pressure separator into the next-lower-pressure separator. The purpose of stage separation is to obtain maximum recovery of liquid hydrocarbons from the fluids coming from the wellheads and to provide maximum stabilization of both the liquid and gas effluents.

Surface Production Facility: The physical installation where fluids coming from the wellhead are separated into three main constituents: water, oil, and natural gas.

Figure 20.1: Purpose of a Surface Production Facility

Usually it is most economical to use three to four stages of separation for the hydrocarbon mixture. Five or six may payout under favorable conditions, when — for example — the incoming wellhead fluid is found at very high pressure. However, the increase in liquid yield with the addition of new stages is not linear. For instance, the increase in liquids gained by adding one stage to a single-stage system is likely to be substantial. However, adding one stage to a three or four stage system is not as likely to produce any major significant gain. In general, it has been found that a three stage separating system is the most cost effective. Figure 20.2 shows this typical configuration.

Figure 20.2: Three-stage Surface Separation Facility

Under the assumption of equilibrium conditions, and knowing the composition of the fluid stream coming into the separator and the working pressure and temperature conditions, we could apply our current knowledge of VLE equilibrium (flash calculations) and calculate the vapor and liquid fractions at each stage. However, if we are looking at designing and optimizing the separation facility, we would like to know the optimal conditions of pressure and temperature under which we would get the most economical profit from the operation. In this context, we have to keep in mind that stage separation aims at reducing the pressure of the produced fluid in sequential steps so that better and more stock-tank oil/condensate recovery will result.

Separator calculations are basically performed to determine:

• Optimum separation conditions: separator pressure and temperature
• Compositions of the separated gas and oil phases
• Oil formation volume factor
• Producing Gas-Oil ratio
• API gravity of the stock tank oil

Let us look at the case of three-stage separation. In general, temperature conditions in the surface separation facility are very much determined by the atmospheric condition and incoming stream temperatures. As for pressures, the very first separator pressure is controlled by the gathering lines coming from well heads, thus there is not much room for playing with pressure in the first separator. The same arguments are valid for the last stage of separation (stock tank), which usually operates at atmospheric conditions. Therefore, we are only left with the middle separator for optimization.

As it turns out, the key to designing a three stage separation system is finding the optimum pressure at which to operate the second separator. The question that we would answer is “what is the pressure that will result in the best quality liquid going out of the stock tank for sales?” We do not want to do this empirically. This is, we do not want to play with the second stage separator pressure in the field, until we ultimately find the optimum condition. What we can do, using our phase behavior knowledge, is to find this optimum middle stage pressure applying our understanding of VLE equilibrium.

Figure 20.3 shows the typical effect of playing with the middle separator pressure on the quality and quantity of produced oil/condensate at the stock tank. Quality and quantity are measured in terms of properties, such as API and Bo, and the overall GOR at the separation facility.

Figure 20.3: Selection of Optimum Middle Separator Pressure

The optimum value of pressure for the middle stage is the one that produces the maximum liquid yield (by minimizing GOR and Bo) of a maximum quality (by maximizing stock-tank API gravity). The smaller the value of GOR and Bo, the larger the liquid yield. The higher the API gravity of the stock-tank fluid, the more profitable its commercialization. From Figure 20.2, we see that this condition is found at neither extreme (low/high) values of middle stage pressure. There is, in fact, an optimal value for middle stage pressure. This is the value we are looking for.

The Phase Behavior model that we have described throughout these series of lectures provides the basic framework for the type of calculations required here. Additionally, we discussed how API and Bo can be calculated using the output of the phase behavior model. While doing the calculations for a 3-stage separating system, keep in mind that we have minimal control over feed pressure, as we do not want to inhibit the well (high-pressure separator). We do not control the sales line pressure (stock-tank pressure) either. The control that we do have is the operating pressure of the middle separator.

Recall that finding the optimum pressure calls for, in part, finding the minimum gas to oil ratio (GOR, SCF/STB). We are dealing, in this case, with total GOR. The total GOR is the cumulative amount of gas from all three separators divided by the amount of liquid/condensate leaving the stock tank. During our discussion on Bo-calculations, we called “nst” the moles of liquid leaving the stock tank per mole of feed entering the separation facility. This number can be obtained by sequentially flashing 1 lbmol of feed through each of the separation stages. Recalling the definition of GOR,

(20.1)

where:

${\left({V}_{o}\right)}_{sr}=\frac{{n}_{st}{\left({Z}_{o}\right)}_{st}R{T}_{st}}{{P}_{st}}$
${\left({V}_{g}\right)}_{sc}=379.4\frac{SCF}{lbmol}\cdot \left({n}_{g}\right)TOTAL=379.4\cdot \left(1-{n}_{st}\right)$ [basis: 1 lbmol of feed]

Therefore,

$GOR=\frac{5.615}{379.4}R\left(\frac{{n}_{st}}{1-{n}_{st}}\right)\left(\frac{{T}_{st}}{{P}_{st}}\right){\left({Z}_{o}\right)}_{sc}\left[SCF/STB\right]$
(20.2)

Usually, the stock tank is considered to operate at standard conditions (psc, Tsc). Now you are ready to make your own surface separation design!

# Compositional Modeling of Gas-Condensate Reservoirs:The Zero-dimensional Approach

Phase behavior is relevant to every aspect of petroleum and natural gas engineering. According to the complexity of the reservoir fluid phase behavior, reservoir modeling is classified under two distinct groups: black-oil simulation and compositional simulation. ; Often times, it is safe to assume that reservoir fluid behavior is only a function of pressure and independent of composition. This simplified behavior is typical of “black-oil” systems. In this case, reservoir hydrocarbon fluids are assumed to be comprised of two components, namely oil and gas. The model allows for certain amount of gas to be in solution with the oil at reservoir conditions. The amount of dissolved gas increases with, and is a sole function of, pressure for conditions below the bubble point. Above the bubble point pressure, the oil-component carries all the available gas in the reservoir, and a “variable bubble point” algorithm is usually implemented to predict conditions for the release of the dissolved gas. For this “black-oil” simplified model to be valid, actual oil and gas phases should maintain a fixed-composition throughout the process simulated in the reservoir. In certain cases, the assumption of fixed oil and gas composition is no longer valid: for instance, depletion of gas-condensate reservoirs and volatile-oil reservoirs, and processes that aim at the vaporization or miscible flooding of the in-situ fluids by fluid injection. More complex fluid behavior requires treating all hydrocarbon phases as nc-component mixtures and thus performing a “compositional simulation.”

Although the need for taking compositional dependence of thermodynamic and hydrodynamic parameters into account in reservoir description has been recognized for a long time, the actual implementation has not been realized until relatively recently. One of the main reasons is due to the lack of simple and reliable methods of predicting phase behavior under the conditions of interest. Two things have happened within the last three decades that have changed the situation: They are (1) the availability of fast and relatively inexpensive computational power to carry out the great number of calculations involved and (2) the development of samples and fairly good equations of state.

Simply put, a compositional reservoir simulator is a dynamic integration of the fluid dynamic porous media model and the phase behavior model, neither of which is subordinated to the other. In fact, in the early days of compositional reservoir simulation, a non-dynamic integration of these two was the norm. The more recent models actually attempt a full dynamic integration. The need for integration had been recognized earlier particularly for gas condensate reservoir, gas cycling processes, and volatile oil systems.

Early developments in reservoir engineering analysis relied on zero-dimensional or tank material balances for the evaluation and forecasting of reservoir performance. In 1936, Schilthuis devised what we now regard as the classical material balance equation. Schilthuis-type material balances are only valid for black-oil systems and are not applicable for reservoir fluids with complex behavior such as gas condensates and volatile oils. Compositional considerations were incorporated into zero-dimensional modeling in the decade of 1950’s with the works of Allen and Row (1950), Brinkman and Weinaug (1957), Reudelhuber and Hinds (1957), Jacoby and Berry (1957), and Jacoby et al. (1959). These can be regarded as the first generation of compositional simulators. Even though zero-dimensional simulations have been largely overthrown by more sophisticated numerical simulation techniques, they are still considered the most simple and fundamental tool available for the analysis of reservoir performance.

The zero-dimensional models make two principal assumptions. The first is to basically ignore the two-way coupling between fluid’s thermophysical properties and the hydrodynamic characteristics. The reservoir is treated as a perfectly mixed tank reactor with uniform properties. We assume no dimensions, and that one single value of pressure and temperature describes the average behavior of the entire reservoir. The second assumption is to neglect the hydrodynamic interactions between the flowing gas and liquid phases. In other words, the zero-dimensional compositional model relies on phase behavior as the most crucial effect controlling the description of recovery performance.

There is no doubt that significant insights, albeit qualitative, are provided by these studies. Zero-dimensional modeling provides a less expensive tool for the engineer to gain some insight into the expected performance of a gas condensate reservoir under depletion. Sometimes we spend a great deal of time dealing with the numerics of the dimensional compositional simulators, and we may forget that phase behavior is the single most important constituent of the depletion performance of, for instance, gas-condensate systems. Nevertheless, the reservoir engineer must keep in mind that this type of analysis does not provide the most accurate reservoir description. For instance, the effect of heterogeneities on reservoir performance cannot be studied with a zero-dimensional model. Again, the goal is to take a look at the qualitative insights that fluid PVT behavior can provide us.

The typical depletion sequence that a zero-dimensional compositional simulator follows is described below. This classical analysis treats gas condensate performance as constant volume depletion (CVD) in a PVT cell. The ultimate output of the model is basically comprised of GOR prediction and the compositions of the gas and condensate surface effluents. We start with a single-phase gas reservoir fluid, of a known composition, at an initial reservoir pressure and temperature. We flash this fluid several times through a series of pressure depletion stages until abandonment conditions are found. Reservoir volume is kept constant throughout the depletion process.

The typical depletion sequence that a zero-dimensional compositional simulator follows are listed below:

1. Calculate the density and molecular weight of the initial reservoir fluid. With this information, and knowing the initial reservoir volume (Vti), calculate the initial amount of gas in place (lbmols). An alternative approach is to start with a fixed amount of lbmols of reservoir fluid, and calculate the corresponding initial reservoir volume (Vti). We assume a volumetric reservoir and thus the initial reservoir volume (Vti) is kept constant throughout the calculations.
2. Using the reservoir fluid composition, distribute the initial moles of reservoir fluids into components, and store such quantities for material balance accounting.
3. Depletion step: Lower the reservoir pressure by a given amount (typically, 200 psi). Flash the reservoir fluid at the new pressure, calculate amount of moles in the gas and liquid phases (“GasT” and “LiqT”) and the densities and molecular weights of each of the phases at the new condition.
4. Expansion: With molecular weight, density, and total molar amount of each of the phases, calculate the new total volume that the fluids occupy at the new condition (Vexp).
5. Calculate the excess volume of fluids, taking the difference between the new volume upon expansion (Vexp) and the reservoir volume (Vti). This represents the volume of fluid that must have been withdrawn by the well in order to reach the newly imposed pressure condition, i.e.

Vws = Vexp – Vti(20.3)
6. Calculate the percentage of liquid in the well stream using mobility ratio considerations. Trial and error procedure is necessary for the liquid accounting. The total amount of liquid available upon depletion (LiqT) must be distributed between the wellstream (LiqWS) and the liquid remaining in the reservoir (LiqR). Additionally, the moles of reservoir liquid “LiqR” can be expressed in terms of oil/condensate saturation. Oil saturations define the mobility of the gas and liquid phases. These interrelations are shown below:
$LiqT=LiqWS+LiqR$
(20.4a)

${S}_{o}=\frac{M{W}_{o}}{{P}_{o}V{t}_{i}}LiqR$
(20.4b)

${\lambda }_{o}=\frac{{k}_{ro}\left({S}_{o}\right)}{{\mu }_{o}};{\lambda }_{g}=\frac{{k}_{rg}\left({S}_{g}\right)}{{\mu }_{g}};{\lambda }_{t}{\lambda }_{o}+{\lambda }_{g};{S}_{g}=1-{S}_{o}$
(20.4c)

${\left({V}_{ws}\right)}_{liq}=\frac{{\lambda }_{o}}{{\lambda }_{t}}{V}_{ws}=\frac{M{W}_{o}}{{P}_{o}}LiqWS$
(20.4d)

LiqWS is a function of mobility ratio, which is a function of LiqR (through So). As a first guess, take LiqR = LiqT and calculate the corresponding LiqWS. Make new guesses (decreasing the value in every new trial) until the liquid balance in Equation (20.4a) is satisfied.

7. Once LiqWS and (Vws)liq are known, obtain the total volume of gas withdrawn from the reservoir by subtracting the liquid volume (Vws)liq from the total wellstream volume (Vws). Express this gas volume in moles using reservoir gas density and molecular weight. Calculate the total number of moles of the wellstream.
8. Material Balance Accounting: Calculate the number of moles of each component remaining in the reservoir. To do this: subtract the number of moles leaving the well stream from the number in the reservoir before flashing for each of the components. Calculate the new overall composition of the reservoir fluid.
9. Calculate the overall composition of the produced well stream, by mixing the composition of gas and liquid coming along.
10. Surface Production Facility: Flash the incoming wellstream composition through the train of separations. Calculate the total amount of gas and liquid leaving the separation facility and GOR. Calculate the percentage of recovery from the reservoir.
11. A depletion loop has been completed. Go back to step 3 until abandonment pressure is reached (typically, 600 psia).
12. Plot liquid production, gas production, GOR, and recovery from the reservoir as a function of pressure depletion, from initial reservoir conditions to abandonment conditions.

The basic VLE calculations required here can be performed using the Peng-Robinson EOS and equilibrium considerations. Liquid viscosities can be calculated through the Lohrenz, Bray, and Clark correlation, and gas viscosities can be calculated using the Lee-Gonzalez correlation. Fluid densities are obtained directly through the Peng-Robinson EOS.

# Action Item

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

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

#### Problem Set

1. API, GOR and Bo are the parameters that we use to optimize a train of separators. Speculate on why all three parameters seem to predict the same optimal separator pressure.
2. Although zero-dimensional modeling can be used to model behavior of gas-condensate systems, there are some pitfalls in the description. What are these pitfalls? In what cases can we overlook them?

# Engineering Applications (II)

### Module Overview

Module Goal: To highlight some of the important applications of phase behavior in production operations.

Module Objective: To highlight the use of phase behavior for the description of gas pipelines, gas metering and hydrate systems.

# Natural Gas Pipeline Modeling

Once natural gas is produced and processed, few to several hundred miles may lie in between it and its final consumers. A cost-effective means of transport is essential to bridge the gap between the producer and consumer. In the technological arena, one of the challenges pertains to the capacity of the industry to ensure continuous delivery of natural gas while its demand is steadily increasing. Thus, it is no wonder that pipelines have become the most popular means of transporting natural gas from the wellhead to processing — and from there to the final consumer — since it better guarantees continuous delivery and assures lower maintenance costs.

Phase Behavior (P-V-T data) is crucial for all our engineering designs. Accurate prediction of the P-V-T properties of natural gases is especially critical when dealing with pipeline design, gas storage, and gas measurement. While describing natural gas pipeline design, it is necessary to distinguish between two cases: the design of pipelines for transportation of regular dry gases (no liquid, single-phase transportation) and the design of pipelines for transportation of wetter gases — where multiphase conditions due to condensate dropout may are possible.

The major variables that affect the design of gas pipelines are: the projected volumes that will be transported, the required delivery pressure (subject to the requirements of the facilities at the consumer end), the estimated losses due to friction, and the elevation changes imposed by the terrain topography. Overcoming such losses will likely require higher pressure than the one available when the gas is being produced. Thus, forcing a given gas rate to pass through a pipeline will inevitably require the use of compressor stations.

Loss in mechanical energy results from moving fluids through pipelines. Energy losses in a pipeline can be tracked by virtue of the pressure and temperature changes experienced by the flowing stream. Design equations relate pipeline pressure drop with the gas flow rate being transported. The following is the general equation for a single-phase gas pipeline flow in steady state:

${q}_{g}=c\left(\frac{{T}_{b}}{{P}_{b}}\right)\left(Eff\right)\sqrt{\frac{1}{f}}\sqrt{\frac{{d}^{5}\left({p}_{1}^{2}-{p}_{2}^{2}\right)}{g\cdot L\cdot T\cdot Z}}$
(21.1)

Note that flow rate is proportional to the inverse square root of compressibility (Z). For near-ideal conditions, the effect of compressibility on flow rate is likely to be small. But for high-pressure flows, Z may deviate greatly from 1. Under these conditions, inaccuracy in the prediction of Z may lead to a substantial error in the calculated flow rate and thus a completely wrong pipeline sizing for design purposes.

Once a pipeline is deployed, it has a more or less a fixed operational region. An upper and lower set of operational conditions allowable within the pipeline (in terms of pressure and temperature) will exist. On the one hand, the upper allowable condition will be set by the pipe strength, pipe material, diameter, and thickness. These will dictate the maximum pressure that the pipe can endure without failure (i.e., maximum operating pressure). On the other hand, maximum pressure and temperature of the compressor station discharge (which feeds the inlet of the pipe) will also contribute to set this upper level. It is clear that the conditions at the discharge of the compressor station cannot go beyond the maximum operating pressure of the pipe — otherwise the pipe will fail. The minimum or lower pressure and temperature condition of the operational region will be assigned by contractual agreement with the end consumer. The foregoing description of the operational region is shown schematically as the shaded area in Figure 21.1.

Figure 21.1: Pipeline operational curve and transported gas phase envelope

In natural gas flow, pressure and temperature changes (P-T trace) may cause formation of a liquid phase owing to partial condensation of the gaseous medium. Retrograde phenomenon — typically found in multi-component hydrocarbon systems — takes place by allowing condensation of the gas phase and liquid appearance even under expansion of the flowing stream. The same phenomenon may also cause vaporization of the liquid phase such that it reenters the gas phase. Liquid and gas phase composition are continuously changing throughout the pipe due to the unceasing mass transfer between the phases. In general, the amount of heavies in the stream determines the extent of the retrograde behavior and liquid appearance. Figure 21.1 shows a P-T trace or operational curve for a given pipeline, which is always found within the pipeline operational region.

Figure 21.1 also shows four typical phase envelopes for natural gases, which differ in the extent of their heavy components. For a given composition, the prevailing pressure and temperature conditions will determine if the fluid state is all liquid (single-phase), all gas (single-phase) or gas-liquid (two-phase). Each envelope represents a thermodynamic boundary separating the two-phase conditions (inside the envelope) from the single-phase region (outside). Each envelope is made of two curves: the dew point curve (right arm, where the transition from two-phase to single-gas occurs) and the bubble point curve (left arm, where the transition from single-liquid to two-phase occurs). Both arms meet at the critical point, which is shown in Figure 21.1. The wetness of the gas is an important concept that helps to explain the different features presented in Figure 21.1. This concept pertains to the amount of heavy hydrocarbons (high molecular weight) that are present in the gas composition. In Figure 21.1, the driest gas — i.e., the least wet — can be recognized as that whose left and right arms are the closest to each other, having the smallest two-phase region (gas A). In this figure, it can be seen that the right arm is extremely susceptible to the presence of heavies in the natural gas composition. Depending on the gas composition, the pipeline operational region can be either completely free of liquid (gas A, the driest) or partially submerged in the two-phase region (gas B, C). If the gas is wet enough, the pipeline will be entirely subjected to two-phase conditions (gas D, the wettest). One may describe the sensitivity of the right arm to heavies as having a hook-seizing effect: the larger the extent of heavies in the natural gas, the more the ‘hook’ is able to seize part of the pipeline operational region. In conclusion, since the operational region is more or less given by contractual and design considerations, the liquid presence in a pipeline is ultimately dictated by the properties of the gas that is being transported.

In the preceding figure, a pipeline handling a dry gas (gas A) will be operating a single-phase mode from its inlet through its outlet. For this case, any of the popular single-phase gas equations (Weymouth, Panhandle type, AGA equation) can be used for design purposes and to help to predict the actual operational curve (P-T trace). If a richer gas comes into the system (gas C), it will show a single-phase condition at the inlet, but after a certain distance the pressure and temperature conditions will be within the two-phase region. The case might also be that the system is transporting a wetter gas (gas D), in which case it would encounter two-phase conditions both at the inlet and at the outlet of the pipe.

Penn State has devoted a great deal of effort in the development of two-fluid models for the description of multi-phase flow condition in natural gas pipelines. In this approach, mass, momentum, and energy equations are solved simultaneously. Some simplifying assumptions are made based on engineering judgment. For instance, the knowledge of averaged flow field characteristics and fluid properties at every point of the pipeline is usually more meaningful than a detailed profile of the said properties within the cross section. Hence, generally speaking, the two-fluid model always deals with conservation equations written only in one dimension for pipeline flow (the direction of the flow along the pipe), employing cross-sectional-averaged values for each term. The use of averaged quantities absorbs the variations across the pipe section. Pressures and temperatures are assumed to be the same in both phases at any given point of the pipe. Additionally, since the main interest is to focus on normal operation conditions, the further simplification of steady state conditions is invoked.

As we have discussed, phase behavior is a crucial component in pipeline design. Not only because we need to account for gas volumetric behavior in the design equations (through, for instance, Z-factor calculations), but also because it provides a means for predicting whether multi-phase flow conditions are to be found. Liquid appearance in natural gas pipelines is as undesirable as it is inevitable. On one side, the fluid phase behavior and prevailing conditions make it inevitable. On the other, the condensate subjects the gas pipe to an increasing and undesirable energy loss. Thus, a proper pipeline design must account for the effect of condensate formation on the performance of the gas line.

# The Hydrate Problem

Natural gas hydrates are solid crystalline compounds of snow appearance with densities smaller than that of ice. Natural gas hydrates are formed when natural gas components, for instance methane, ethane, propane, isobutene, hydrogen sulfide, carbon dioxide, and nitrogen, occupy empty lattice positions in the water structure. In this case, it seems like water solidifying at temperatures considerably higher than the freezing point of water.

Gas hydrates constitute a solid solution—gas being the solute and water the solvent—where the two main constituents are not chemically bounded. Figure 21.2 presents a typical phase diagram for a mixture of water with a light, pure hydrocarbon (HC), similar to that presented by McCain (1990).

Figure 21.2: Phase Diagram for a Water/Hydrocarbon (HC) System

There are a number of points on the diagram in Figure 21.2 that are noteworthy. First of all, hydrate formation is clearly favored by low temperature and high pressure. The three-phase critical point is point C on the diagram that represents the condition where the liquid and gas hydrocarbon merge into a single hydrocarbon phase in equilibrium with liquid water. Point Q2 is the upper quadruple point, where four phases (liquid water, liquid hydrocarbon, gaseous hydrocarbon, and solid hydrate) are found in equilibrium. Point Q1, the lower quadruple point, typically occurs at 32 °F (ice freezing point) where four phases (ice, hydrate, liquid water, and hydrocarbon gas) are found in equilibrium. In this context, phases are not pure as they contain some amount of the other substances according to their mutual solubility.

For practical applications, the most important equilibrium line is the Q1Q2 segment. It represents the conditions for hydrate formation or dissociation, a critical piece of information for most industrial applications where hydrates are involved. When we focus on this zone, the phase behavior of water/hydrocarbon system is simplified to the schematics shown in Figure 21.3.

Figure 21.3: Phase Behavior of Water/Hydrocarbon System (Q1Q2 segment)

Phase Behavior thermodynamics is usually invoked for the prediction of the Q1Q2 hydrate formation/dissociation line. The first two methods of prediction were proposed by Katz and coworkers, and are known as the Gas Gravity Method (Katz, 1945) and the Ki-value Method (Carson and Katz, 1942). Both methods allow calculating the P-T equilibrium curves for three phases: liquid water, hydrate and natural gas. These methods yield initial estimates for the calculation and provide qualitative understanding of the equilibrium; the latter method being the more accurate of the two. The third method relies on Statistical Mechanics for the prediction of equilibrium. It is recognized as the most accurate of all three-phase calculations as it is more comprehensive and detailed.

The key circumstances that are essential for hydrate formation can be summarized as:

1. Presence of “free” water. No hydrate formation is possible if “free” water is not present. Here, we understand the importance of removal of water vapor from natural gas, so that in case of free water occurrence there is likelihood of hydrate formation.
2. Low temperatures, at or below the hydrate formation temperature for a given pressure and gas composition.
3. High operating pressures.
4. High velocities, or agitation, or pressure pulsations, in other words turbulence can serve as catalyst.
5. Presence of H2S and CO2 promotes hydrate formation because both these acid gases are more soluble in water than the hydrocarbons.

The best and permanent remedy for the hydrate formation problems is the dehydration of the gas. Sometimes, it is quite possible that hydrates will form at the well site or in the pipeline carrying natural gas to the dehydration unit, so that the need for well head techniques arises. At well site, two techniques are appropriate:

1. Heating the gas stream and maintaining flow lines and equipment at temperature above the hydrate point,
2. In cases where liquid water is present and the flowlines and equipment cannot be maintained above hydrate temperature, inhibiting hydrate formation by injecting additives that depress both hydrate and freezing temperatures.

The most common additives are methanol, ethylene glycol, and diethylene glycol. Methanol injection is very beneficial in cases where a low gas volume does not permit the dehydration processing. It is also extremely useful in cases where hydrate problems are relatively mild, infrequent, or periodic, in cases where inhibitor injection is only a temporary phase in the field development program, or where inhibition is done in conjunction with a primary dehydration system.

# Gas Metering

Gas measurement is another area of hydrocarbon engineering where accurate prediction of the P-V-T properties of the working fluid is especially critical. One of the most widely used meters used in the measurement of gas flow is the orifice meter. Orifice meters are classified as inferential meters because the gas volume is calculated from readings of pressure variation as the gas passes through an orifice, and it is not obtained by a direct reading.

The orifice meter is arranged so that the flowing gas is constricted at a particular location by a thin orifice plate very accurately gauged and calibrated so as to be in a concentric position in the pipe. The reduction of the cross section of the flowing gas stream in passing through the orifice increases the velocity head at the expense of the pressure head, and the reduction in pressure between the taps is measured by manometers (or a recording meter). A typical orifice meter is shown in Figure 21.4.

Figure 21.4. Orifice Meter

Among the advantages of using oriface meters for gas measurement purposes are the following facts:

• They are simple in design and have no moving parts.
• They are relatively accurate.
• They are easy to install and maintain.
• They cover a wide range of capacity.
• They represent a low cost.
• There is a great deal of experience in their use.

Among the disadvantages of oriface meters are the following facts:

• They represent an intrusive measurement technique and a flow restriction that translates into a large energy loss.
• The orifice hole can be eroded by sand or corrosive fluids.
• The hole may be obstructed by wax or hydrate.

Among the advantages of using orifice meters for gas measurement purposes are the fact that they are simple in design and have no moving parts, relatively accurate, easy to install and maintain, cover a wide range of capacity, represent a low cost, and there is a great deal of experience in their use. Among the disadvantages of orifice meters, are the fact that they represent an intrusive measurement technique and a flow restriction that translates into a large energy loss, the orifice hole can be eroded by sand or corrosive fluids, and the hole may be obstructed by wax or hydrate.

Bernoulli’s equation is then used as the basis for correlating the increase in velocity head with the decrease in pressure head. In the calculation of gas flow rate using an orifice meter, two quantities must be measured: the static pressure (i.e. the line pressure) and the differential pressure (i.e. the pressure drop across the orifice plate). The following is the basic equation for gas flow through an orifice meter:

${q}_{g}=\frac{\pi }{4}C{Y}_{1}{d}^{2}\sqrt{\frac{2\Delta p}{\left(1-{\beta }^{4}\right)\rho }}=\frac{\pi }{4}C{Y}_{1}{d}^{2}\sqrt{\frac{2\Delta pZRT}{\left(1-{\beta }^{4}\right)p\left(MW\right)}}$
(21.2)

In Equation (21.2), flow rate is a function of gas compressibility factor (Z). Again, for high-pressure flows, an error in the compressibility factor could result in an erroneously calculated flow rate. If you have some error on Z-factor, this automatically translates into error in the gas meter. Accurate phase behavior prediction techniques are a must in gas metering.

In the Natural Gas Industry, the point of gas exchange between the buyer and the seller is called custody transfer. During custody transfer operations, accurate measurements of the quantity and quality of the exchanged gas are of crucial importance because of its economical implications. Economic transactions are based on volumetric rate measurements, which are regulated to be made at the same base conditions. Industry base conditions or standard conditions (SC) are usually taken as P = 14.7 psia and T = 60.0 °F. A low percent inaccuracy in the Z-factor calculation of a gas in transfer can easily translate into thousands of dollars of losses on a daily basis! In fact, flow rate estimations can prove extremely sensitive to values of compressibility factor. This is why the gas industry does not accept Z-factor predictions with a range of uncertainty larger than + 0.01 % for custody transfer operations.

# Action Item

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