You must log in to edit PetroWiki. Help with editing

Content of PetroWiki is intended for personal use only and to supplement, not replace, engineering judgment. SPE disclaims any and all liability for your use of such content. More information

PEH:Thermodynamics and Phase Behavior

Jump to navigation Jump to search

Publication Information


Petroleum Engineering Handbook

Larry W. Lake, Editor-in-Chief

Volume I – General Engineering

John R. Fanchi, Editor

Chapter 7 – Thermodynamics and Phase Behavior

R.T. Johns, U. of Texas at Austin

Pgs. 333-369

ISBN 978-1-55563-108-6
Get permission for reuse

Phase behavior describes the complex interaction between physically distinct, separable portions of matter called phases that are in contact with each other. Typical phases are solids, liquids, and vapors. Phase behavior plays a vital role in many petroleum applications, such as enhanced oil recovery, compositional simulation, geochemical behavior, wellbore stability, geothermal energy, environmental cleanup, multiphase flow in wellbores and pipes, and surface facilities.

Thermodynamics, which is central to understanding phase behavior, is the study of energy and its transformations. Using thermodynamics, we can follow the energy changes that occur during phase changes and predict the outcome of a process. Thermodynamics began as the study of heat applied to steam power but was substantially broadened by Gibbs in the middle to late 1800s. Gibbs’ most significant contribution was the development of phase-equilibrium thermodynamics applied to multicomponent mixtures, particularly the concept of chemical potential.[1] The concept of chemical potential leads to the result that, at equilibrium, the chemical potential of each component must be the same in all phases (μiL = μiV).

Phase-equilibrium thermodynamics seeks to determine properties such as temperature, pressure, and phase compositions that establish themselves once all tendencies for further change have disappeared. This chapter reviews the fundamentals of phase-equilibrium thermodynamics used in petroleum applications, especially those that require liquid-vapor phase behavior.[2][3][4][5][6][7] The next chapter in this section of the Handbook illustrates phase behavior through diagrams.

Fundamental Ideas and Problem Statement

Fig. 7.1 is a schematic showing a closed container of liquid and vapor. Given constant and known temperature, pressure, and overall compositions (zi where i = 1, … , nc) at equilibrium, the fundamental task is to quantify the molar fractions of the phases (L, V) and compositions of the vapor (yi where i = 1, … , nc) and liquid phases (xi where i = 1, … , nc) that form at equilibrium. The phases are assumed to be homogeneous, in which intensive parameters such as pressure, temperature, density, viscosity, and phase compositions are uniform throughout the phase. (Thus, gravity effects are not typically considered.) Intensive properties are those that are independent of the amount of the phases (e.g., phase density, pressure, and temperature). Alternatively, extensive properties depend on the amount of the phases (e.g., total volume and moles of liquid). Intensive properties can be determined as the ratio of two extensive properties; for example, molar density is the number of moles divided by the total volume.

The overall compositions and phase compositions in Fig. 7.1 are written as mole fractions, which are defined by




where RTENOTITLE for the container, RTENOTITLE for liquid, and RTENOTITLE for vapor. The relative amounts of the phases are defined by the phase mole fractions,




where L + V = 1. The phase molar fractions are not saturations, although they could be converted to saturations from the phase densities. The molar fractions of the phases are related to the overall and phase compositions by


Thus, once the overall compositions and phase compositions are known, the phase molar fractions, L and V, are also known.

Gibbs Phase Rule and Duhem’s Theorem

The Gibbs phase rule and Duhem’s theorem assure us that the problem illustrated in Fig. 7.1 can be solved. Ideas and theories from thermodynamics are based on observations. Gibbs, for example, observed that the equilibrium intensive state of the system is fully known once the pressure, temperature, and phase compositions are specified. The number of intensive properties that we would like to know is, therefore, 2 + ncnp, where np is the number of phases (for vapor/liquid equilibrium, np is two). These intensive properties can only be determined if a sufficient number of equations are available or if some of them are explicitly specified. An inventory of equations shows that there are np summation equations (i.e., the phase mole fractions for each phase sum to 1.0) and nc(np - 1) equilibrium relations, for a total of np + nc(np - 1) equations. The equilibrium relations could be given as K-values RTENOTITLE, which relate the component liquid and vapor mole fractions or, as described later, chemical potential criteria for equilibrium (i.e., μiL = μiV).

The Gibbs phase rule says that the degrees of freedom are 2 + nc - np, which is the difference between the number of required intensive properties (unknowns) and the number of relations (equations). The Gibbs phase rule is only practically useful for a small number of components but does offer significant insight into the maximum number of phases that can form as well as how many intensive properties can be independently specified.

For example, suppose that only one phase (np = 1) is present at equilibrium in a system containing a pure fluid (nc = 1). The Gibbs phase rule says that only two intensive properties can be specified (degrees of freedom are two). We cannot specify three or more intensive properties for this case, but we are free to choose which intensive properties are set. Typically, we would choose temperature and pressure. The choice of intensive properties is not completely arbitrary, for only properties related to an individual phase can be selected. Thus, properties such as the overall density of the two-phase system or the phase molar fractions, L and V, cannot be used to reduce the degrees of freedom.

Suppose next that three equilibrium phases exist in the pure fluid (i.e., the triple point). For this case, the degrees of freedom are zero, and no intensive properties can be specified. That is, the intensive properties, such as temperature and pressure, are determined and are not arbitrary at the triple point. Four phases in equilibrium with each other are not allowed by the Gibbs phase rule (neither are they observed experimentally).

Duhem’s theorem is another rule, similar to the phase rule, but it specifies when both the extensive and intensive states of the system are determined. The theorem states that for any closed system containing specified moles of nc components (from which the overall compositions can be calculated), the equilibrium state is completely determined when any two independent properties are fixed. The two independent properties may be either intensive or extensive; however, the maximum number of independent intensive properties that can be specified is given by the Gibbs phase rule. For example, when the degrees of freedom are one, at least one of the two variables must be extensive. When the degrees of freedom are zero, both must be extensive. Thus, the combination of the Gibbs phase rule and Duhem’s theorem shows that the extensive and intensive state of the two-phase problem in Fig. 7.1 can be determined when the temperature, pressure, and moles of all the components are specified. For a pure component system, the intensive and extensive state of the system is determined when the temperature, pressure, and the total number of moles are given.

Equilibrium, Stability, and Reversible Thermodynamic Systems

Thermodynamics is a macroscopic viewpoint in that it concerns itself with the properties of a system, such as temperature and density. Thermodynamics predicts the nature of a new equilibrium state—not the rate at which that state is reached. One of the characteristics of equilibrium is that the thermodynamic properties are time invariant. Furthermore, once equilibrium is reached, the process or pathway that led to equilibrium cannot be determined.

The equilibrium state is always time invariant, whether it is dynamic or static. A dynamic equilibrium process is a steady-state process, in which the properties change spatially but not temporally. A static process, while having the appearance of reaching a static state on a macroscopic scale (such as that shown in Fig. 7.1), is anything but static on a microscopic scale. Molecules from the liquid phase continue to move into the vapor phase and vice versa, but the rates of energy and mass transfer are equal, giving the appearance of equilibrium on a macroscopic scale. Indeed, this is exactly the definition for equilibrium embodied by the chemical potential criterion μiL = μiV.

The criterion of time invariance is a necessary, but not sufficient, condition for equilibrium. Some systems can exist in metastable states that are time invariant. For example, at the Earth’s surface, diamonds are in a metastable state of pure carbon, whereas graphite is the equilibrium state. Fig. 7.2 illustrates the concept of equilibrium vs. metastable or unstable (nonequilibrium) states by considering a ball rolling down a hill into a valley. When the ball is on the side of the hill, it is unstable and will roll down the slope because of gravitational forces; this is an unstable process. The ball, however, if initially trapped in a small depression on the side of the hill, will not roll down the hill; this is a metastable state. Lacking any additional energy, the ball will stay in the metastable position. If the depression is removed, or the ball is slightly moved, the ball will roll down the hill until it reaches the lowest position, which corresponds to the lowest gravitational potential energy or equilibrium position.

Later on, we will find that a definition for equilibrium is when the Gibbs free energy of the system is the lowest value possible, and this is how one can recognize unstable or metastable states from the true equilibrium. Generally, equilibrium states that arise naturally are stable to small disturbances. On the flip side, metastable equilibrium states, which are not stable to small disturbances, do not occur often in nature. Our mathematical description of equilibrium, however, will exhibit these unstable and metastable states, so we must be able to recognize them.

Processes of interest to us are often not time invariant, and it would appear that equilibrium thermodynamics is not very useful. For example, we typically run transient simulations to estimate the recovery of reservoir oil by injection of a gas. The concept of local equilibrium and reversibility are used to overcome this apparent limitation of thermodynamics. Equilibrium at a point in a reservoir, termed local equilibrium, often applies when internal relaxation processes are rapid with respect to the rate at which changes are imposed on the system. That is, equilibrium thermodynamics can be applied over small volumes of the reservoir, even though pressure and other gradients remain in the reservoir. In reservoir simulation, the small volumes are gridblocks, although the size of the gridblocks must be sufficiently small so that good accuracy is obtained.

