pvt

PVT Analysis

conservation laws. If we have a mixture of N-components, each of which has ni moles, then the total moles is:

	N
(5.7)	nT ni

i 1

We can now define the mole fraction for the ith component as:

n

(5.8) zi ni

T

If the mole weight of each species in our mixture is Mwi, then the mixture mole weight is:

	N
(5.9)	M wm zi M wi

i 1

5.1.2Deficiencies in the Ideal Gas Law

The two principle deficiencies in the Ideal Gas Law are:

Prediction of non-physical zero volume at zero temperature, and

No account of second phase.

The problem of zero volume is easily corrected. The ideal gas law has an implicit assumption that the molecules occupy zero volume: instead, we will assume they occupy a [molar] volume of b, namely:

(5.10) Vm Vm b

This so-called hard-sphere approximation has the effect of defining the maximum packing possible as the fluid pressure is raised infinitely.

The presence of a second [liquid] phase is handled by adding another term to the pressure to account for attractive forces between molecules. van der Waals proposed:

	(5.11)	p p	a
			V 2
			m

Where a is a second constant which will be determined shortly.

5.1.3 The Real Gas Law

The Ideal Gas law now becomes the Real Gas Law, where:

Roxar Oxford

61

12/12/12

PVT Analysis

(5.12)

pV ZRT

The dimensionless Z-factor encapsulates the departure from non-ideal behaviour. We saw in section 3.4 that the Z-factor of simple gases and mixtures of gases can be predicted from charts or correlations fitted to those charts. For mixtures involving heavier hydrocarbons, we will resort to Cubic equations to find the deviation factor, Z.

5.2 Cubic EoS

With the two corrections proposed, we have the van der Waals Equation of State, first discussed in 1873:

			a
	(5.13)
		p V 2		Vm b RT
			m

This can be re-arranged to give:

(5.14)	p	RT		a
		Vm b		Vm2

All the popular EoS used in petroleum engineering calculations are modifications of this equation. We will study the essential features of this equation before we consider the modern EoS. Søreide gives an excellent account of the development of cubic EoS.

5.2.1 Van der Waals EoS

Equation (5.14) can be expanded in volume to give:

(5.15)	V 3	b	RT V 2			a	V ab 0
			p			p	p

If we define the following three terms:

(5.16)	A	ap	B	bp	Z	pV
		RT 2		RT		RT

Now (5.15) can be re-written as:

(5.17)

Z 3 B 1 Z 2 AZ AB 0

Let us now look at again at the shape of the p-V curve for a pure component.

Roxar Oxford

62

12/12/12

PVT Analysis

Figure 37: p-V Behaviour for Pure Component with Cubic EoS Behaviour.

Three isotherms [constant temperature lines] are shown at temperatures T1 < Tc < T2. The highest temperature line at T2, IJ, is characteristic of ideal gas behaviour in that pV constant: at any given pressure, the fluid can have only 1 unique volume solution.

The lowest temperature line at T1, ABDE has three characteristic parts to it. Section AB is typical of liquid behaviour: large pressure changes give corresponding small changes in volume. Section DE is typical of vapour behaviour: small pressure changes give corresponding large changes in volume. Line BD corresponds to a point on the VLE line. Point B is the liquid volume VL and point D is the vapour volume, VV. The cubic EoS approximates the true behaviour in the 2-phase region, shaded, by predicting three real

roots at B, D and C, however, the root at C is unphysical since dpdV 0 .

As the temperature is increased from T1, the points B and D come together [along with the spurious point C] until at T = Tc, there are three real equal roots at point G – the Critical Point.

This condition of three real equal roots can be written in mathematical form as:

(5.18)	f Z Z Zc 3	0
This can be expanded as:
(5.19)	f Z 3 3Zc Z 2	3Zc2 Z Zc3	0

Comparing the coefficients of (5.17) and (5.19), we see:

(5.20)	3Z	c	B	c	1	3Z 2	A	Z 3	A B	c
						c	c	c	c

Roxar Oxford

63