The concept of reversibility of a process is also important. A reversible process proceeds in sufficiently small steps so that it is essentially in equilibrium at any given time (i.e., the process at a point in the reservoir proceeds in a succession of local equilibrium steps). A process is reversible when its direction can be reversed at any point by an infinitesimal change in external conditions.

The concept of reversibility is, in a sense, the temporal equivalent to the spatial concept of local equilibrium. Thus, the concepts of local equilibrium and reversibility allow the application of equilibrium thermodynamics to real systems, which are invariably nonequilibrium at large scales. For most cases, very little accuracy is lost in making such assumptions.

Fundamental Equations

Relatively few ideas and equations are used to solve the phase behavior problem illustrated in Fig. 7.1. The most fundamental idea in thermodynamics is the conservation of total energy, which is termed "the first law of thermodynamics." The first law is based on our every day observation that for any change of thermodynamic properties, total energy, which includes internal, potential, kinetic, heat, and work, is conserved. The second fundamental idea in thermodynamics is the total entropy balance or "the second law of thermodynamics." Entropy is a thermodynamic property that expresses the unidirectional nature of a process and, in some sense, is "nature’s clock." For example, a cup of hot coffee at room temperature cools down instead of heating up. The conservation of total mass is also used to constrain thermodynamic processes.

These equations are applied to a thermodynamic system. A thermodynamic system is defined as that part of the universe we are considering—for example, the inside of the container in Fig. 7.1. Everything else is called the surroundings. A system may be related to its surroundings in a variety of ways, depending on whether mass or energy (in the form of heat or work) is exchanged (see Fig. 7.3). When no heat or mass is transferred, and no work is done on or by the surroundings, the system is referred to as an "isolated" system. When only energy is exchanged between the system and surroundings, the system is "closed." Last, the system is "open" when both mass and energy are exchanged between the system and its surroundings. No work is allowed on or by an isolated system, and its boundaries are therefore rigid.

A thermodynamic state is given by its thermodynamic properties (e.g., pressure, density, enthalpy, temperature, internal energy, entropy, and other properties). All of these are state functions that depend only on the present state reached (point conditions)—not the path that the system took to reach that state. For example, if methane is heated, compressed, and then returned to its initial volume and temperature, the methane will have exactly the same pressure as before, independent of how it was heated or compressed. The usefulness of state functions is the simplest possible path can be selected for the calculation of the change in a state function; that is, we would likely choose a reversible path that consists of isothermal or isobaric steps.

In contrast to state functions such as entropy or pressure, heat and work are not thermodynamic properties but depend on the nature or path of the process that the system undergoes. A different path will give a different amount of work and heat.

First Law of Thermodynamics

We begin with the first law of thermodynamics applied to an open thermodynamic system. As illustrated in Fig. 7.3, an open system allows mass and energy to flow into or out of the system. We make the following assumptions and definitions:

  1. Mass flows into or out of the system along one boundary of the system. The mass flow rate into the system is positive, whereas flow rates out of the system to the surroundings are negative.
  2. Mass can carry internal energy into or out of the system. We neglect kinetic and potential energy carried by the mass. This is often a good assumption when the fluid is not moving near the speed of sound, the change in height over the system is not large, or the system temperature variations are not large.
  3. The only types of work that are present are expansion/compression of the system and flow work. The boundaries of the system can expand or contract. Thus, work can be done by the system on the surroundings or vice versa. Work is positive when the surroundings do work on the system (i.e., the system contracts). The mass that enters or exits the system also does work—sometimes called flow work or pressure work.
  4. We neglect potential energy and kinetic energy changes within the system.
  5. Energy in the form of heat might enter or leave the system across the system boundaries. Heat transfer is positive when heat is exchanged from the surroundings to the system.

Before proceeding, we must define internal energy. Internal energy of a substance is the sum of the potential energy arising from chemical bonds of atoms and electrons and the sum of the kinetic energy of the atoms and molecules. The microscopic kinetic energy is sometimes called thermal energy, which is proportional to temperature.

With this definition, a total macroscopic energy balance in the system, at an instantaneous point in time, gives




and nU, the total internal energy, is equal to the total energy within the system by assumption four previously discussed. The property, U, is the molar internal energy (total energy/mole). Eq. 7.3 shows that when work or heat is added to the system, the molecular activity increases, causing the total internal energy to increase; that is, RTENOTITLE.

The term on the right side of Eq. 7.2 contains three terms: mass influx into the system that carries energy; heat transfer into the system; and compression work done by the surroundings on the system. Because we neglect potential and kinetic energy of the mass that flows into the system (assumption two), the energy associated with the mass influx into the system is simply Uṅ, where is the molar flow rate. Based on assumption one, there is only one molar flow rate into the system. The rate of heat flow from the surroundings across the system boundaries into the system is given by RTENOTITLE. Compression (or expansion) of the system boundaries causes work on the system denoted by Ẇ. Substitution of these terms into Eq. 7.2 gives


where the left side refers to energy within the system and the right side to energy that flows across the system boundaries into the system.

The two types of work considered are expansion/compression work and flow work (assumption three). From physics, work is performed whenever a force acts over a distance. Thus, the differential mechanical work that results from a differential displacement RTENOTITLE is given by RTENOTITLE. For expansion/compression work, the external force is equal to an external pressure supplied by the surroundings multiplied by the corresponding area along the boundary of the system.

In Cartesian coordinates, dW = - pext(Axdx + Aydy + Azdz) , where the external pressure is constant along the boundary of the system. Ax is the area normal to the x-coordinate that is being displaced, and so forth. The minus sign indicates that work is positive if the displacement is negative (i.e., an external force compresses the system). The expression for the differential work can be simplified further as RTENOTITLE, where nV is the total volume and V is the molar volume (volume/mole). The rate of work is RTENOTITLE. This equation applies to any arbitrarily shaped system.

For example, consider a rectangular box that expands differentially into the surroundings on three sides, as illustrated in Fig. 7.4. Here, Axyz, Ayxz, and Azxy, where the differential cross terms are neglected. The differential work is, therefore, dW = -pext(yzdx + xzdy + xydz) , provided the external pressure is the same on all faces of the box. The differential displacement volume is equal to d(nV) = d(xyz) = yzdx + xzdy + xydz, which gives the desired result. In this example, the differential work is negative because the system does work on the surroundings. The example also illustrates that even though the system expands into the surroundings, the work is always related to the external pressure. If the external pressure is zero, no work will be done by the system because the surroundings will offer no resistance.

Flow work is done by mass that enters or exits the system. A flowing fluid element does work on the fluid ahead of it, and the fluid behind it does work on that fluid element. Flow work, for example, turns the turbine shaft of a hydroelectric power plant. For one-dimensional (1D) inflow or outflow of fluid, the instantaneous rate of flow work is RTENOTITLE, where u is the velocity of the fluid; ρ is the molar density (i.e., inverse of V); and the molar flow rate is = ρAu. Flow work is positive when the fluid is entering the system; that is, the surroundings do work on the system.

For reversible displacements, the pressure in the system must equal the external pressure, p = pext, supplied by the surroundings, so the system and surroundings are always in equilibrium. With the assumption of reversibility, the total rate of work (expansion/compression work plus flow work) becomes RTENOTITLE. Eq. 7.4 is then written as


Eq. 7.5 can be simplified by defining the enthalpy, H = U + pV. The definition for enthalpy is defined strictly for mathematical convenience. For liquids and solids at low pressures, we often take H = U because the product pV is small compared to U (the molar volume of the condensed phases is also small). Combining the first and third terms on the right side of Eq. 7.5 gives


Eq. 7.6 can be written in thermodynamic shorthand as


For closed systems, dn = 0 (the total moles in the system is constant), and Eq. 7.7 becomes


For isolated systems, dn = 0, dQ = 0, and d(nV) = 0; therefore, Eq. 7.7 reduces to


which shows that the total internal energy of an isolated system is constant (nU=constant).

Second Law of Thermodynamics

Conservation of total mass and energy are insufficient to solve many phase-equilibrium problems. Processes that satisfy these conservation equations may not be physically possible; that is, the process of a cold cup of coffee spontaneously heating up on your dinner table would satisfy the first law of thermodynamics but has a near zero probability to occur. Processes have a natural direction to them in that spontaneous processes tend to dissipate gradients (e.g., Darcy’s law for pressure gradients and Fick’s law for concentration gradients) in the system until equilibrium is reached. A system that is not subject to forced flows of mass or energy from its surroundings will evolve to a time-invariant state that is uniform or composed of uniform subsystems—the equilibrium state. The second law of thermodynamics introduces a new thermodynamic property, entropy, and provides a mathematical statement that describes this unidirectional nature of processes.

The second law also has implications for the efficiency of processes. Heat and work are not of the same quality in that work can be efficiently converted to thermal energy (e.g., frictional heat losses), but thermal energy can be only partially converted into mechanical energy (e.g., steam power plants). Thus, work is a more valuable form of energy than heat—work has a high quality. Furthermore, energy at higher temperatures is more useful than energy at lower temperatures. For example, the ocean contains an immense amount of energy, but it is not very useful because of its low temperature. Energy is degraded when heat transfers from one system to another of lower temperature. Entropy is a measure of the energy degradation or disorder of the system.

Entropy is a thermodynamic property just like temperature and pressure. Entropy is a state function, in which changes during a reversible process in a closed system are given by the ratio Q/T. Entropy increases as T decreases or Q increases.

Entropy is related to the likelihood that equilibrium will be reached. Entropy is best understood by examining a very simple example at the microscopic scale. Fig. 7.5 shows the initial state of a hypothetical closed system that contains four molecules. The system is initially partitioned into two halves, such that the molecules from one half cannot move into the other half. One molecule is in the left subsystem and three are in another, thus, the pressures are not initially the same. Each subsystem has only one possible configuration—the initial state. Thus, the subsystems are well ordered, and entropy is initially small.

When the partition is removed, however, the molecules from each subsystem are free to move into the other half of the system, and a total of 16 different configurations are possible, as shown in Fig. 7.5. Because each of these configurations is equally likely, the probability that the system will be found in its original configuration is only 1/16 (i.e., 2–4). The system is more likely to contain two molecules in each subsystem (probability of 6/16 or 3/8), which is the equilibrium state—the most disordered state. Entropy is related to the maximum number of possible configurations of the system, and thus, the entropy after the partition is removed has increased. The configurations could also be arrangements of energy quanta, instead of molecular arrangements, as in this example.

Although the original configuration for four molecules is not improbable, real systems contain many more molecules. For example, if one mole of a gas were present in the system (6.0 × 1023 molecules), the likelihood that the system would be found in its initial state would be very unlikely (2−6 × 1023). However, the probability that the system would contain a similar number of molecules on each side would be near 1.0 (i.e., the pressure would be equal throughout the closed system).

The steps to write the entropy balance for an open system are similar to those for the first law of thermodynamics. We allow for the following:
  1. Mass flows into or out of the system at only one location on the boundary of the system. The mass flow rate into the system is positive, whereas flow rates out of the system to the surroundings are negative.
  2. Mass can carry entropy into or out of the system. The rate of entropy transfer into the system by mass flow is given by Sṅ.
  3. Energy in the form of heat may enter or leave the system across the system boundaries at a specified exterior temperature, T. Heat transfer is positive when heat is exchanged from the surroundings to the system. Entropy, Q/T, is transferred from the surroundings to the system during heat transfer. The rate of entropy transfer into the system by heat transfer is therefore RTENOTITLE.
  4. The temperature at the boundary or exterior of the system is equal to the temperature within the system (i.e., the heat transfer process is reversible).
  5. Unlike total energy or mass, entropy is generated within a system that is not in equilibrium. Entropy generation is related to irreversibilities in the system such as temperature gradients, pressure gradients, or concentration gradients. The second law of thermodynamics states that entropy generation is always positive.

With these assumptions and definitions, the entropy balance is


where the term on the left is the change in total entropy within the system (S is entropy/mole), and the first two terms on the right are the net rate of total entropy transported into or out of the system by mass and heat transfer. The last term is the internal generation rate of entropy within the system. In thermodynamic shorthand, the entropy balance can be written as


For closed systems (dn = 0), the entropy balance becomes RTENOTITLE. For isolated systems (dn = 0 and dQ= 0), the entropy balance is


In an isolated system at equilibrium, the total entropy cannot change with time. Thus, the generation of entropy must be zero at equilibrium (dSG = 0), which we stated previously. Away from equilibrium, entropy generation is positive (dSG > 0), which implies that entropy, in an isolated system, increases with time and reaches a maximum at equilibrium.

Why is entropy generation positive away from equilibrium? Consider an isolated system that is composed of two open subsystems A and B.[4] Heat is exchanged only from the high temperature subsystem A to the low temperature subsystem B. Thus, the rate of heat transfer is RTENOTITLE, where h is a constant heat-transfer coefficient. Each subsystem is well mixed so that T is always uniform [i.e., no internal gradients exist, and subsystems A and B must pass through a succession of equilibrium states (i.e., the process is reversible and GA = GB = 0)]. From Eq. 7.9 and the expressions for the rate of heat transfer, the entropy balance for subsystem A is RTENOTITLE, and for subsystem B, RTENOTITLE.

For the isolated system not in equilibrium, the entropy balance is RTENOTITLE. Furthermore, the total entropy, nS, is an extensive property so that nS = nASA + nBSB. The entropy generation term of the combined system is therefore given by RTENOTITLE. Because absolute temperatures are positive and the heat transfer coefficient is positive, this result demonstrates that entropy generation is always positive away from equilibrium. At equilibrium TA = TB, and the entropy generation term is zero, as stated before. Furthermore, the rate of entropy generation is proportional to the square of the gradients (temperature difference in this case). If the temperature gradients are kept small with time (infinitesimal equilibrium steps), entropy generation remains near zero.

A process is reversible if it occurs with small gradients (i.e., consists of a succession of infinitesimal equilibrium steps). Thus, dSG = 0 for a reversible process, and from Eq. 7.10,


For a closed system that undergoes a reversible process, Eq. 7.12 reduces to dQrev = Td(nS).

In summary, the second law of thermodynamics states that the rate of entropy generation within a system must be greater than or equal to zero. A process for which the rate of generation of entropy is always zero is a reversible process. A large rate of entropy generation corresponds to greater process irreversibilities.

Fundamental Property Relations and Equilibrium Conditions

We would like to calculate thermodynamic properties and the equilibrium state from a simplified mathematical model, called an equation-of-state (EOS). To do this, we need equations that relate thermodynamic quantities in terms of pressure, molar volume, and temperature data (PVT data), and we want to eliminate any path dependence by eliminating all properties that are not state functions. Substitution of Eq. 7.10 into Eq. 7.7 by elimination of dQ (a path dependent quantity) and selection of a reversible path (such that dSG = 0) gives


All of the properties in Eq. 7.13 are state functions; thus, Eq. 7.13 is independent of the path or process. After combining like terms, Eq. 7.13 becomes


where GH - TS is defined as the molar Gibbs energy. For a closed system (dn = 0), Eq. 7.14 becomes


Eqs. 7.14 and 7.15 are examples of fundamental property relations. Other fundamental property relations are possible. For example, differentiation of the definition for total Gibbs energy gives d(nG) = d(nH) - Td(nS) - (nS)dT. Similarly, differentiation of the total enthalpy gives d(nH) = d(nU) + pd(nV) + (nV)dp. Substitution of these relations into Eq. 7.14 (or Eq. 7.15) by elimination of the enthalpy term gives the fundamental property relation for the total Gibbs free energy of a closed system as


or, for an open system,


Equilibrium Criteria for Single-Component Liquid/Vapor Systems. Consider an isolated system of a pure fluid with two phases, vapor and liquid. Initially, the temperature, pressure, and other properties of the two phases are not in equilibrium. Fig. 7.6 illustrates the composite isolated system in which each phase is treated as a subsystem.

We begin by writing the differential entropy change from Eq. 7.14 for each open subsystem. The vapor phase equation is


and for the liquid phase,


The change in the total entropy of the isolated system can be written as the summation of Eqs. 7.18 and 7.19. The result is


For an isolated system, the change in the total internal energy is zero (see Eq. 7.8), as is the change in the total mass and volume. Thus, dnL = -dnV, d(nU)L = -d(nU)V, and d(nV)L = -d(nV)V. The differential entropy change for an isolated system at equilibrium must also be zero (see Eq. 7.11). Eq. 7.20 becomes


Because changes in internal energy, volume, and mass of the liquid phase can be arbitrarily set (i.e., are independent), we must have at equilibrium that RTENOTITLE. Thus, at equilibrium, TL = TV, pL = pV, and GL = GV. The first two equilibrium criteria are obvious. The equilibrium condition that the Gibbs free energy of the phases is equal is not as obvious.

Other systems lead to similar equilibrium conditions. For example, for a closed system at constant pressure and temperature (dp = 0, dT = 0) the fundamental property relation from Eq. 7.16 becomes d(nG) = 0. Thus, the equilibrium criterion here is that the Gibbs free energy must be a minimum. This criterion also leads to the equality of the Gibbs free energy of both phases at equilibrium, GL = GV.

Fugacity of a Pure Fluid. Fugacity criterion is often used as a substitute for the Gibbs free- energy criterion. The definition for fugacity comes from an analogue with ideal gases that is derived for a closed system under isothermal conditions. Eq. 7.16 for an isothermal process (dT = 0) is


For an ideal gas, RTENOTITLE, and Eq. 7.21 becomes


Fugacity is defined by analogy for a fluid that is not ideal. That is, we define the fugacity, f, based on a comparison with Eq. 7.22, which is written as


Eq. 7.23 shows that the value for fugacity is whatever is required to give the correct behavior of the real fluid. More exactly, fugacity measures how the Gibbs free energy of a real fluid deviates from that of an ideal gas. Fugacity has units of pressure, and for an ideal gas the fugacity is equal to the pressure (compare Eqs. 7.23 and 7.22).