12/12/12

12/12/12

versus temperature for each

PVT Analysis

Note we have applied (5.17) at the critical point, hence the use of the subscript ‘c’. Simple algebra then shows that:

(5.21)	AvdW	27	BvdW	1	Z vdW	3
	c	64	c	8	c	8

Substituting the first two item of (5.21) into (5.20) we have:

(5.22)				RT 2				RT
	ac A pc						bc B pc
				c				c
For the van der Waals EoS:
(5.23)	vdWA		27		vdWB		1
			64				8

As we shall see, the values of these magic numbers – the Omega-A and Omega-B – along with the critical Z-factor depend on the form of the EoS.

5.2.2 Redlich-Kwong Family of EoS

Modifying the pressure correction (5.11) generates the RK family of EoS:

	(5.24)	p	RT		a
			(V b)		V (V b)

This gives rise to the following equation in Z:

(5.25)

Z 3 Z 2 A B B2 Z AB 0

Comparing coefficients at the critical point, it can be shown that:

	RK		1	RK		21 3 1	RK		1
(5.26)	A		9 21 3 1	B		3	Zc		3

In a further modification, RK changed the a-coefficient such that:

(5.27)

a ac (T )

The new term, , is a temperature dependent correction to the critical value, ac. For the original RK EoS:

(5.28)	RK	T 0.5
		r

5.2.2.1 Zudkevitch Joffe RK EoS

Rather than accepting the A , B as fixed constants, Zudkevitch and Joffe suggested that they could be functions of temperature, or more particularly, reduced temperature. Originally, this consisted of setting up a table of A , B

Roxar Oxford

64

PVT Analysis

phase, which were adjusted to match pure component saturated density and fugacity data. In a modification of their original proposal, they suggested just the A , B for the liquid phase be determined and these values be applied to the vapour phase also.

The model is often referred to as the ZJRK EoS.

5.2.2.2 Soave RK EoS

Soave improved the ability of the original RK EoS to predict pure component VLE behavior. This he achieved by making the parameter introduced in (5.23) not just a function of reduced temperature but of Acentric factor15, , also.

Soave conducted a series of experiments on light hydrocarbons at varying temperatures. By requiring 1 at Tr 1 he was able to generate the following expression:

(5.29)	0.5	1 m 1 Tr 0.5
The term, m,	is	to include the	dependence on . Soave proposed two possible

correlations for m( ). A more detailed study by Graboski and Daubert [see Søreide] suggested:

(5.30)

m 0.47979 1.576 0.1925 2 0.025 3

The Soave RK or SRK EoS has proved to be one of the two most successful EoS used in the upstream petroleum industry.

5.2.3 Peng-Robinson EoS

The PR EoS is other most successful EoS used in the upstream petroleum industry. The PR EoS is:

(5.31)	p	RT				a(T )
		(V b)			V (V	b) b(V b)
Again it can be shown that:
(5.32)	PRA	0.457235...				BPR 0.077796...		ZcPR 0.307401...

Note these are not simple values but must be determined from solving cubic equations. The form of the a-coefficient dependence is the same as for the SRK, namely (5.27) and (5.29) except (5.30) becomes:

(5.33)

m 0.379642 1.48503 0.164423 2 0.016666 3

15 The acentric factor measures the nonsphericity of a molecule. Generally, , increases with molecular weight.

Roxar Oxford

65

12/12/12

PVT Analysis

Theoretically, the PR is a better predictor of liquid volumes because its critical Z-factor is lower. In practice, the critical Z-factor of hydrocarbons has a maximum of about 0.29 for Methane and decreases with increasing mole weight: C20 has a Zc 0.20. We will return to the issue of estimating liquid volumes when we consider the 3-parameter correction.

5.2.4 The Martin’s 2-Parameter EoS

The SRK and PR EoS are so similar, a common form is often used due to Martin:

	(5.34)	p	RT		a(T )
			(V b)		(V m b)(V	m	b)
					1	2

The coefficients (m1, m2) are given by:

EoS	m1	m2
SRK	0	1
PR	1+ 2	1- 2

The form of the cubic in Z is:

(5.35)	Z 3 E			Z 2			E Z E	0	0
			2				1
where:
	E2 1 m1 m2 B 1
(5.36)	E m m				2		m m		B2 m m		B A
	1		1				1	2	1	2
	E	0	m m			B2 (B 1) AB
			1	2

Using this form, the fugacity coefficient for the ith component in an N-component mixture [see Appendix A.2.1], which plays a key role in the Flash [see section 6], becomes:

B

1

A

A

B

Z m

B

(5.37)

ln

ln Z B

i

Z 1

2

i

i

ln

2

m

m

A

i

B

2

B

B

Z m B

1

1

We will see shortly when we look at multi-component mixtures how the (A,B) for a mixture are computed from the component (Ai,Bi).

5.2.5 Other Cubic EoS

One of the obvious weaknesses of the PR and SRK EoS is they only have two degrees of freedom, i.e. (A,B), implying a fixed value for critical Z-factor. To rectify this shortcoming, a number of authors have added a third parameter whose value can be set by requiring that the theoretical Z-factor match that observed experimentally.

Roxar Oxford

66

12/12/12

PVT Analysis

The better known examples of this class of EoS are due to Usdin and McAuliffe (UM):

	(5.38)	p	RT			a
			(V b)		V (V d)
	and Schmidt and Wenzel (SW):
	(5.39)	p	RT			ac (T )
			(V b)		V 2	(1 3 )bV 3 b2

See Søreide for details.

Although the prediction of volumetric properties, in particular liquid densities, is improved with this case of EoS, the VLE behaviour between PR and SW is broadly similar. This improved behaviour is not really justified compared with the extra complexity of the equation, and hence computational cost in assembling the fugacity coefficient and its derivatives required in the Flash calculation. This is especially so when we consider Volume Translation in section 5.4.

5.3 Multi-Component Systems

The b-coefficient was introduced to account for the fact that molecules are not pointobjects. Clearly, a satisfactory mixing rule for the b [or B] coefficient is:

	N
(5.39)	b xi bi

i 1

where xi is the mole fraction of the ith component and:

RT

(5.40) bi bi p ci

ci

It is generally assumed the component bi [and ai] take the values predicted as if the fluid were a single component.

The a-coefficient was introduced to account for interactions between molecules. The traditional and most widely used mixing rule is:

	N	N
(5.41)	a xi x j aij
	i 1	j 1

The usual expression for the coefficient aij is:

(5.42)

aij ai a j 1 2 (1 kij )

The kij are called Binary Interaction Parameters (BIPS) and were introduced to account for deviations between theoretical and measured behaviour of binary mixtures: they have been empirically determined for the common EoS, i.e. PR and SRK.

BIPS between a component and itself are always zero. BIPS between the inorganics [N2, CO2 and H2S] and hydrocarbons are usually non-zero. Some authors have suggested that

Roxar Oxford

67

12/12/12

PVT Analysis

all hydrocarbon-hydrocarbon BIPS should be zero: some have suggested all bar Methanehydrocarbon should be zero.

5.4 Volume Translation

As we saw when we developed the 2-parameter EoS, the critical Z-factor of the popular EoS was fixed at constant valuables irrespective of the component or mixture of interest. A suggestion due to Peneloux et al. called Volume Translation or Volume Shift corrects this problem and vastly improves the prediction of volumetric properties such as liquid densities.

Essentially, the difference between theoretical and measured molar volume, on a component-by-component basis is computed by:

(5.43)	c	i	V EoS		V Obs
				i	i
where:
(5.44)	V Obs			M wi
		i		ref
				i

iref is the reference density of the ith component at standard conditions and ViEoS is the molar volume calculated from the EoS for the pure component, again at standard conditions. The coefficients ci are usually expressed in terms of dimensionless coefficients si by:

c

(5.44) si bi

i