We showed that at equilibrium for a pure fluid GL = GV. By integration of Eq. 7.23 under isothermal conditions, the Gibbs free-energy criterion implies that the fugacity of the liquid and vapor phases must also be equal at equilibrium. That is, at equilibrium for a pure fluid,


We would like an expression for fugacity in terms of our convenient quantities of pressure, molar volume, and temperature, so that an EOS can be used. Substitution of Eq. 7.21 into Eq. 7.23 gives


Subtraction of RTENOTITLE from both sides and some algebraic rearrangement gives


Finally, integration from a reference state of zero pressure (ideal gas state) to the actual pressure gives


From the definition of fugacity, RTENOTITLE (i.e., the fugacity is equal to the pressure for an ideal gas), we have


where RTENOTITLE is known as the fugacity coefficient, and RTENOTITLE is the compressibility factor. The fugacity coefficient is therefore equal to 1.0 for an ideal gas. Eq. 7.25 requires knowing the compressibility factor as a function of pressure.

Models for compressibility factor, such as a cubic EOS, however, are typically not explicit functions of pressure. A more convenient form would be to transform the integral with respect to pressure to one with respect to volume. Eq. 7.25 can be transformed to


The importance of Eq. 7.26 is that the fugacity can be calculated if the molar volume, temperature, and pressure are known over the full range of molar volumes from V to ∞ . Typically, sufficient laboratory data (p, V, T) is not available, and mathematical models, such as cubic EOS, are used. We will use Eqs. 7.24 and 7.26 in Sec. 7.4 to calculate the intensive and extensive state of a pure fluid at equilibrium using a cubic EOS.

Equilibrium Criteria for Multicomponent Liquid/Vapor Systems. The procedure to determine the equilibrium criterion for multicomponent systems is similar to that used for pure fluids. We consider a closed system with a multicomponent mixture of n moles as illustrated in Fig. 7.1. Transfer of mass from one phase to the other is allowed, but the overall system is closed, such that the overall composition of the system is constant. Given the overall compositions (zi) , pressure, and temperature, we seek to determine the amount of liquid and vapor present at equilibrium, as well as the component mole fractions for the phases (xi and yi).

As before, the closed system consists of two subsystems, the liquid and vapor phases (see Fig. 7.6). The primary difference between the derivation for pure fluids and the derivation for multiple components is that the fundamental property relations for the open system must be modified to include mass transfer of different components. That is, we must compute the change in the total Gibbs energy of the liquid phase as small amounts of each component (dni) are transferred from the vapor phase to the liquid phase (or vice versa for the vapor phase). For example, Eq. 7.17 becomes


where RTENOTITLE is the molar Gibbs free energy added to the liquid phase when dni moles are added to it. The partial molar Gibbs energy RTENOTITLE is also named the chemical potential, RTENOTITLE. The chemical potential measures how much Gibbs energy is added to a mixture when dni is added to it.[1] Thus, Eq. 7.27 is commonly written as


for the liquid phase, and


for the vapor phase. As for pure fluids, these two equations are added to obtain the differential total Gibbs energy of the entire closed system. Because the differential total Gibbs free energy of the closed system must be zero when pressure and temperature are constant, we obtain


Conservation of mass requires that any component that enters a phase must have come from the other phase so that dniL=-dniV, and upon substitution into Eq. 7.28,


Because the dni are independent and arbitrary, we must have at equilibrium


Eq. 7.29 says that at equilibrium the chemical potential of a component in the liquid phase must be equal to the chemical potential of the same component in the vapor phase. This equilibrium criterion reduces to GL = GV for the case of a pure fluid.

Fugacity of a Component in a Mixture. The equilibrium criterion expressed as component fugacities is often used instead of chemical potentials. The reason for this is primarily one of convenience because component fugacity has units of pressure. Just as for pure fluids, the fugacity of a component is defined as an analogue to an ideal gas mixture.

Consider an ideal gas mixture at a temperature T. The pressure for n moles is RTENOTITLE. In this mixture, each component has ni moles. If ni moles of each component in this mixture occupy the same total volume alone at the same temperature, the pressure would be RTENOTITLE. Division of this result by the pressure gives the partial pressure of a component in an ideal gas mixture. That is, pi = yip, where RTENOTITLE is the vapor molar fraction of each component. The sum of the partial pressures equals the pressure RTENOTITLE.

For an ideal pure gas at constant temperature, we had RTENOTITLE (see Eq. 7.22). It follows, therefore, that the partial molar Gibbs energy of a component should be evaluated at the partial pressure, or


For a real mixture (not an ideal gas or solution), the component fugacity is defined by analogue with Eq. 7.30.


where RTENOTITLE is the fugacity of component i in a mixture. The component fugacity for real fluids is sometimes referred to as a corrected partial pressure. Comparison of Eqs. 7.30 and 7.31 show that for ideal mixtures, RTENOTITLE. From the integration of Eq. 7.31 and the use of the equilibrium criteria of Eq. 7.29, we obtain the equilibrium criteria for component fugacities as


Eq. 7.32 is often used instead of the equality of chemical potentials to determine equilibrium.

To calculate component fugacities of a real mixture, we subtract the chemical potential for component i in an ideal gas (Eq. 7.30) from both sides of Eq. 7.31. The result is


where RTENOTITLE is the component fugacity coefficient. Eq. 7.33 is used to calculate the deviation of the component fugacity from ideal behavior (this is also known as the residual partial Gibbs energy of component i). Integration of Eq. 7.33 from zero pressure to the actual pressure gives RTENOTITLE, where the chemical potential is zero and the component fugacity coefficient is 1.0 at zero pressure (the mixture is ideal at zero pressure). From Eqs. 7.21 and 7.22 and the definition of fugacity, we obtain using calculus:


where RTENOTITLE is the partial molar compressibility factor.

Eq. 7.34 is similar in form to Eq. 7.25 for a pure fluid. Table 7.1 compares the fundamental equilibrium equations for pure and multicomponent fluids. Because cubic EOS represent Z as an explicit function of V and not Z as a function of p, Eq. 7.34 is often rearranged to


Eq. 7.35 shows that the fugacity of a component in a mixture can be calculated when the molar volume, temperature, pressure, and compositions are known over the full range of molar volumes from V to ∞ . Sufficient laboratory data is typically not available for the integration and an EOS must be used. We will use Eqs. 7.32 and 7.35 in Sec. 7.5 to calculate the intensive and extensive state of a multicomponent mixture at equilibrium using a cubic EOS.

Volumetric Properties of Pure Fluids

The phase behavior of a typical pure fluid can be represented on a pressure-temperature diagram, as illustrated in Fig. 7.7. From the Gibbs phase rule (Sec. 7.2.1), the number of degrees of freedom are 3 - np, which means that one, two, or three phases can be present at equilibrium. For simplicity, these phases are shown as a solid, liquid, and a vapor, although numerous additional solid and liquid phases are possible, as long as no more than three of those phases coexist at equilibrium at the same temperature and pressure. Water, for example, has nine different solid phases, each of which has a different crystalline structure.[6] At the pressure and temperature of the Earth’s surface, however, we experience only one solid, liquid, and vapor phase of water.

According to the Gibbs phase rule, there are no degrees of freedom when three phases are in equilibrium. This necessarily implies that three phases must be in equilibrium only at one temperature and pressure; this is the triple point indicated in Fig. 7.7. That is, we are not free to choose the temperature and pressure at which three phases can coexist. The temperature and pressure of the triple point are determined.

For two phases, however, the degrees of freedom are one, and we can set either the temperature or the pressure, but not both. Once the temperature is specified, the pressure is determined. Thus, two equilibrium phases are represented by curves on a pressure-temperature diagram. The sublimation curve gives the locus of points where solid and vapor are in equilibrium; the melting or fusion curve gives the locus of points where solid and liquid are in equilibrium; and the vaporization or saturation curve gives the locus of points where the liquid and vapor are in equilibrium. The pressure along the vaporization curve is called the vapor or saturation pressure. For example, the pure fluid at temperature T1 in Fig. 7.7 has a vapor pressure of RTENOTITLE.

A single solid, liquid, or vapor phase can exist over a range of temperatures and pressures (degrees of freedom are two). This is indicated in Fig. 7.7 by the three single-phase regions of solid, liquid, and vapor.

When the temperature and pressure along the vaporization curve is increased to the critical pressure, pc, and temperature, Tc, the interface between the liquid and vapor phases becomes indistinct. This occurs at the critical point shown in Fig. 7.7. Beyond the critical point, fluids become "supercritical" and no phase interface is visible. For example, the vapor at point "A" in Fig. 7.7 would become more liquid-like as the pressure and temperature are varied along the semicircle to point "B," a liquid. The fluid would never exhibit two phases during the process, as identified by an interface. Because there is no actual boundary that defines the supercritical region, fluids should be described in terms of liquid-like or vapor-like.

One important characteristic of critical fluids is that thermodynamic properties approach each other; that is, the densities/viscosities of the vapor and liquid phases become the same at the critical point. We often inject supercritical fluids, such as carbon dioxide, into reservoirs so that, as the supercritical fluid and the reservoir fluid mix, phase interfaces disappear, facilitating production of hydrocarbon components previously not recovered (see the chapter on steam injection in the Reservoir Engineering and Petrophysics section of this Handbook). Thus, an accurate characterization of critical behavior is very important.

A pressure-temperature diagram contains no information about molar volume—just phase boundaries. Consider an experiment in a closed system that is held at constant temperature (ignoring the solid phase). The experiment is designed to measure pressure as the volume is changed and is generally referred to as a PVT experiment. Fig. 7.8 illustrates several phase behavior states that occur in a PVT experiment. The chapter on phase diagrams in this section of the Handbook outlines similar experiments for multicomponent mixtures. The steps in the PVT experiment for a pure fluid are listed next.
  1. Place a known amount of a single-component fluid in a constant temperature PVT cell at vapor state one, illustrated in Fig. 7.8. The molar volume of the fluid can be calculated from RTENOTITLE. Measure the pressure.
  2. Inject a volume of mercury into the cell. Calculate VV again, and monitor the corresponding pressure change. Because vapor is highly compressible and the mercury is highly incompressible, the pressure change with decreasing molar volume of the vapor phase will be relatively small. The fluid remains a vapor.
  3. Continue to inject mercury until the vapor becomes saturated at state two. At this pressure (the vapor pressure, RTENOTITLE, in Fig. 7.7), a small amount of liquid forms (for multicomponent mixtures, this is the dew point). Calculate the saturated VV, and measure the pressure.
  4. Continue to inject mercury to state three, where the molar fractions of vapor and liquid are about equal. The pressure is constant at the vapor pressure RTENOTITLE, corresponding to temperature T1. The pressure is constant because the temperature is constant (i.e., the degrees of freedom from the Gibbs phase rule are one for two phases in equilibrium; thus, the pressure is fixed if the temperature is constant). Because the pressure does not change with decreasing molar volume, the compressibility of the closed system is infinite.
  5. Continue to inject mercury to state four until only a small drop of vapor remains (for multicomponent mixtures, this is the bubble point). Calculate the saturated liquid molar volume, VL.
  6. Continue mercury injection to state five and calculate VL, while monitoring the pressure. During this step, the fluid remains a liquid, and the pressure rises very quickly for small changes in the molar volume. The pressure rises quickly with decreasing molar volume because liquids are nearly incompressible compared to vapors.

The PVT process just described generates an isotherm on a pressure-volume diagram, as illustrated in Fig. 7.9. States one through five correspond to those in the isothermal PVT experiment. A vapor pressure dome or two-phase envelope outlines the two-phase region. Along the left boundary of the dome, the liquid is completely saturated, for any reduction in pressure would cause a bubble of vapor to form. Along the right boundary of the dome, the vapor is saturated. To the left of the dome the liquid is subcooled, whereas to the right of the dome, the vapor is superheated. Two metastable states have been observed within the two-phase region: superheated liquid and subcooled vapor. These states are unlikely to exist in a reservoir where nucleation is always present (see Firoozabadi[2] on nucleation).

From the Gibbs phase rule, the pressure within the two-phase region along the T1 isotherm remains at its corresponding vapor pressure as molar volume changes from states two to four (i.e., VV to VL). VV in Fig. 7.9 is the molar volume of the saturated vapor, whereas VL is the saturated liquid molar volume. VC is the molar volume of the fluid at its critical point.

The critical point, in which the molar volumes of the liquid and vapor become equal, is at the apex of the dome in Fig. 7.9. The critical point must be at the apex because the isotherms in the two-phase region are horizontal, as required by the Gibbs phase rule. Thus, the slope and the inflection point of the isotherm at the critical temperature must be zero. For isotherms where T > TC, the molar volume (inverse of molar density) changes continuously from vapor-like values to liquid-like values as the pressure increases.

The isothermal compressibility of the fluid is given by the inverse of the slopes of the isotherms. For isotherms where T < TC, the slope is small in the superheated vapor region, indicating that the fluid compressibility is large. The slope is large in the subcooled liquid region, indicating that the fluid is nearly incompressible. At the boundary of the two-phase dome, the isothermal compressibility is discontinuous and is equal to infinity (zero slope). Compressibility is infinite in the two-phase region because, as the system volume is decreased, some of the vapor molecules, which occupy more space, are condensed into the denser liquid phase in which they occupy less space.

For a two-phase mixture at constant temperature, such as that shown by state three in Fig. 7.9, the molar phase volumes, VV and VL, do not change as the volume of the closed system changes. This occurs because the pressure is constant in the two-phase region for a pure fluid. Thus, as the volume is changed, the molar densities of the vapor and liquid phases do not change, but only the molar fractions (or amounts) of the phases change. The overall density of the two-phase mixture will change as the closed system is compressed or expanded.

For example, consider state three in the two-phase region at a temperature of T1. The vapor molar density is RTENOTITLE (at state two), whereas the liquid molar density is RTENOTITLE (at state four). The mole fraction of vapor at state two is V. Extensive parameters, such as the total molar volume of the two-phase mixture, can be calculated by


For total molar enthalpy,


where HL is the molar enthalpy of the saturated liquid, and HV is the molar enthalpy of the saturated vapor. Eqs. 7.36 and 7.37 can be solved for V as


Eq. 7.38 is known as a lever rule. The fluid quality is the molar volume fraction V as a percentage. An illustration of quality lines is given in Fig. 7.10. A quality of 100% is a saturated vapor. State three in Fig. 7.10 has a quality of about 60%.

Phase Behavior Models of Pure Fluids

An accurate characterization of phase behavior is critical to the prediction of oil recovery. Often, sufficient PVT experimental data is not available, and mathematical models that are "tuned" to experimental data are needed. EOS calculations are used for this purpose. EOS models are typically easy to implement in a numerical simulator.

One of the first EOS was Boyle’s and Charles’s law. These laws were combined into the ideal gas equation we use today, pV = RT. The ideal gas equation is generally satisfactory for vapors at pressures below a few atmospheres. Numerous other types of EOS were developed through the years,[5] but the most widely used EOS type in the petroleum industry is the cubic EOS.

Van der Waals developed the first cubic EOS. Unlike the ideal gas equation, which is limited to low pressure vapors, the van der Waals EOS attempted to provide good phase behavior estimates for both liquids and vapors by using only one equation. He also developed the principle of corresponding states, which is frequently used in the petroleum industry today. Numerous cubic EOS models are available today that give better accuracy than the van der Waals’ EOS. The two most widely used cubic EOS are the Peng-Robinson EOS[8] and modified versions of the Redlich-Kwong EOS.[9][10]

Prediction of the phase behavior of real reservoir fluids is difficult because of the complex interaction of molecules. Intermolecular forces of attraction and repulsion determine thermodynamic properties for any mixture of molecules. The attraction forces allow fluids to form liquid and solid phases, whereas repulsions are responsible for resistance to compression.

The accuracy of any EOS depends on its ability to model the attractions and repulsions between molecules over a wide range of temperatures and pressures. EOS models are empirical in that they do not attempt to model the detailed physics but only the cumulative effect in terms of a small number of empirical parameters. Generally, EOS models are more accurate when attractions are small, which explains why water, a polar substance, is more difficult to model.

Ideal Gas Equation. The simplest and most fundamental EOS is the ideal gas equation, in which the pressure, volume, and temperature of a fluid are related by


As stated previously, the behavior of a gas may be approximated by Eq. 7.39 if the pressure is relatively low (i.e., less than a few atmospheres). A gas is ideal if molecular interactions are negligible, something that could only occur at zero pressure. Thus, molecular interactions are negligible at zero pressure; therefore, thermodynamic properties, such as the molar internal energy of an ideal gas, are only a function of temperature.

Real Fluid Equation. Most fluids are not ideal. Real fluids, whether vapor or liquid, can be defined by the compressibility factor


Intensive thermodynamic properties, such as the molar internal energy, are a function of both temperature and pressure for a real fluid. A common mistake is to believe that Eq. 7.40 is only applicable to gases. Eq. 7.40 will be used later in this chapter to represent the behavior of any phase, whether liquid or vapor.

Coefficient of Isothermal and Isobaric Compressibility. For a single-phase fluid, we often employ very simple equations-of-state that describe the relationships between pressure, temperature, and molar volume. As stated previously, molar volume is a state function—it is not determined by the process or path taken to that state. Thus, for every temperature and pressure, there corresponds only one value of molar volume, and hence, we can write an equation that describes differential changes in molar volume as


Division of Eq. 7.41 by the molar volume gives


where RTENOTITLE is the isothermal compressibility of the fluid, and RTENOTITLE is the isobaric compressibility of the fluid. The minus sign is introduced, per convention, to make the compressibilities positive.

Eq. 7.42 describes the fractional change in the molar volume for small changes in temperature and pressure. Because reservoirs are usually thought to be isothermal, the isothermal compressibility is most often used in reservoir engineering. It is often written in terms of density or formation volume factor (FVF) as


Several special cases of Eq. 7.42 are useful to consider. First, when the volume under consideration is closed (dV = 0), Eq. 7.42 gives RTENOTITLE. Under the assumption of constant compressibilities, this result can be integrated to RTENOTITLE, which gives the pressure change in a closed volume resulting from a temperature change.