where bi is the usual b-coefficient of the EoS. Since the ci are proportional to the bi, the linear mixing rule (5.39) is applicable to compute c-coefficients for mixtures:

	N
(5.45)	c xi ci

i 1

The corrected volume for a mixture is then calculated from:

(5.45)

Vmcorr Vm2P cm

where Vm2P is the molar volume predicted by the original 2-parameter EoS.

When this corrected volume appears in the Flash equations such as (A22), the same constant c appears on both sides of the equality. Therefore, the constant can be subtracted and we are left with the original Flash equations. This simplicity of treatment has meant a near universal take-up of the Peneloux model.

Roxar Oxford

68

12/12/12

PVT Analysis

6.Flash Calculations

Suppose we have an N-component fluid system that’s set of physical properties and moles are known. At some pressure and temperature, (p,T), we want to know if the mixture is 1-phase or 2-phase and if it’s 2-phase, what are the vapour and liquid mole fractions and the compositions of those phases: this is the role of the Flash calculation.

There are a number of ways to formulate this problem – all of which demand a high degree of mathematical rigour. The details of some of the mathematics used, in particular the role of Classical Thermodynamics, have been set out in the Appendix A.

Under isothermal conditions, we can find the state of a system by minimizing the Gibbs- Free-Energy (GFE), G. The GFE our N-component system is:

		N
(6.1)	G(0)	ni i

i 1

where ni and i are the moles and chemical potential of the ith component. Chemical potential plays a role similar to that of pressure and temperature. In the absence of gravity, a pressure difference will cause a fluid to flow. A temperature difference will allow heat to conduct. Chemical potential causes components to diffuse from regions of high to low chemical potential.

If our N-component system can form 2-phases, then the GFE will be given by:

		N	N
(6.2)	G(2)	niL iL niV iV
		i 1	i 1

where the moles in the liquid and vapour phases, (niL,niV) must satisfy the conservation of moles [mass] constraint:

(6.3)

ni niL niV

If G(2) is the minimum GFE, then we require:

(6.4)	G(2)	0
	n jV
Differentiating (6.2) gives:
(6.5)	G(2)	jV	N	iV
	n jV		jL niV	n jV
			i 1

	N	iL
	niL
		n jV
	i 1

Roxar Oxford

69

12/12/12

PVT Analysis

The second term uses (6.3), which implies n jL n jv 1 since the moles of feed, n, are

constant. The third and fourth terms are identically equal to zero because of the GibbsDuhem relationship16. With these conditions, (6.5) becomes:

(6.6)

iV iL 0

i 1, , N

From Appendix A.2 we see that chemical potential is related to fugacity coefficient, i, by:

(6.7)

ij i0 (T ) RT ln p RT ln xij RT ln ij

where i0 is the reference temperature, xij is the mole and (6.7) gives:

or ideal gas chemical potential, which is only a function of fraction of the ith component in the jth phase. Combining (6.6)

(6.8)

ln yi ln iV

ln xi ln iL 0

Using the Martin’s generalized 2-parameter EoS, we can now compute fugacity coefficients for the PR or SRK EoS and hence, in principle, solve (6.8) to calculate the phase split (xi,yi).

6.1 Successive Substitution (SS) Method

The K-value of the ith component was previously defined by (4.1) as:

y

(6.9) Ki xi

i

Substituting this into (6.8) gives:

(6.10)

ln Ki ln iL ln iV

Given a suitable set of initial K-values, we have the start of a process by which we can perform the Flash calculation. The most commonly used initial K-values come from the correlation due to Wilson:

p

ci

T

(6.11)

K

i

exp

5.3727 1

i

1

ci

p

T

which like the Hoffman et al. K-values, (4.2) and (4.3), give a good approximation for pressures less than 2000 psia and temperatures less than 200 oF.

16 Pressure, temperature and the chemical potentials are not independent. At constant pressure and temperature, the Gibbs-Duhem relationship is

N

ni d i 0 .

i 1

The proof requires several pages of mathematics and the interested reader is referred to Firoozabadi.

Roxar Oxford

70

12/12/12