Second, when the process is isothermal, Eq. 7.42 reduces to RTENOTITLE. This result can be integrated under the assumption of constant compressibility to obtain


where subscript o is a reference state. Eq. 7.43 is often referred to as the EOS for a constant compressibility fluid. It can be rewritten in a more familiar form as RTENOTITLE or RTENOTITLE. This result is often used in well testing.

Third, when the compressibility of the fluid is both constant and small, and the pressure change is not too large, the exponential term in Eq. 7.43 can be simplified so that V = Vo(1 - cΔp) or ρ = ρo(1 + cΔp) . This result is often referred to as the EOS for a fluid with a slight, but constant, compressibility.

Last, substitution of the definition of compressibility factor (Eq. 7.40) into Eq. 7.42 gives


Eq. 7.44 is exact, whether the fluid is a vapor or liquid. For an ideal gas, the compressibility factor is constant and equal to 1.0, and Eq. 7.44 reduces to RTENOTITLE. Thus, the isothermal compressibility for an ideal gas is inversely proportional to pressure.

Cubic Equations-of-State. The van der Waals EOS was the first EOS capable of representing both the liquid and vapor. The van der Waals EOS, however, is not used because of its lack of accuracy near critical points. The Peng-Robinson EOS and a modified version of the Redlich-Kwong EOS are generally used (see Table 7.2 for comparison). Nevertheless, the simplicity of the van der Waals EOS makes it useful in demonstrating several key concepts.

As for other cubic EOS, the van der Waals EOS describes the pressure as a function of molar volume and temperature, which is written as


where b is the repulsion parameter, and a is the attraction parameter. The first term on the right side of Eq. 7.45 attempts to represent the pressure deviation from an ideal gas that results from molecules occupying and competing for space. The effective volume available for movement of the molecules (on a molar basis) is Vb—not V, as it would be for an ideal gas. Thus, b represents the smallest possible volume that one mole of molecules can occupy (no space would exist between the molecules).

As V approaches b, the denominator in the first term on the right side becomes small so that pressure increases very rapidly. Because b is based on the effective molecule size, the value for b will change with the nature of the pure fluid and will determine the lower boundary for the region of interest on a pressure-volume diagram (i.e., the physical region of interest is only where V > b).

The second term on the right side of Eq. 7.45 accounts for the attractive forces between molecules. The attractive forces will be proportional to the square of the number of molecules present and, thus, RTENOTITLE on a macroscopic scale. The proportionality constant, a, depends on the nature and strength of the forces between the molecules and, therefore, the fluid type. As V becomes large, the contribution of the attractive forces becomes small, and the second term on the right side of Eq. 7.45 becomes negligible. Under these conditions, the van der Waals EOS approaches the ideal gas equation (Eq. 7.39). This is also true for the other EOS in Table 7.2.

The van der Waals EOS can be rewritten in cubic form as


or in terms of the compressibility factor,


Determination of a and b from van der Waals EOS. We previously demonstrated that an isotherm on a pressure-volume diagram must have zero slope at the critical point. There is also an inflection point there. Thus, any cubic EOS must be constrained by


Eq. 7.47 consists of two equations that can be solved simultaneously to determine the unknowns, a and b. This procedure can be somewhat tedious, especially for more advanced cubic EOS.

A simpler method is to recognize that Eq. 7.47 implies that only one volume root of the cubic EOS exists at the critical point. Mathematically, this can be expressed as (V - VC)3 = 0, which upon expansion is RTENOTITLE. Comparison of this expansion to Eq. 7.48 term by term at the critical point gives RTENOTITLE, RTENOTITLE, and RTENOTITLE. These three relationships are easily solved to obtain a and b in terms of only the critical pressure and temperature, as well as to determine the relationship between critical parameters. That is, RTENOTITLE; RTENOTITLE; and RTENOTITLE. A similar procedure can be used for any cubic EOS.

The result shows that the critical compressibility factor is constant and is independent of the fluid. This is true for the other cubic EOS models in Table 7.2. For the van der Waals EOS, we obtained RTENOTITLE. In reality, the critical compressibility factor is not constant for different fluids and is generally smaller than 0.3. For example, water has a critical compressibility of about 0.23, and that of carbon dioxide is 0.27. The poor match of the van der Waals EOS at the critical point explains why it is not used today. Both the Redlich-Kwong EOS and the Peng-Robinson EOS have critical compressibilities closer to measured values (see Table 7.2).

A constant critical compressibility factor means that only two of the three critical properties can be satisfied at the critical point. For example, if critical pressure and temperature are specified, the critical volume will not be correctly predicted. For the van der Waals EOS, RTENOTITLE. Thus, when the critical pressure and temperature are specified, the critical volume or critical density of the fluid is likely in error. In general, liquid densities predicted by cubic EOS exhibit greater error than do vapor densities. Firoozabadi[2] gives a good description of how density predictions can be improved using volume translation parameters.

We determined a and b only at the critical point. Their values could be different away from the critical point and could be functions of temperature. Therefore, in general, RTENOTITLE, and RTENOTITLE, where α and β are functions of temperature; these functions must approach 1.0 at the critical point. For most cubic EOS, β is taken to be 1.0, and α is adjusted to give the correct vapor pressure (see Table 7.2).

Construction of Pressure-Volume Diagram From Cubic EOS. Once a and b are defined, a cubic EOS can be used to generate the pressure-volume diagram, as illustrated in Fig. 7.9. From the Gibbs phase rule, there is only one degree of freedom within the two-phase region for pure fluids. We assume that temperature has been specified but not the pressure. Our goal, therefore, is to solve for the vapor pressure, given temperature and the critical properties of the fluid.

The procedure used to estimate a and b ensures that the cubic EOS gives the experimentally measured shapes of the isotherms for temperatures greater than or equal to the critical point. Fig. 7.11 illustrates the isotherms generated from a cubic EOS. As shown by the loop in the curve, cubic EOS can have three roots for isotherms below the critical temperature. The loop is not physical because the pressure must be constant in the two-phase region.

The nonphysical condition must be removed to achieve the correct physical response in the two-phase region. That is, the loop within the two-phase region must be discarded and replaced by the correct vapor pressure. The procedure is relatively simple. Within the two-phase region, there are three roots along an isotherm. At the vapor pressure, the root with the largest molar volume is taken to be the molar volume of the saturated vapor VV, whereas the smallest root is VL. The middle root is discarded because choosing that root would lead to unstable phases (see Sec. 7.2.2). The middle root is clearly nonphysical in that it is located on the isotherm where RTENOTITLE (i.e., pressure increases as density decreases).

Principle of Corresponding States. Correlations for reservoir fluids, such as the generalized compressibility factor charts for natural gases, use reduced temperature, pressure, and volume, where RTENOTITLE, RTENOTITLE, and RTENOTITLE. A cubic EOS shows why these parameters give good correlations. For example, substitution of the reduced parameters into the van der Waals EOS (Eq. 7.45), along with the definitions of a and b, gives


Eq. 7.48 is dimensionless and is often called the reduced form of the van der Waals EOS. The reduced form leads directly to the principle of corresponding states. The two-parameter principle of corresponding states says that all fluids, when compared at the same reduced temperature and pressure, have approximately the same compressibility factor, and all deviate from ideal-gas behavior by about the same degree. The reduced compressibility factor is given by RTENOTITLE. Because ZC is constant for a cubic EOS, the compressibility factor is constant for the same reduced temperature and pressure (reduced volume is related to reduced pressure and temperature through Eq. 7.48).

The principle of corresponding states is a powerful idea even though it is only qualitatively correct. Experiments show that ZC is not constant for different fluids. Nevertheless, it demonstrates that, to obtain reasonable estimates of fluid properties, only the reduced pressure and temperature must be known. This is why most fluid correlations use reduced temperature and pressure.

In reality, fluids can deviate from the principle of corresponding states. Pitzer noted that a plot of RTENOTITLE vs. RTENOTITLE for simple fluids (molecules that are roughly spherical in shape such as the Noble gases) collapse onto a straight line (see Fig. 7.12). The parameter, RTENOTITLE, is the reduced pressure at the vapor pressure. Other more complex and nonspherical molecules such as hydrocarbon-chained molecules, however, do not plot on that same line. Thus, Pitzer defined an additional correlation variable called the acentric factor, where


The acentric factor measures the deviation of complex fluids from the simple fluids at a reduced temperature of 0.70 (see Fig. 7.12). Hydrocarbons with longer chains generally have greater acentric factors. For example, methane has an acentric factor of 0.008, while n-butane has an acentric factor of 0.193.

Because the acentric factor is simple to measure, it is often used to improve phase-behavior calculations from cubic EOS. The three-parameter principle of corresponding states is that a fluid will have about the same compressibility factor as another fluid, if the reduced pressure, reduced temperature, and acentric factors are similar.

Calculation of Vapor Pressure. Although the shape of an isotherm from a cubic EOS can be made qualitatively correct, the problem remains that the vapor pressure is unknown for a given temperature. The vapor pressure is determined using the equilibrium criterion of Eq. 7.24. For example, substitution of the Soave Redlich-Kwong EOS into Eq. 7.25 and subsequent integration gives the fugacity as a function of pressure, molar volume, and temperature. That is,


The fugacity of the vapor is computed using the molar volume of the vapor phase, VV, whereas the liquid fugacity is determined using VL. Thus, the fugacity for the vapor phase from the Soave Redlich-Kwong EOS is


where RTENOTITLE; and for the liquid phase,



The problem of calculating the vapor pressure reduces to finding the pressure that gives the required phase molar volumes so that the fugacities of the phases are equal. Fig. 7.13 illustrates the procedure. The procedure works well as long as the initial guess for the pressure is in the range of the cubic EOS where three roots exist (i.e., the pressure is within the loop of the cubic EOS). If the pressure is above the critical pressure, only one root exists for the molar volume. This is also true if the pressure is below the minimum value of the loop (the minimum pressure of the loop could be negative).

For a pure fluid, the vapor pressure can also be determined graphically with Maxwell’s equal area construction. Fig. 7.11 shows that the vapor pressure is the pressure required so that the areas bounded by the vapor pressure line and the loop of the cubic EOS must be equal. The equal area construction results from the equality of Gibbs energy (or fugacities). This method is less accurate but serves as a useful check to the calculated vapor pressure. Firoozabadi[2] outlines Maxwell’s construction method in detail for pure fluids and mixtures.

Example Calculation of Two-Phase Envelope. This section demonstrates a calculation of vapor pressure and the two-phase envelope for a pure fluid using the procedure outlined in Fig. 7.13. Propane is selected as the pure fluid at a temperature of 40°C (313°K). We also select the Soave-Redlich-Kwong EOS to model the phase behavior. The properties for propane are a critical temperature of 370°K; a critical pressure of 42.5 bars; and an acentric factor of 0.152. The gas constant in consistent units is 83.1 cm3-bar/mol/K.

With these values, the parameter "a" in Table 7.2 is calculated to be 9.51E6 cm6-bar/mol2, and parameter "b" is 62.7 cm3/mol. The value for α from Table 7.2 is found to be 1.05 (the reduced temperature at 40°C is 0.913). Fig. 7.14 shows the isotherm at 40°C calculated with the SRKEOS.

Based on the calculated isotherm, we select an initial value of 10 bars for the vapor pressure. Any value within the S-loop of the isotherm would be satisfactory as an initial guess for the vapor pressure. A vapor pressure of 13.8 bars is calculated with the iteration procedure of Fig. 7.13. The calculated vapor pressure is the pressure at which the fugacities of the vapor and liquid phases are equal (illustrated in Figs. 7.11 and 7.14 by the Maxwell equal-area rule). The fugacities are 11.3 bars at the vapor pressure, which are calculated with Eqs. 7.49 and 7.50. The equilibrium liquid molar volume is 105 cm3/mol, and the vapor molar volume is 1462 cm3/mol.

Fig. 7.14 also illustrates the phase behavior with the SRKEOS at a higher temperature of 70°C. The calculated vapor pressure at this temperature is 26.2 bars. The equilibrium liquid molar volume is 128 cm3/mol, and the vapor molar volume is 691 cm3/mol. Thus, as the temperature is increased, the size of the two-phase region shrinks. Fig. 7.14 shows the two-phase envelope that connects the liquid and vapor molar volumes. At the critical temperature, the two-phase region disappears.

The values for vapor pressure and molar volumes are calculated parameters only. Using the critical temperature and pressure, the critical volume from the SRKEOS is approximately 241 cm3/mol (ZC = 1/3 always for the SRKEOS). The actual critical volume from experimental data is 200 cm3/mol, which is about 20 percent smaller than the calculated value. Firoozabadi[2] outlines a more complex method to improve the calculated match of cubic EOS to experimental data.

Volumetric Properties of Mixtures

The interaction of the different molecules in a mixture causes behavior not observed in pure fluids. The chapter on phase diagrams in this section of the handbook describes the volumetric behavior of mixtures. Sec. 7.5 presents the basic procedure to predict the equilibrium phase behavior of mixtures by a cubic EOS. More detailed information can be found in many sources, including Firoozabadi[2] and Whitson and Brule. [7]

The thermodynamic properties of a mixture can be calculated with the same EOS for a pure fluid, with some modifications. The primary difference is that the mixture molar volume for a phase is calculated with EOS constants and temperature-dependent functions of the phase molar composition, either xi or yi. For example, the Soave Redlich-Kwong EOS written for a mixture is


where subscript m indicates a mixture property. The mixture properties are calculated with mixing rules that are often linear or quadratic functions of the phase mole fractions. For example, for the liquid phase, the mixing rule for the product, , is often the quadratic equation,


where RTENOTITLE. The parameters kij are called binary interaction parameters. Binary interaction parameters are constants that are determined by fitting the cubic EOS to experimental PVT data. The mixing rule for is theoretically justified from virial EOS, which are discussed in several sources.[2][3][4][5][6][7] For bm, the linear relationship, RTENOTITLE, is often used.

For equilibrium calculations, the fugacity of every component in each phase must be calculated. Eq. 7.35 is used for this purpose. For example, substitution of the Soave Redlich-Kwong EOS into Eq. 7.35 gives the fugacity of a component in the liquid phase, which is written as


where RTENOTITLE. A similar equation is written for the vapor phase, where xi is replaced by yi, and superscript L is replaced by V.

Procedure for Equilibrium Calculations of a Mixture

The procedure for equilibrium calculations of a potential two-phase mixture is more complex than that of a pure fluid. For an equilibrium flash calculation, the pressure and temperature and overall mole fractions are specified (i.e., pressure and temperature are now independent, as specified by the Gibbs phase rule). The general procedure for a flash calculation is discussed next.

  1. Make an initial guess of the K-values, where RTENOTITLE. When the guess of the K-values is near the equilibrium solution, the procedure will converge rapidly. If the guess is not good, the procedure might not converge at all. Most EOS programs use some empirical correlation to estimate the phase mole fractions based on K-values. The Wilson equation[11] is often used, where RTENOTITLE.
  2. Calculate xi and yi with the Rachford-Rice procedure.[12] Once the K-values for each component are specified, the Rachford-Rice procedure is used to estimate the phase mole fractions. A material balance on each component gives zi = Lxi + (1 - L)yi, where L is the mole fraction liquid (see Eq. 7.1). Substitution of yi = Kixi into the material balance equation gives, upon rearrangement, RTENOTITLE or alternatively RTENOTITLE. Substitution of these expressions into the function RTENOTITLE gives RTENOTITLE. This is a nonlinear equation that can best be solved by a simple Newton-Raphson iteration, where for each iteration the new value of the liquid mole fraction is found from RTENOTITLE. For the first iteration, choose L = 0.5 and iterate until RTENOTITLE.
  3. Calculate the cubic EOS parameters (e.g., am and bm). This step is very straightforward and depends on the selected EOS and its associated mixing rules. The critical temperatures, pressures, and acentric factors for each component are needed to calculate the EOS parameters.
  4. Solve the cubic EOS for the phase molar volumes VmL and VmV. This step requires solution of the cubic EOS for the compressibility factor, Z, of the vapor and liquid (or alternatively for VV and VL). Because the compositions of the vapor and liquid are different, two separate solutions for the roots of the cubic EOS are required. A cubic equation-solver or iteration method should be used to obtain the roots of the EOS.

The procedure for this step is more complex than for a pure fluid because six roots of the cubic EOS are calculated (i.e., three roots for the liquid and three for the vapor). The middle root for the vapor and liquid are discarded because that solution leads to unstable phases, similar to pure fluids. One of the remaining two liquid roots is paired with one of the other vapor roots to calculate component fugacities and equilibrium. If the wrong root pairing is selected, the solution could be false in that an unstable or metastable solution could be obtained. The correct equilibrium solution is the one that minimizes the total Gibbs energy compared with the other possible root-pairings. Two sources[2][7] provide a good description of how to select the liquid and vapor roots so that the total Gibbs energy of the two-phase mixture is minimized. For most cases, the correct root for the liquid is the one that gives the smallest molar volume, and the correct root for the vapor is the one that gives the largest molar volume.

A couple of sources[2][7] also examine using stability analyses to determine whether a mixture will form three phases instead of just one or two phases. The chapter on phase diagrams in this section of the handbook discusses the formation of three equilibrium phases in CO2/crude oil systems.

  1. Calculate the component fugacities of each component in each phase, RTENOTITLE and RTENOTITLE. The selected cubic EOS is used to determine an expression for the fugacity of a component in a phase (see Eq. 7.51 for example).
  2. Check to see if equilibrium has been reached. A good criterion is RTENOTITLE for all components. If the criteria are satisfied, equilibrium has been obtained. The correct equilibrium solution is found when RTENOTITLE for each component. Because the solution is never found exactly, we accept the solution if RTENOTITLE for each component. The tolerance of 10-5 can be decreased if better accuracy is required.
  3. If the criteria have not been satisfied, the K-values should be updated and steps two through six repeated. This step is also very important; it affects both the rate of convergence and whether the iteration converges at all. One procedure that works well is the simple successive substitution scheme that relies on the fact that RTENOTITLE and RTENOTITLE for each component. Therefore, RTENOTITLE. At equilibrium, the component fugacities are equal so that RTENOTITLE. We can use this ratio to estimate new K-values from the old ones. That is, RTENOTITLE. Once the new K-values are determined, steps two through six are repeated until convergence in step six is achieved. Convergence from successive substitutions can be slow near the critical region. Other methods may be required when convergence is slow.[2]

Characterization of In-Situ Fluids

Phase behavior calculations require that all components and their properties be specified. Crude oils, however, typically have hundreds of components, making the EOS procedure in Sec. 7.5 computationally intensive. Thus, components are often lumped into pseudocomponents to approximate the in-situ fluid characterization. The characterization usually takes the following three steps:

  1. The hydrocarbon components in the in-situ fluid are analyzed using analytical techniques, such as chromatography or distillation. New analytical techniques often give a reliable analysis for hydrocarbon components up to C30, instead of the traditional C7. Properties for hydrocarbon components greater than C30 are reported as a C30+ fraction.
  2. The measured components are separated and lumped into a minimum number of pseudocomponents. The chosen number of pseudocomponents is often a result of the measured fluid characterization and the degree of accuracy required (see step three). The properties and selection of the pseudocomponents are determined using a variety of methods as reported in Whitson and Brule.[7] The required pseudocomponent properties are those needed for the cubic EOS calculations, such as critical temperature, pressure, and acentric factor.
  3. The pseudocomponent properties are adjusted to match all available phase behavior data (e.g., PVT reports). This process, which typically uses nonlinear regression, is known as EOS tuning. EOS tuning is needed because of the inherent uncertainty in the properties estimated from step two, especially for the heavier components. Binary interaction parameters are typically the first parameters to be adjusted, although all of the parameters may need some tuning. The number of pseudocomponents may need to be increased from step two to obtain a good fit of the calculated phase behavior with the measured phase behavior data.

The selection of pseudocomponents and their property values are likely not unique, as is often the case when numerous model parameters are estimated by fitting measured data with nonlinear regression. Care should be taken to avoid estimates in the pseudocomponent properties that are unphysical and to reduce the number of parameters. Furthermore, the final EOS characterization is most accurate in the range of the measured phase behavior data. Phase behavior data should be collected that covers, as much as possible, the conditions that occur in the reservoir. The characterization should be updated when new data becomes available.

Finally, fluid characterizations may vary from one location in the reservoir to another. In such cases, multiple EOS characterizations might be required. Compositional variations can occur for a variety of reasons. For example, gravity can cause vertical compositional gradients, where heavier components become more concentrated at greater depths. Several sources[2][13][14] provide examples of gravitational concentration gradients. Variations caused by thermal gradients are also discussed in Firoozabadi.[2]


a = constant parameter in cubic EOS in Table 7.2, pressure-volume2, Pa-m6
A = area normal to specified direction, m2
b = constant parameter in cubic EOS in Table 7.2, volume/mole, m3/mole
B = formation volume factor of the fluid, volume/volume
c = isothermal compressibility of a fluid, 1/pressure, 1/Pa
cp = isobaric compressibility of a fluid, 1/pressure, 1/Pa
f = fugacity of a pure fluid, pressure, Pa
RTENOTITLE = fugacity of a component in a mixture, mole2-pressure/mole2, Pa
F = external force on one side of system, energy/length, J/m
RTENOTITLE = external force vector of surroundings on system, energy/length, J/m
G = molar Gibbs free energy, energy/mole, J/mole
h = heat transfer coefficient, energy/temperature/time, J/(Kelvin-sec)
H = molar enthalpy of fluid, energy/mole, J/mole
k = binary interaction parameter, dimensionless
Ki = K-value of ith component, yi/xi, dimensionless
RTENOTITLE = displacement vector of system, length, m
L = liquid mole fraction, moles liquid/total moles, dimensionless
M = net mass transferred, mass, moles
n = total moles of all components, moles
nc = number of components
np = number of phases
p = pressure, force/area, Pa
Q = net heat transferred, energy, J
R = gas constant, pressure-volume/temperature/mole, Pa-m3/(Kelvin-mole)
S = molar entropy of fluid, entropy/mole, J/(Kelvin-mole)
t = time, seconds
T = temperature, Kelvin
u = velocity of fluid, length/time, m/sec
U = molar internal energy, energy/mole, J/mole
V = vapor mole fraction, moles vapor/total moles, dimensionless or molar volume of fluid, volume/mole, m3/mole
W = net work performed, energy, J
x = x-coordinate, length, m
xi = mole fraction of ith component in liquid, moles ith component in liquid/total moles liquid, dimensionless
y = y-coordinate, length, m
yi = mole fraction of ith component in vapor, moles ith component in vapor/total moles vapor, dimensionless
z = z-coordinate, length, m
zi = overall mole fraction of ith component, moles ith component/total moles, dimensionless
Z = compressibility factor of a fluid, dimensionless
α = temperature dependence function in cubic EOS in Table 7.2, dimensionless
β = temperature dependence function in cubic EOS, typically set to 1.0, dimensionless
μi = chemical potential of ith component, energy/mole, J/mole
ρ = molar density of fluid, moles/volume, mole/m3
φ = fugacity coefficient for a pure fluid, pressure/pressure, dimensionless
RTENOTITLE = fugacity coefficient for a component in a mixture, mole2-pressure/mole2-pressure, dimensionless
ω = acentric factor, dimensionless


A = open subsystem A
B = open subsystem B
C = state is at critical point
ext = external to system
G = generated quantity within system
i = ith component
j = jth component
L = liquid
m = mixture
o = reference state
rev = reversible process
R = reduced parameter, ratio of quantity/critical value, dimensionless
T = total
v = vapor pressure
V = vapor
x = direction is along x-coordinate
y = direction is along y-coordinate
z = direction is along z-coordinate


ig = ideal gas
L = liquid
v = vapor pressure
V = vapor


  1. 1.0 1.1 Gibbs, J.W. 1961. The Scientific Papers of J. Willard Gibbs. H.A. Bumstead and R.G. Van Name, eds. New York City: Dover.
  2. 2.00 2.01 2.02 2.03 2.04 2.05 2.06 2.07 2.08 2.09 2.10 2.11 Firoozabadi, A. 1999. Thermodynamics of Hydrocarbon Reservoirs. 355. New York City: McGraw-Hill Book Co. Inc.
  3. 3.0 3.1 Prausnitz, J.M., Lichtenthaler, R.N., and de Azevedo, E.G. 1999. Molecular Thermodynamics of Fluid-Phase Equilibria, third edition. New Jersey: Prentice Hall.
  4. 4.0 4.1 4.2 Sandler, S.I. 2000. Chemical and Engineering Thermodynamics, third edition. New York City: John Wiley & Sons.
  5. 5.0 5.1 5.2 5.3 5.4 Smith, J.M., Van Ness, H.C., and Abbott, M.M. 2001. Chemical Engineering Thermodynamics, sixth edition. 787. New York City: McGraw-Hill Book Co. Inc.
  6. 6.0 6.1 6.2 Walas, S.M. 1985. Phase Equilibria in Chemical Engineering. 671. Boston, Massachusetts: Butterworth Publishings.
  7. 7.0 7.1 7.2 7.3 7.4 7.5 Whitson, C.H., and Brule, M.R. 2000. Phase Behavior, Vol. 20. Richardson, Texas: Monograph Series, SPE.
  8. Peng, D. and Robinson, D.B. 1976. A New Two-Constant Equation-of-State. Industrial & Engineering Chemistry Fundamentals. 15 (1): 59.
  9. Redlich, O. and Kwong, J.N.S. 1949. On the Thermodynamics of Solutions, V. An Equation-of-State. Fugacities of Gaseous Solutions. Chemical Reviews 44 (1): 233.
  10. Soave, G. 1972. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Sci. 27 (6): 1197-1203.
  11. Wilson, G.M. 1969. A Modified Redlich-Kwong Equation-of-State, Application to General Physical Data Calculations. Paper 15c presented at the AIChE Natl. Meeting, Cleveland, Ohio, 4–7 May.
  12. Rachford Jr., H.H. and Rice, J.D. 1952. Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium. J Pet Technol 4 (10): 19, 3. SPE-952327.
  13. Sage, B.H. and Lacey, W.N. 1939. Gravitational Concentration Gradients in Static Columns of Hydrocarbon Fluids. Transactions of the AIME 132 (1): 120-131. SPE-939120-G.
  14. Schulte, A.M. 1980. Compositional Variations Within a Hydrocarbon Column Due to Gravity. Presented at the SPE Annual Technical Conference and Exhibition, Dallas, Texas, 21–24 September. SPE-9235-MS.

SI Metric Conversion Factors

bar × 1.0* E + 05 = Pa
ft × 3.048* E – 01 = m
ft2 × 9.290 304* E – 02 = m2
ft/sec × 3.048 E – 01 = m/s
°F (°F – 32)/1.8 = °C
kelvin (K – 273.15) = °C


Conversion factor is exact.