Small-Group Discovery Experience (SGDE) The Phase Behavior Model (PBM) Execution Plan – Semester 2, 2016
The main objective of this special assignment is to introduce you to the application of computers in phase equilibrium calculations. We are using a modular approach for this project in which each module forms a stand-alone sub-project. The idea of this is to help you collaborating in a team work and achieving partial goals as the team go along. For this, each group will submit a 20-page (max) report.
The PBM Final Comprehensive Report
The final comprehensive report will contain, as a minimum:
1. Definition of the Problem and Objectives.
2. Description of the Technical Approach.
3. Description of the Algorithm and Flow Charts (a flow chart must be included for
each subroutine)
4. Testing of the Code.
5. Detailed discussion of the model and analysis of the results.
6. Suggestions/Limitations that have been found.
7. Summary.
8. Teamwork evaluation.
Max 1-page explaining each team member role and the tasks performed by individuals
9. Appendix: the Matlab Code (excluded from max 20 pages)
Note: It is optional to implement Peng-Robinson EOS OR Soave RK EOS; however you’ll get 10+ bonus point by implementing both EOSs.
An efficient code development is expected in this assignment. Code listings must be well documented and attached to the final report. All codes must be well documented; with
generous and detailed comments every so many lines. As a rule of thumb, for each 100 lines of code, you should have 30 lines of comments. The developed model should readily admit modification and addition. This is essential for future improvements.
If you are using function approach in your coding, you need to have a main code which calls other functions and print out the results.
Submission:
Submission:
o Submit your report as one file (Microsoft word) through myUni (one file
per group)
o Email your “Matlab code” to abbas.zeinijahromi@adelaide.edu.au (one file per group)
Submission Date: October 28, 2016
Group rules:
Self-enroll group of maximum 5 students
Deadline to enroll in groups September 9, 2016
o Students that have not enrolled in any group, will be assigned to a random group after the deadline
No change in group enrolment is accepted.
o Agree with your group before enrolling
It is a team work and the tasks can be divided between the team members; however, all members must participate in the whole project and be fluent on the project execution.
The Step-by-Step Guide for Building a Phase Behavior Predictor
Most analyses in reservoir and production engineering applications require predicting the behavior of complex mixtures of hydrocarbons in the liquid and gas state. For such applications, a reliable Phase Behavior Model that could perform such a prediction as a function of temperature and pressure becomes of paramount importance. This is especially true when the compositional behavior of petroleum fluids is to be modeled at its fullest. The objective of this project is to build the tool for doing this type of predictions. This document provides a step-by-step procedure for a modularized approach. If the
2
modules are done in sequential order, the end-result should be a phase behavior prediction package.
MODULE 1
Compressibility Factor Prediction THE OUTCOME: Z-root Finder Subroutine
The volumetric behavior of multi-component systems is described by Equations of State (EOS). EOS are analytical expressions that provide functional relationships between pressure, temperature and volume of a system—i.e., the PVT relationship or volumetric behavior of the system.
The first step in the building of a phase behavior model is to choose the EOS that best suits the purpose of the simulation. Peng-Robinson EOS (Peng and Robinson, 1976) and the Soave RK EOS (Soave, 1972) are, by far, the most popular cubic EOS currently in use—although the extension of the PBM to other equations of state follows similar logical steps. For the case of the PR EOS, Peng and Robinson proposed the following equation of state for multi-component hydrocarbon systems in 1976:
p RT (a)m (1a) ~~2 ~2
vbm v 2bmvbm
where:
p = pressure (psia)
T = temperature (R)
v~ molar volume (ft3/lbmol)
[am and bm, the mixture parameters, are calculated as a function of the properties of the pure components with the following mixing rules:
a nc nc cc aa (1 ) (1b) mij i j ij
ij
nc
bm cb (1c)
3
ii i
where:
oci 0.5 a R 2T 2 1 m 1 T 2
for all i; i = 1,2,…, nc
for all i; i = 1,2,…, n
i aipci i ri b o RTci
i bip
c
m i
if 0.49 iii
0.3746401.54226 0.269922
0.379642 1.48503 0.164423 2 0.016666 3 if 0.49
iiii i = Pitzer’s acentric factor of the i-th component in the mixture.
The most practical form of Equation 1a is the one obtained when expressed in its cubic version. This can be done in terms of either the specific volume ( v~ ) or the compressibility
~
factor ( Z Pv / RT ). The latter is the more useful, hence we recast Equation 1 as follows:
Z3aZ2bZc0 (2) 111
where:
1
c (ABB2 B3 ) 1
A ( a )m P R2T2
B bm P RT
a1 (1B)
b A3B2 2B
4
ci
(a)i = dimensional attraction parameter for the i-th component in the mixture,
bi = dimensional co-volume parameters for the i-th component in the mixture,
aio = attraction parameter constant for the i-th component in the mixture. The original Peng- Robinson EOS actually uses ai = 0.45724 for all substances. That is the value that should be used if no other information for ai is available.
bio = co-volume parameter constant for the i-th component in the mixture. The original Peng- Robinson EOS actually uses bi = 0.07780 for all substances. That is the value that should be used if no other information for bi is available.
R = universal gas constant (10.73 psia ft3 / lbmol R).
Tci = absolute critical temperature of the i-th component (R).
pci = critical pressure of the i-th component (psia).
Tri T = dimensionless reduced temperature. ci
T = absolute temperature (R).
nc = number of components in the hydrocarbon mixture.
ij = binary interaction coefficient between the i-th and j-th components, symmetric in i and j with ii = 0.
Equation 2 can be solved either numerically (using an iterative technique such as Newton-Raphson procedure) or analytically (analytic cubic solver), ), as explained later in this the handout.
GENERALIZED FORM OF CUBIC EOS’s
Several authors (see, for example, Martin (1979), Coats (1985), Danesh (2001)) have shown that all cubic EOS’s can be represented by single general formulations. One of the most useful generalized formulations is the one presented by Coats3, outlined below.
Z-factor form of the generalized EOS:
where:
Z3 (m m 1)B1Z2 12
AmmB2 (mm)B(B1)ZABmmB2(B1)0 121212
Anc ncccA, i j ij
ij
A (1 )(AA)0.5, ij ijij
2p Ao1m(1T0.5) ri, (SRKandPRonly)
i ai i ri
T2 ri
pri = p / pci,
Tri = T / Tci,
ij = binary interaction coefficient between the i-th and j-th components.
nc BciBi ,
i1 Bo pri,
i bi Tri
5
Whenever needed, the different equations of state are obtained by substituting the proper definitions of “m1”, “m2”, ai , and bi for the particular EOS. These definitions are given below for some of the most traditional EOS.
m 0.481.574 0.1762 iii
0.379642 1.48503 0.164423 2 0.016666 3 if 0.49
In the case of the RK EOS, use the following definition for Ai:
A o pri (RKonly) i ai T2.5
ri
SUBROUTINE PROGRAMMING – SUBROUTINE ZROOT
Develop a subroutine (ZROOT) to solve the roots of the Z-cubic Peng-Robinson EOS OR Soave RK EOS by implementing the analytical and/or numerical for Z-root polynomial calculations. Write a computer program with the mixture characterization as its input (i.e., composition and properties of each of its components) that calls the subroutine ZROOT to provide the compressibility factor(s) at a given pressure and temperature condition and the selected EOS. Use double precision for all variables carrying real numbers.
INPUT
* Gas Characterization (ci, Pci, Tci, i, MWi, Vci)
* Pressure TE*STTemIpTer!ature
Let us consider the following reservoir fluid characterization:
SRKEOS: PR EOS:
m i
0.3746401.54226 0.269922 if 0.49
iii
iiii
6
Equation of State
m1
m2
oa i
bo i
RK EOS Redlich Kwong
0
1
0.4274802
0.08664035
SRK EOS Soave RK
0
1
0.4274802
0.08664035
PR EOS
Peng Robinson
1+ 2
1- 2
0.457235529
0.077796074
PROGRAM MAIN
Call ZROOT
OUTPUT
Compressibility Factors Z (Zmax and Zmin)
Use the previous data as your input and validate your program by matching the following values, obtained using the PR EOS and SRK EOS:
7
Component
Composition (Molar), ci
Pci (psia)
Tci (oR)
i
MWi (lbm/lbmol)
Vci (ft3/lbm)
N2
0.0200
493.10
227.49
0.040
28.013
0.0510
C1
0.8866
666.40
343.33
0.008
16.043
0.0988
C2
0.0492
706.50
549.92
0.098
30.070
0.0783
C3
0.0246
616.00
666.06
0.152
44.097
0.0727
n-C4
0.0098
550.60
765.62
0.193
58.123
0.0703
n-C5
0.0098
488.60
845.80
0.251
72.151
0.0675
ij’s
N2
C1
C2
C3
n-C4
n-C5
N2
0.0000
0.0180
0.0390
0.0460
0.0470
0.0480
C1
0.0180
0.0000
0.0050
0.0100
0.0145
0.0182
C2
0.0390
0.0050
0.0000
0.0017
0.0032
0.0048
C3
0.0460
0.0100
0.0017
0.0000
0.0012
0.0024
n-C4
0.0470
0.0145
0.0032
0.0012
0.0000
0.0008
n-C5
0.0480
0.0182
0.0048
0.0024
0.0008
0.0000
P (psia)
T (F)
Peng-Robinson EOS
Soave-Redlich-Kwong EOS
A
B
Zmax
Zmin
A
B
Zmax
Zmin
1000
-55
0.5643
0.1079
0.4147
–
0.5232
0.1201
0.4428
–
900
-60
0.5234
0.0983
0.4296
–
0.4860
0.1095
0.4569
–
750
-70
0.4638
0.0840
0.4792
–
0.4319
0.0936
0.5074
–
1000
-70
0.6184
0.1120
0.3127
–
0.5758
0.1247
0.3404
–
500
-100
0.3750
0.0607
0.5735
0.1587
0.3522
0.0676
0.5985
0.1776
You should realize that the values in the columns vci (critical volumes, ft3/lbm) and MWi (molecular weights, lbm/lbmol) would not be needed for the specific calculations we are performing here. Nevertheless, your program must be able to read them and keep them stored for further developments.
Note that, using the composition molar fraction (ci) as an input for calculating the coefficients of the cubic equation of state gives roots of EOS with no physical meaning. If the vapor/liquid molar fraction (fng) was known before, we would be able to calculate vapor and liquid mole fraction (xi and yi) and consequently compressibility factors of liquid or vapor (as we will see in the next sections).
References
Coats, K.H.: Simulation of Gas Condensate Reservoir Performance, SPE Paper 10512, Journal of Petroleum Technology, October 1985.
Danesh, A: PVT and Phase Behaviour of Petroleum Reservoir Fluids, Elsevier, 2001.
Martin, J.J.: Cubic Equations of State—Which?, Industrial and Engineering Chemical Fundamentals, v. 81,
May 1979.
Peng, D. and Robinson, D. B.: A New Two-Constant Equation of State, Industrial and Engineering Chemical Fundamentals, v. 15, n. 1, 59-64, 1976.
Redlich, O., and Kwong, J.N.S.: On the Thermodynamics of Solutions: V – An Equation of State. Fugacities of Gaseous Solutions, Chem. Review (1949), v. 44, p. 233.
Soave, G.: Chem. Eng. Sci. (1972), v. 27, p. 1197.
8
9
HOW DO WE SOLVE A CUBIC EQUATION?
We have up to three methods at our disposal when it comes to solving a cubic equation such as x3 + ax2 + bx + c = 0. These are:
1. Analytical scheme
Given the cubic polynomial with real coefficients: x3 + ax2 + bx + c = 0, the first step towards solving the equation is to calculate the parameters:
Q a2 3b and R 2a3 9ab27c 9 54
Let M = R2 – Q3 be the discriminant. We then consider the following cases:
a)
If M < 0 (R2 < Q3), the polynomial has three real roots.
For this case, compute arccos( R ) and calculate the three distinct real roots as:
Q3
x2(2 Qcos2)a 33
x3(2 Qcos2)a 33
Note that x1, x2, x3 are not given in any special order; and that has to be calculated in radians.
If M > 0 (R2 > Q3), the polynomial has only one real root. Compute:
S3 R M T3 R M
and calculate the real root as follows: x S T a3
x1(2 Qcos)a 33
b)
Two complex roots (complex conjugates) may be found as well. However, they are of no interest for our purposes, and thus no formulas are provided. Such formulas can be found in the following suggested readings:
– W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Matlab,
2nd Edition, Cambridge Univ. Press, (1992), p. 179
– Spiegel, M., Liu, J., Mathematical handbook of formulas and tables, 2nd Edition, Schaum’s
Outline Series, McGraw Hill, p. 10
Sometimes, the equations for S and T listed above cause problems while programming. This usually happens whenever the computer/calculator performs the cubic root of a negative quantity. If you want to avoid such situation, you may compute S’ and T’ instead:
10
4.
S’sign(R)3abs(R) M
T’ Q / S’ (making T’=0 when S’=0) where: abs(R)=Absolute value of R
sign(R) = (+1) or (-1) if R is positive or negative respectively. It may be defined as: sign(R)=R/Abs(R)
and then the real root is: x1 S’T’ a3
Keep in mind the following useful relationships among the roots of any cubic expression: x1 + x2 + x3 = -a
x1x2 + x2x3 + x3x1 = +b x1 x2 x3 = -c
Self-exercise: Test your understanding on cubic roots calculations by analytical means by solving the following examples.
Solve the following cubic expressions:
x3-6×2 +11x-6=0 x1=1 x2 = 2 x3 = 3
x3+7×2 +49x+343=0 x1 = -7 x2 = 0 + 7 i
x3 = 0 – 7 i
x3+2×2 +3x+4=0 x1=-1.65063
x2 = -0.174685 + 1.54687 i
x3 = -0.174685 – 1.54687 i Solve the following cubic EOS in terms of volume (v),
v3 – 7.8693 v2 + 13.3771 v – 6.5354 = 0 One Phase (just one real root) v1 = 5.7357
v3 – 15.6368 v2 + 30.315v –14.8104 = 0 Three real roots,
v1 = 0.807582 (smallest, usually the liquid phase)
v2 = 1.36174 (rejected)
v3 = 13.4675 (largest, usually the gas phase) Solve the following cubic EOS in terms of volume (Z) (assume a pure component calculation):
Z3 –1.0595Z2 +0.2215Z–0.01317=0 Onephase, z=0.8045
Z3 – Z2 + 0.089 Z – 0.0013 = 0 Two Phases, zliquid = 0.0183012
zx = 0.0786609 (useless)
zvapor = 0.903038
Numerical Scheme
The Newton-Raphson method provides a useful scheme to solve for non-explicit variables from any form of implicit non-linear equation (not only cubic ones). Newton-Raphson is an iterative procedure with a fast convergence, although not always capable of providing an answer—because a first guess close enough to the actual answer must be provided.
2.
11
If we wanted to solve for “x” in any equation of the type f(x)=0, the method would provide a new estimate (new guess) closer to the actual answer, based on the previous estimate (old guess). This is written as follows:
xnewxold f(xold) f'(xold )
Considering a cubic equation f(x) = x3 + ax2 + bx + c = 0, the previous equation takes the form:
x n e w x o l d x o3 l d a x o2 l d b x o l d c 3xo2ld 2axold b
The iterations are kept until no significant improvement for “xnew” is achieved, i.e., | xnew – xold | < tolerance. An educated guess must be provided as the starting value for the iterations. If you are solving a cubic equation in Z (compressibility factor), some authors recommend taking Z= b P / RT as the starting guess for the compressibility of the liquid phase and Z=1 for the vapor root.
3. Semi-analytical Scheme
If you used the previous numerical approach and you have already calculated one of the roots of the cubic expression, the semi-analytical scheme can give you the other two real roots (in the case that they exist). By using the relationships given before, with the value ‘x1’ as the root already known, the two other roots are calculated by solving the system of equations:
x2 +x3 =-a–x1 x2x3 = -c/x1
which leads to a quadratic expression and where x1, x2, and x3 are the three roots of the cubic polynomial and a, b, c the coefficients of the cubic polynomial. This procedure can be reduced to the following steps:
1. Let x3 + ax2 + bx + c = 0 be the original cubic polynomial and “E” the root which is already known (x1=E). Then, we may factorize such cubic expression as:
(x - E) (x2+Fx+G) = 0
where: F = a + E G = -c / E
2. Solve for x2, x3 by using the quadratic expression formulae, x1 = E
MODULE 2
The Objective Function (Vapor-liquid equilibrium (VLE) Flash)
THE OUTCOME: The VLE Subroutine
The first problem we would have to address for our Phase Behavior Model is whether a given mixture at certain pressure and temperature conditions will exist as a single phase or will split into two phases (gas and liquid). If it does split, we would want to know the relative quantities of liquid and gas (fng and fnl) that coexist in equilibrium at such conditions. The equilibrium constants (Ki) play a key role during two-phase splitting prediction problems. The equilibrium constant, Ki, is defined as the ratio of the mole fraction of the i-th component in the gas phase (yi) to the mole fraction of the same component in the liquid phase (xi), i.e.,
K yi (3)
xi
From a molar material balance applied to a two-phase system in equilibrium and the definition of Ki [ef. Equation 3], we can derive the following expression (Rachford and Rice, 1952):
12
i
c (K 1) i11fng(Ki 1)
nc
g(f)ii 0 (4)
ng
Equation (4) is widely used and known as the Rachford-Rice Objective Function. Note that if all the Ki’s values were known, the only unknown left to solve for (in equation 4) is fng (the molar fractional amount of vapor that is found in the system).
SOLVING THE OBJECTIVE FUNCTION FOR fng
Since Equation (4) is a non-linear equation in one variable, Newton Raphson (NR) procedure is usually implemented. In general, Newton Raphson is an iterative procedure with a fast convergence rate. The method calculates a new estimate (fngnew) closer to the real answer than the previous guess (fngold), as follows:
g( f old )
f new f old ng (5)
ng ng g'( f old ) ng
In this iterative scheme, convergence is achieved when
fnewfold (6)
where is a small number (=1.0 x 10-10 is usually adequate). After solving for
liquid molar fraction and composition of each of the phases can be calculated as follows:
Vapor Phase Composition:
Where ci mole fraction of component i in the total mixture
x ci (7d)
13
ng ng
Liquid Molar Fraction: Percentage of Liquid: Percentage of Vapor:
Liquid Phase Composition:
f nl 1 f ng (7a) %L 100* fnl (7b) %V 100* fng (7c)
fng , the
i
1 fng (Ki 1)
yi Ki xi Ki ci (7e)
1 fng(Ki 1) Beyond Newton Raphson – A Bisection Strategy
Newton Raphson is an iterative procedure with a fast convergence rate but have some convergence problem (while g (fng) is so close to the zero).
Bisection method has a slower convergence rate comparing to the Newton Raphson method, but it’s a safe method and always converges to the solution. The procedures and criteria are as follows;
First, we should check if it’s necessary to look for exact fng or not. It means we don’t need to enter iteration procedure in cases of single phase liquid or single phase vapor. Since we know that Rachford and Rice equations are continuous and have a decreasing behavior, we can predict if our fng values will have a physical meaning or not (i.e. 0 <
𝑓 < 1) such that; 𝑛𝑔
If at fng = 0 g (fng) < 0 we are in single phase liquid If at fng =1 g (fng) >0 we are in single phase vapor
To find out that Newton Raphson method is converging to the answer or not, we should check the following conditions for fng in each iteration step;
1<𝑓<1
(1 − 𝐾𝑖 𝑚𝑎𝑥) 𝑛𝑔 (1 − 𝐾𝑖 𝑚𝑖𝑛)
where,
1 / (1-Ki min) is first positive asymptote of Rachford and Rice curve 1 / (1-Ki max) is first negative asymptote of Rachford and Rice curve
If our next guess of fng is not in the expected range we will have convergence problem and therefore the iteration procedure should switch to Bisection method as follows;
(𝑓 +𝑓 )
𝑓𝑛𝑒𝑤 = 𝑛𝑔
if g (fngnew) >0 then we update fng L value as; fng L = fngnew if g (fngnew) <0 then we update fng U value as; fng U = fngnew
and we can use Upper and Lower initial values of fng U =1 and fng L=0 .
Finally, we go back to first step and calculate fngnew until we reach the adequate minimum error ε as follow;
(𝑓 +𝑓 )
<𝜀
𝑛𝑔 𝑈
𝑛𝑔 𝐿
2
14
𝑛𝑔 𝑈
𝑛𝑔 𝐿
where,
ε =1.0 x 10-10
2𝑛
DO WE REALLY KNOW Ki’s?
The Ki value of each component in a real hydrocarbon mixture is a function of pressure, temperature and also of the composition of each of the phases. Since the compositions of the phases are not known beforehand, equilibrium constants are not as well. If they were known, the splitting calculation would be performed in a straightforward manner. This is because once the equilibrium constants of each component of the mixture are known for the given pressure and temperature of the system, both gas and liquid molar fractions, fng, and fnl can be calculated by solving the Rachford-Rice Objective Functions.
Numerous correlations have been developed throughout the years to estimate the values of Ki’s for each hydrocarbon component. In the following sections, we will develop a rigorous thermodynamic model for predicting equilibrium constants. In the meantime, we will use the popular Wilson’s empirical correlation as our first approximation. This correlation gives the value of Ki as a function of reduced conditions (Pri, Tri: reduced pressure and temperature respectively) and Pitzer’s acentric factor and is written as:
K 1 EXP5.37(1 )(1 1 ) (8) iPiT
15
ri ri SUBROUTINE PROGRAMMING – SUBROUTINE VLE
Develop a MATLAB subroutine (VLE) to solve for the gas molar fraction, fng, from the Rachford-Rice Objetive Function once all the values of Ki’s are given.
Use Ki‘s calculated with Wilson’s correlation (equation 8.)
INTEGRATE THEM!
Let us combine what we have done so far. In the program MAIN we had, and before calling ZROOT, let us call the VLE subroutine to find out how much gas and liquid we have in
INPUT
Equilibrium Constants (Ki’s)
SUBROUTINE VLE
OUTPUT
* Vapor Molar Fraction (fng) * Composition of the Gas
and Liquid Phase
the mixture. After doing this, we will have not only know how much liquid and gas we have, but also the composition of those phases.
Now we have three distinct compositions: the overall composition (ci), the gas composition (yi) and the liquid composition (xi). What we will do now is, instead of sending ‘ci’ to the ZROOT subroutine as we did before, we will be sending the corresponding ‘yi’ or ‘xi’ (one at the time) to obtain the compressibility factor of the gas and liquid (Zg and Zl) respectively.
What if you get more than one root for either the gas or liquid phase? The general rule of thumb here is to take the largest root if you are dealing with a gas phase and the smallest for liquids. Nevertheless, if you want to be completely rigorous, you are to choose the root that minimizes the free Gibbs Energy in either case. In cases where it may not be straight forward to identify a fluid as vapor or liquid, and when three roots are found for such systems, the intermediate root is ignored and the one which gives the lower Gibbs energy from the other two is selected. The selection of the smaller root identifies the fluid as the liquid-like, whereas that of the larger root indicates a vapor-like fluid. As described by Danesh (p. 176), the system molar Gibbs energy difference at two roots Zmax and Zmin is determined as:
A ln(Zmin m1B)(Zmax m2B) B(m m) (Z mB)(Z mB)
dG (Z Z )ln(ZminB) RT max min Z B
16
2 1 min 2 max 1
If the above is positive, Zmin is selected, otherwise, Zmax is selected. A, B, m1 and m2 are
max
defined in the Z-root calculation procedure.
TEST IT !
Using the same reservoir fluid, validate your program by matching the following values:
P (psia)
T (F)
fng
Phase
Peng-Robinson
Soave Redlich-Kwong
A
B
Z
A
B
Z
1000
-55
0.7774
Vapor
0.4646
0.0996
0.6086
0.4264
0.1109
0.6428
Liquid
0.9924
0.1367
0.2193
0.9430
0.1522
0.2468
900
-60
0.7977
Vapor
0.4319
0.0908
0.6261
0.3970
0.1011
0.6590
Liquid
0.9748
0.1277
0.1975
0.9292
0.1422
0.2226
750
-70
0.8177
Vapor
0.3833
0.0777
0.6566
0.3535
0.0865
0.6870
Liquid
0.9243
0.1124
0.1659
0.8851
0.1252
0.1872
1000
-70
0.6288
Vapor
0.5018
0.1029
0.5481
0.4622
0.1146
0.5815
Liquid
0.8454
0.1274
0.2258
0.7982
0.1419
0.2533
500
-100
0.8109
Vapor
0.3079
0.0559
0.7083
0.2865
0.0623
0.7326
Liquid
0.7419
0.0809
0.1128
0.7149
0.0901
0.1274
Furthermore, for following T and P condition if you use Newton Raphson method, it cannot converge to the correct solutions, but using Bisection method iteration, the correct results will be as follows;
17
P (Psia)
T (F)
fng
Phase
Peng-Robinson
Using Bisection method
Soave-Redlich-Kwong Using Bisection method
A
B
Z
A
B
Z
50
-100
0.971396
Vapor
0.034192
0.005824
0.971244
0.031973
0.006486
0.974097
Liquid
0.245992
0.01432
0.016512
0.243572
0.015947
0.018573
250
-100
0.91414
Vapor
0.159283
0.028318
0.860121
0.148435
0.031537
0.873438
Liquid
0.63892
0.051806
0.064369
0.624709
0.057696
0.072617
300
-100
0.89875
Vapor
0.18943
0.03387
0.83147
0.17645
0.03772
0.8472
Liquid
0.68266
0.058876
0.074549
0.66564
0.065569
0.084143
References
Rachford, H. H. and Rice, J. D.: Procedure for Use of Electrical Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium, Journal of Petroleum Technology, Technical Note 136, 19, October 1952.
MODULE 3
The IsoFugacity Criteria
THE OUTCOME: The Fugacity Subroutine
The Succesive Acceleration Method (SSM)
As we discussed in the previous section, reliable values for the equilibrium constants (Ki’s) must be obtained before we can solve the Rachford-Rice Objective Function (equation 4). Generally, first estimates for equilibrium constants are calculated by using Wilson’s empirical equation (equation 8).
However, Wilson’s correlation yields only approximate values for equilibrium ratios. Here what we do is to apply thermodynamic equilibrium considerations in order to obtain more reliable predictions for Ki’s.
For a system to be in equilibrium, any net transfer (of heat, momentum, or mass) must be zero. For this, all potentials (Temperature, Pressure and Chemical Potential) must be the same in all the phases. Provided that the temperature and pressure of both phases are the same, a zero net transfer for all components in the mixture results when all chemical potentials are the same. A restatement of this results in all fugacities for all components in each phase being equal. This is, the thermodynamic phase equilibrium equation:
fli fgi (9) where: fli = fugacity of the i-th component in the liquid phase
fgi = fugacity of the i-th component in the vapor phase
Since fugacity is a measure of the potential for transfer of a component between two phases, equal fugacities of a component in both phases results in a zero net transfer.
18
THE FUGACITY SUBROUTINE
At this point the importance of the fugacity to our phase behavior model has become clear. Fugacity calculations are necessary in order to implement the thermodynamic phase equilibrium condition represented by equation (9).
The fugacity coefficient is defined as the ratio of the fugacity of a material to its partial pressure. For a two-phase system which consists of N-components, the fugacity coefficient of component ‘i’ for both vapor and liquid phases can be written, respectively, as:
19
gi fgi yiP
fli
li
i = 1,2,...., nc (10)
lnln(ZB) i
nc
A 2 Ai j c j
B Z m B B
xiP
where fgi and fli are fugacity of component i in vapor and liquid phases, respectively. For the case of the generalized cubic EOS, the fugacity coefficient of the i-th component in a mixture is given by:
j1 (m1m2)B A
iln 2 i(Z1) BZm1B B
where:i = fugacity coefficient of the i-th component (dimensionless),
nc nc AcicjAij ,
ij
A (1 )(AA )0.5 , ij ijij
2p
Ao1m(1T0.5) ri , (SRKandPRonly)
i ai i ri T2 ri
... (11)
nc BciBi ,
i1
Bo pri. i bi Tri
Note that the Z-factor of the liquid and gas phases are needed in order to calculate of the fugacity coefficients in the liquid and gas phases, respectively. Please note that ci=xi for a liquid fugacity calculation; and ci=yi for a vapor fugacity calculation.
THE SUCCESIVE SUBSTITUTION METHOD (SSM)
Equation (9) shows that the fugacity of the components in each phase can be introduced to develop a criterion for thermodynamic equilibria. This is the fact that SSM (Successive Substitution Method) takes advantage of. It can be shown that Ki-values are related to fugacity coefficients ‘gi’ and ‘li’ as follows:
f(xP)yf
Ki li li i i li (12)
20
gi fgi(yiP) xifgi
Using the above equation, the correction step formulation to find K-values in SSM is
written as:
ynf n Kin1 i li
xi fgi f n
(13)
What SSM does is to basically update all previous equilibrium ratios (Ki) using the fugacities predicted by the equation of state (equation 11). This iteration method requires an initial estimation of Ki-values—for which the Wilson’s correlation is used. It can be easily concluded that the convergence criteria will be satisfied whenever the fugacity ratios of all the components in the system are close to unity. Such condition is achieved when the following inequality is satisfied:
nf 2
li 1 1014 (14)
Kin1 Kin li fgi
i fgi
SUBROUTINE PROGRAMMING –- SUBROUTINE FUGACITY SUBROUTINE SSM
a) Develop a Matlab subroutine to calculate fugacity coefficients for all the components of a mixture using the fugacity equation of the generalized cubic EOS. It is advisable to sketch the flow diagram of the algorithm for this routine before actually programming it.
b) Implement the SSM procedure in your Matlab code. Use Wilson’s equilibrium ratios for your initial estimate of Ki’s and iterate until the iso-fugacity criteria is satisfied.
TEST IT !
Using the same reservoir fluid, validate your program by matching the following values, utilizing the PR EOS and SRK EOS:
P T Soave Redlich-Kwong Peng-Robinson
(psia) (F) fng Zl Zv Ki fng Zl Zv Ki
1000 -55
900 -60
750 -70
1000 -70
500 -100
0.8774 0.2496
0.8741 0.2255
0.8714 0.1899
0.7593 0.2565
0.8454 0.1292
0.5646 C2
0.6005 C2
0.6506 C2
0.4613 C2
0.7184 C2
0.8778
0.8732
0.8692
0.7546
0.8404
0.2220
0.2001
0.1680
0.2302
0.1139
0.5285 C2
0.5656 C2
0.6187 C2
0.4278 C2
0.6936 C2
N2 3.0972 C1 1.3998
0.3667 C3 0.1411
nC4 0.0560
nC5 0.0223
N2 3.6788 C1 1.4805 0.3158 C3 0.1045
nC4 0.0359
nC5 0.0123
N2 4.7138 C1 1.6034 0.2525 C3 0.0669
nC4 0.0185
nC5 0.0051
N2 2.0841 C1 1.1806 0.4465 C3 0.2224
nC4 0.1134
nC5 0.0577
N2 6.7553 C1 1.6976 0.1502 C3 0.0260
nC4 0.0048
nC5 0.0009
N2 2.9885 C1 1.3807
0.3779 C3 0.1506
nC4 0.0621
nC5 0.0257
N2 3.5381 C1 1.4571 0.3256 C3 0.1118
nC4 0.0400
nC5 0.0143
N2 4.5084 C1 1.5722 0.2604 C3 0.0718
nC4 0.0208
nC5 0.0060
N2 2.0242 C1 1.1702 0.4606 C3 0.2368
nC4 0.1249
nC5 0.0660
N2 6.3868 C1 1.6547 0.1554 C3 0.0282
nC4 0.0055
nC5 0.0011
21
MODULE 4
Property Predictions
THE OUTCOME: The Molecular Weight Subroutine The Density Subroutine
The Liquid Viscosity Subroutine The Gas Viscosity Subroutine
Once SSM has been implemented, your program will be able to perform more reliable flash calculations. With the new values of the vapor and liquid compositions, we are ready to predict some important properties of both the liquid and vapor phases. Some of the most relevant are:
a) Molecular weight of the liquid and gas phases.
The molecular weight (MWa) of each of the phases is calculated as a function of the molecular weight of the individual components (MWi), once the composition of the gas (yi) and the liquid (xi) are known:
(15a) (15b)
22
n
i1
MWg
MWl xMW
n
i1
ii
yMW ii
b) Density of the liquid and gas phases
The density of the phase ‘a’ is calculated using its compressibility factor (Za) as predicted by the Peng Robinson equation of state. From the real gas law, the density is expressed as:
MWl ; MWg (16,17) l v
LV
Expression (16 and 17) are applicable for both the gas and liquid density. In either case, the proper value for MWg or MWl and Zg or Zl has to be used.
c) Viscosity of the gas, using Lee-Gonzalez-Eakin method.
Lee, Gonzalez and Eakin (1966) presented the following correlation for the calculation of the viscosity of a natural gas:
yv
g 1104kvEXPxv g (18)
9.40.02MW T1.5 where: kv g
62.4
20919MWg T yv 2.40.2xv
xv 3.5 986 0.01MWg T
In this expression, temperature is given in (oR), the density of the fluid (g) in lbm/ft3 (calculated at the pressure and temperature of the system), and the resulting viscosity is expressed in centipoises (cp).
d) Viscosity of the liquid, using Lohrenz, Bray & Clark correlation.
Lohrenz, Bray and Clark (1964) proposed an empirical correlation for the prediction of the viscosity of a liquid hydrocarbon mixture from its composition. Such expression, originally proposed by Jossi, Stiel and Thodos (1962) for the prediction of the viscosity of dense-gas mixtures, is given below:
(19a)
23
4
* 10.10230.023364 0.0585332 0.0407583 0.00937244 1104
l m r r r r where:
l=liquid viscosity (cp)
*=viscosity at atmospheric pressure (cp) m=mixture viscosity parameter (cp-1) r=reduced liquid density (unitless)
The coefficients in expression (19a) are presented as they were proposed by Jossi, Stiel and Thodos (1962), without the typo that is found in Lohrentz et al. (1964) original paper.
The following expressions give the definitions of the parameters needed to employ expression (19a). They are presented in the English system of units, although the original paper used scientific units.
As for the viscosity of the mixture at atmospheric pressure (*), Lohrentz et al. (1964) suggested to use the Herning & Zipperer equation:
x * MW iii
* i (19b) xi MWi
i
where:
xj = mole composition of the i-th component
MWi = molecular weight of the i-th component (lbm/lbmol)
iviscosity of the i-th component at low pressure (cp). For this parameter, Stiel & Thodos correlations are recommended:
34105T0.94 * ri
i i
forTri<1.5 ri for Tri > 1.5
(19c) (19d)
24
5 0.625 17.7810 4.58T 1.67
i*
where Tri is the reduced temperature of the i-th component (T/Tci) and i is the
i
viscosity parameter of the i-th component given by:
5.4402T 1 / 6 ci
(19e)
As for the mixture viscosity parameter m), Lohrentz et al. (1969) presented the following expression:
5.4402T1/6
m pc (19f)
where:
Tpc = pseudocritical temperature (oR)
Ppc = pseudocritical pressure (psia)
MWl=liquid mixture molecular weight (lbm/lbmol)
And finally, for the calculation of the reduced density of the liquid mixture (r), from Lohrentz et al. (1969) work it is derived:
i
MW P2 / 3 i ci
MW P2/3 l pc
l l V (19g)
r MWpc pc l
where pc is the mixture pseudocritical density (lbm/ft3) and Vpc pseudocritical volume of the mixture per unit mole (ft3/lbmol). All mixture pseudocritical properties are calculated using Kay’s mixing rule, as shown:
Tpc=xiTci
Ppc=xiPci
Vpc=xiVci
where Tci is given in oR, Pci in psia, and Vci in ft3/lbmol. Note that if the pure component critical volumes are given in a mass basis (ft3/lbm), each of them has to be multiplied by the corresponding molecular weight. Lorentz et al. (1969) also proposed a correlation for the estimation of critical volumes of C7+ fractions.
Selected References:
Lee, A.L., Gonzalez, M.H., Eakin, B.E.: The Viscosity of Natural Gas, SPE Paper 1340, Journal of Petroleum Technology, 997-1000, August 1966.
Lohrenz, J., Bray, B.G., Clark, C.R: Calculating Viscosities of Reservoir Fluids from their compositions, SPE Paper 915, Journal of Petroleum Technology, 1171-1176, October 1964.
SUBROUTINE PROGRAMMING –- SUBROUTINE MWEIGHT SUBROUTINE DENS
SUBROUTINE VISCG SUBROUTINE VISCL
Provide subroutines to your Matlab code for the prediction of the following properties, using the suggested approaches (Note: you can alternatively use other methods presented during the lectures)
a) Molecular weight of the liquid and gas phases,
b) Density of the liquid and gas phases, using the compressibility factors of
each phase
c) Viscosity of the gas, using Lee-Gonzalez-Eakin method,
d) Viscosity of the liquid, using Lohrenz, Bray & Clark correlation.
(19h) (19i) (19j)
is the
25
TEST IT !
Using the same reservoir fluid characterization, validate your program by matching the following values, which use the PR EOS and SRK EOS:
26
SRK EOS
P (psia)
T (F)
fng
MW (lb/lbmol)
Density (lbm/ft3)
Viscosity (cp)
1000
-55
0.8774
Vapor
17.4895
7.1281
0.0135
Liquid
26.7524
24.6672
0.0538
900
-60
0.8741
Vapor
17.3242
6.0500
0.0124
Liquid
27.6528
25.7171
0.0590
750
-70
0.8714
Vapor
17.1302
4.7189
0.0111
Liquid
28.7567
27.1383
0.0676
1000
-70
0.7593
Vapor
17.4776
9.0538
0.0156
Liquid
22.2438
20.7265
0.0398
500
-100
0.8454
Vapor
16.8038
3.0277
0.0093
Liquid
28.5814
28.6447
0.0816
PR EOS
P (psia)
T (F)
fng
MW (lb/lbmol)
Density (lbm/ft3)
Viscosity (cp)
1000
-55
0.8778
Vapor
17.5298
7.6333
0.0140
Liquid
26.4946
27.4691
0.0746
900
-60
0.8732
Vapor
17.3556
6.4344
0.0127
Liquid
27.3645
28.6831
0.0857
750
-70
0.8692
Vapor
17.1512
4.9683
0.0113
Liquid
28.4179
30.3130
0.1055
1000
-70
0.7546
Vapor
17.5124
9.7827
0.0166
Liquid
22.0466
22.8898
0.0490
500
-100
0.8404
Vapor
16.8127
3.1377
0.0094
Liquid
28.1661
32.0097
0.1421
MODULE 5
Facing Convergence Problems around CP
THE OUTCOME: The Accelerated Substitution Method (ASSM)
BONUS (+5): Instead of implementing the ASSM algorithm stated in part A, you can alternatively implement a more effective ones (part B). Find details at the end of the module. The ASSM algorithm in part A is simpler to implement, but it fails to perform in several instances.
When the system is close to the critical point and fugacities are strongly composition- dependent, a slowing-down of the convergence rate of the SSM (Succesive Substitution Method) is to be expected. In an attempt to avoid slow convergence problems, some methods have been proposed. Among the most popular ones are the Minimum Variable Newton Raphson (MVNR) Method and the Accelerated and Stabilized Successive Substitution Method (ASSM).
Part A : Rinses et al. (1981) ASSM algorithm
The ASSM is basically the accelerated version of the SSM procedure and thus follows a similar theory. Such procedure is implemented to accelerate the calculation of finding Ki- values, especially in the region close to critical point where the use of the SSM alone will not be efficient. The ASSM technique was presented by Rinses et al. (1981) and consists of the following steps:
1. Use the SSM technique to initiate the updating of the Ki-values the first time.
2. Check all following criteria at every step during iterations using the SSM:
27
nc
(Rrnew 1)2 i
0.8 (Rrold 1)2
1
| fngnew – fngold | < 0.1
1
nc
i
nc
1
0 < fngnew < 1 (20)
These criteria assure that you have sufficiently close proximity to the conditions that ensure efficiency of the method. Rri is the ratio of liquid fugacity to gas fugacity of the i-th component and ‘fng’ is molar fraction of the two-phase system.
Rr fli (21)
fgi
3. If the system satisfies ALL above criteria, the iteration technique is then switched from the SSM to the ASSM. Otherwise, SSM is used for the update of the Ki-values. To update Ki-values by implementing ASS Method, the following expressions are used:
Ki new= Kiold Rri i
where i =[ ( Rriold –1) / (Rriold - Rrinew ) ]
4. Once all the criteria in step (2) are satisfied, skip step (2) for the subsequent iterations and use the ASSM technique to update Ki-values until convergence is attained, unless it does not give acceptable new estimates (as stated next).
5. When ASSM is used, it must always be tested to show that it leads to an improved solution (i.e., that it brings fugacity ratios closer to unity). If not, it must be rejected and switched back to SSM.
Even though we are outlining Risnes et al. version of the accelerated Successive Substitution Method, there are several other published algorithms whose main purpose have also been accelerating the successive substitution method. Fussel and Yanosik (1978), Michelsen (1982) and Mehra et al. (1983) are examples of such attempts. Risnes et al. version is the easiest and most straightforward to implement but it is subjected to limitations. For the purpose of this project, you may select any of these algorithms (or any other that you may find more efficient) and verify its performance and limitations.
Since the Rinses et al. (1981) ASSM method might fail to perform in several instances, we represent different and more effective ASSM algorithms presented by Mehra et al. (1983) in the following section.
105 (Rrnew 1)2 103 i
28
i
Part B : Mehra et al. (1983) ASSM algorithms
Summaries of proposed algorithms by Mehra et al. (1983) are as following and for more details you can see the reference.
As we mentioned in Module 3, in the successive substitution method, the next step is to correct the Ki according to;
𝐾𝑛𝑒𝑤 =𝐾𝑜𝑙𝑑 𝑅𝑟𝜆𝑖 𝑖𝑖𝑖
where Rri fli gi
Also, the equilibrium state is defined as the state of minimum Gibbs free energy among the manifold of states of constant temperature and pressure (Callen, 1960). For a vapor- liquid system, the Gibbs free energy is given by:
𝑔 =𝜕(𝐺/𝑅𝑇)/𝜕𝑛 =𝑙𝑛𝑓 −𝑙𝑛𝑓 𝑖 𝑔𝑖 𝑔𝑖 𝑙𝑖
and the above equations can be re-written as;
∆𝑙𝑛𝐾𝑖 = − 𝜕(𝐺/𝑅𝑇)/ 𝜕𝑛𝑔𝑖
or, in vector notation,
̅̅̅̅̅̅̅
∆𝑙𝑛 𝐾 = − 𝑔̅ = − ∇𝐺/𝑅𝑇
Mehra et al. (1983) modified the successive substitution step of above equation by introducing a step length λ , i.e.
∆𝑙𝑛 𝐾 = − 𝜆 𝑔̅
Gradient methods are accelerated by choosing the step length λ, in some optimal manner. There are several alternatives depending on the choice of the minimization criterion. Here, three alternatives proposed for accelerating the successive substitution method for flash calculations all give significant improvement in the convergence rate. Of these three, the third is simplest.
For the ASSM algorithm was presented by Mehra et al. (1983.) these steps should be followed;
1. The same as before, we use the SSM technique to initiate the updating of the Ki- values the first time.
2. Check the following criteria to reduce the consequent erratic behavior of the algorithms;
29
where,
∆𝑙𝑛𝐾𝑖 = 𝑙𝑛𝐾𝑖𝑛𝑒𝑤 − 𝑙𝑛𝐾𝑖𝑜𝑙𝑑
where,
ḡn = gradient of Gibbs free energy = 𝑙𝑛(𝑓 /𝑓 )
𝑔𝑖 𝑙𝑖
λn = acceleration parameter and calculated as; For first alternative, we have;
𝑔̅ 𝑇 𝑈−1 𝑔̅
𝜆=[ 𝑛−1 𝑛−1 ]𝜆
𝑛 𝑔̅𝑛−1𝑇 𝑈−1 (𝑔̅𝑛−1 − 𝑔̅𝑛) 𝑛−1
|∆𝑙𝑛𝐾𝑖| = |𝜆 𝑔̅𝑖| ≤ 6
0<𝑓𝑛𝑒𝑤 <1 𝑛𝑔
𝑎𝑙𝑙 𝑖
3. For updating Ki, we can use three algorithms introduced by Mehra et al. as following which all three are presented in the form of;
∆𝑙𝑛 𝐾 = − 𝜆𝑛 𝑔̅𝑛
∆𝑙𝑛𝐾𝑖 = equilibrium ratio deference = 𝑙𝑛𝐾𝑖𝑛𝑒𝑤 − 𝑙𝑛𝐾𝑖𝑜𝑙𝑑
30
For second alternative, we have;
𝑔̅ 𝑇 (𝑔̅
𝜆𝑛 =[ 𝑛−1 𝑛−1
(𝑔̅ −𝑔̅ )𝑇 (𝑔̅
𝑛−1 𝑛 𝑛−1 𝑛
For third alternative, we have;
𝑔̅ 𝑇𝑔̅
𝜆=[ 𝑛−1 𝑛−1 ]𝜆
− 𝑔̅ ) 𝑛
]𝜆𝑛−1
where,
T = Transpose of matrix
U = transformation matrix and calculated as;
or simply;
where,
𝑢𝑖𝑗 =( 𝑛 )[−1+𝛿𝑖𝑗 𝑐𝑖] 𝐿 𝑉 𝑥𝑖 𝑦𝑖
𝑛 𝑔̅𝑇(𝑔̅−𝑔̅)𝑛−1
−1 −1 𝑢𝑖𝑗 = 𝜕𝑙𝑛𝐾𝑖/𝜕𝑛𝑖𝑔 = −[𝐿 + 𝑉
−𝑔̅ )
𝑛−1 𝑛−1
𝑛
] + 𝛿𝑖𝑗[1/𝑛𝑖𝑔 + 1/𝑛𝑖𝑙]
δ = Kroneker delta which is 1 if i = j and 0 in other cases.
L=liquidmole= 𝑛× 𝑓 𝑛𝑙
V=vapormole=𝑛× 𝑓 𝑛𝑔
note that we assume n = 1 mole in calculations.
In all three algorithms, only positive values are expected for λ and if a negative value arises the absolute value is taken. Also the algorithms can be initiated with λ = 1. In the PBM module, all alternative of Mehra et al. (1983) algorithm for ASSM loop can be included along with the option for the Rinses et al. (1981) algorithm.
4. Note that, after entering the ASSM loop we should stay in the loop until the end. There is an exception as it’s stated in step 5;
5. During ASSM loops it must always be tested to see if it leads to an improved solution (i.e., that it brings fugacity ratios closer to unity), i.e;
|(𝑅𝑟–1)𝑛𝑒𝑤|< |(𝑅𝑟–1)𝑜𝑙𝑑|
If not, it must be rejected and switched back to SSM.
6. The SSM/ASSM iteration will update all previous equilibrium ratios Ki using the fugacities predicted by the equation of state and we continue this iteration loop until fugacity ratios of all the components became close to one and we reach the adequate minimum error as follow;
nf 2
li 1 1014
Example A: Mehra et al. Algorithm
Here is the example of using Mehra et al. Algorithm and PR EOS for updating of Ki-values. As we said, we use Willson correlation for initial Ki guess and after that, we use the SSM technique to initiate the updating of the Ki-values the first iteration. For example, at P=700 psia and T=-75 F, using the SSM technique, Ki-values of Willson ( K(0) ) and the first SSM iteration ( K(1) ) and relative g1 values are as follows;
31
i fgi
K (0)=
6.920359 1.71031 0.080742
0.00962 0.001397 0.000225
K (1)=
4.9063 1.5887 0.229 0.056821 0.014869 0.0038811
g1 =ln (fgi/ fli) =
0.34394 0.073757 -1.0425 -1.7761 -2.3651 -2.8484
and as we said, the algorithms can be initiated with λ1 = 1 and also we define gT as follows; g1Tij= g1 ji = [ 0.34394 0.073757 -1.0425 -1.7761 -2.3651 -2.8484 ]
For updating Ki, we use third algorithm introduced by Mehra et al. as follows;
∆𝑙𝑛 𝐾 = − 𝜆𝑛 𝑔̅𝑛
𝐾𝑛𝑒𝑤 = exp(−𝜆 𝑔̅ +𝑙𝑛𝐾𝑜𝑙𝑑) 𝑖𝑛𝑛𝑖
0.023216 0.00019143 -0.040955 -0.073526 -0.10502 -0.13772
4.7885
1.5884 0.23904 0.061371 0.016599 0.0044835
and after 14 iteration we will find the Ki values as follows;
32
g2 =ln (fgi/ fli)=
K (2)=exp(-λ2 g2+ln K(1) ) =
Note that we calculate λ2 from;
𝑔̅ 𝑇 𝑔̅
𝜆2 =[ 1 1 ]𝜆1=1.0476
𝑔̅ 𝑇 (𝑔̅ − 𝑔̅ ) 112
K=
4.8253
1.5948 0.23894 0.061221 0.016528 0.0044575
TEST IT !
Test your program by running the following isothermal experiment:
T ( F)
P (psia)
-75
700
-75
800
-75
900
-75
920
-75
950
-75
970
-75
1000
-75
1010
-75
1020
For each set of conditions, you should be able to print number of iterations, gas and liquid phase compositions, liquid molar fraction, densities and viscosities.
The results for above Temperatures and Pressures, using PR/SRK EOS’s are as follow;
PT (psia) (F)
700 -75 800 -75 900 -75 920 -75 950 -75 970 -75
1000 -75 1010 -75 1020 -75
700 -75 800 -75 900 -75 920 -75 950 -75 970 -75
1000 -75 1010 -75 1020 -75
PT psia F
700 -75 800 -75 900 -75 920 -75 950 -75 970 -75
1000 -75 1010 -75
1020 -75
700 -75 800 -75 900 -75 920 -75 950 -75 970 -75
1000 -75 1010 -75
1020 -75
fng
0.864 Vapor 0.825 Vapor 0.759 Vapor
0.74 Vapor 0.705 Vapor 0.675 Vapor 0.613 Vapor 0.583 Vapor 0.538 Vapor
0.864 Liquid 0.825 Liquid 0.759 Liquid
0.74 Liquid 0.705 Liquid 0.675 Liquid 0.613 Liquid 0.583 Liquid 0.538 Liquid
fng
0.86711 Vapor 0.82882 Vapor 0.76491 Vapor 0.74626 Vapor 0.71198 Vapor 0.68324 Vapor 0.62535 Vapor
0.599 Vapor
0.56439 Vapor
0.86711 Liquid 0.82882 Liquid 0.76491 Liquid 0.74626 Liquid 0.71198 Liquid 0.68324 Liquid 0.62535 Liquid
0.599 Liquid 0.56439 Liquid
MW Density*
ASSM # of iterations
(lb/lbmol)
(lbm/ft3)
MW Density* (lb/lbmol) (lbm/ft3)
17.0622 4.3565 17.0935 5.545 17.1989 7.2569 17.2394 7.7176 17.3257 8.5413 17.4113 9.2208 17.6296 10.598 17.7537 11.237
Viscosity* (cp)
0.011 0.012 0.013 0.014 0.015 0.016 0.018 0.019
ASSM #of iterations,
17.942 12.08
0.02
0.057 0.045 0.043 0.039 0.037 0.033 0.031
0.029
con con verge verge
26.0397 25.195 23.2648 22.336 22.6998 21.651 21.8366 20.509 21.2425 19.645 20.2862 18.072 19.9262 17.404
19.5097 16.561
18 7 9 16 16 8 15 15
Peng-Robinson EOS Rinses et al.
Mehra et al. (1983)
(1981)
Viscosity* Total ASSM (cp) # # of of
ite. ite.
Total # of iterations
Algorithm # 123123
17.0802 4.5725 0.011 19 0 8 12 12 7 11 11 17.1148 5.876 0.012 25 0 9 15 16 8 14 15 17.2253 7.7747 0.014 22 0 11 13 13 10 12 12
17.267 8.2879 0.015 42 0 13 18 14 12 17 13
17.355 9.2062 0.016 59 0 13 20 19 12 19 18 17.4412 9.9641 0.017 77 0 17 18 21 16 17 20 17.6578 11.502 0.019 140 0 20 23 22 19 22 21 17.7801 12.22 0.021 47 10 19 35 17 18 34 16 17.9665 13.179 0.022 43 14 23 27 28 22 26 27
28.4656 30.678 0.112 19 0 8 12 12 7 11 11 25.7306 28.063 0.083 25 0 9 15 16 8 14 15 23.0332 24.772 0.059 22 0 11 13 13 10 12 12 22.4875 23.983 0.054 42 0 13 18 14 12 17 13 21.6548 22.668 0.048 59 0 13 20 19 12 19 18 21.0816 21.668 0.044 77 0 17 18 21 16 17 20 20.1551 19.835 0.038 140 0 20 23 22 19 22 21 19.8037 19.047 0.036 47 10 19 35 17 18 34 16 19.3932 18.039 0.033 43 14 23 27 28 22 26 27
Soave-Redlich-Kwong EOS Rinses et al.
Mehra et al. (1983)
(1981)
Total ASSM # # of of ite. ite.
Total # of iterations,
Algorithm # 123123
16 8 8 14 15 7 13 14
18 7 42 0 47 0 54 0 69 0
137 0 193 0
Not Not
10 16 18 9 15 17 13 19 23 12 18 22 11 17 20 10 16 19 13 19 20 12 18 19 16 21 18 15 20 17 18 20 22 17 19 21 19 25 23 18 24 22
24 49 31 23 48 30 28.821 27.461 0.07 16 8 8 14 12 7 13 11
42 0 47 0 54 0 69 0
137 0 193 0
Not Not
con con verge verge
11 19 13 10 18 12 13 17 14 12 16 13 13 19 19 12 18 18 17 21 21 16 20 20 20 20 22 19 19 21 19 25 17 18 34 16
23 49 28 22 26 27
Algorithm #
Algorithm #
33
References
Risnes, R., Dalen, V., and Jessen, J. I.: Phase Equilibrium Calculations in the Near-Critical Region, Elsevier Sequoia, SA, Proceedings of the European Simposium on EOR, 329-350, Lausanne, 1981.
Fussel, D.D. and Yanosik, J.L.: An Iterative Sequence for Phase Equilibria Calculations Incorporating the Redlich-Kwong Equation of State, SPE Journal, p. 173-182, June 1978.
Michelsen, M.L.: The Isothermal Flash Problem: Part II, Phase Split Calculation. Fluid Phase Equilibria, v. 9, p. 21-40, 1982.
Mehra, R.K., Heidemann, R.A., and Aziz, K: An Accelerated Successive Substitution Algorithm. Canadian J. of Chemical Engineering, v. 16, p. 590-596, August 1983.
Curtis H. Whitson et al, “Phase Behavior”, SPE Monograph series, Volume 20, 2000
34
MODULE 6
Handling 1- in a VLE Calculation THE OUTCOME: The Phase-Stability Method
Interesting enough, one of the most difficult aspects of making VLE calculations is not the two-phase splitting calculation itself; but knowing whether or not a mixture will actually split into two (or even more) phases for a pressure and temperature condition.
A single-phase detection routine has to be simultaneously introduced at this stage to detect whether the system is in a true-single-phase condition at the given pressure and temperature or it will actually split into two-phases. Several approaches may be used here: Bring-Back technique outlined by Risnes et al. (1981), Phase Stability Criteria introduced by Michelsen (1982), among others.
MICHELSEN’s STABILITY TEST
Michelsen (1982) suggests creating a second-phase inside any given mixture to then verify whether such a system is stable or not. It is the same idea behind the bring-back procedure (Risnes et al., 1981) but this test additionally provides straightforward interpretation for the cases where trivial solutions are found (Ki’s 1). The test must be performed in two parts, considering two possibilities: that second phase can be either vapor-like or liquid-like. For this test to be successful, the root that minimizes Gibbs free energy should be selected whenever the cubic polynomial yields three roots.
1. Calculate the mixture fugacity (fzi) using overall composition zi. If two roots are found, use the one that yields the minimum Gibbs free energy. Use Danesh method for your root selection (Danesh book, p. 176)
2. Create a vapor-like second phase,
a- Use Wilson’s correlation to obtain initial K-values. b- Calculate second-phase mole numbers, Yi:
Yi ziKi (22) c- Obtain the sum of the mole numbers,
35
n
SV Y (23)
i i
d- Normalize the second-phase mole numbers to get mole fractions: yi SYi (24)
v
e- Calculate the second-phase fugacity (fyi) using the corresponding EOS and the previous composition [equation (24)].
f- Calculate corrections for the K-values:
Rfzi 1 (25)
fyi Sv
K(n1) K(n)R (26) iii
g- Check if:
g.1) Convergence is achieved:
n R 12 11010 (27) i
i
g.2) A trivial solution is approached:
n lnK21104 (28)
36
i
i i
-If a trivial solution is approached, stop the procedure.
-If convergence has not been attained, use the new K-values and go back to step (b).
3. Create a liquid-like second phase,
Follow the same previous steps by replacing equations (22),(23),(24), and (25) by (29), (30), (31), and (32) respectively.
Yi zi /Ki (29) SL n Yi (30)
i
xi SYi (31) L
Ri fxi SL (32) fzi
INTERPRETATION:
Mixture is stable (i.e., the single-phase condition prevails) if: -BothtestsyieldsS<1(SL <1andSV <1),
- or both tests converge to trivial solution,
- or one test converges to trivial solution and the other gives S < 1.
Only one test indicating S > 1 is sufficient to determine that the mixture is unstable and that the two-phase condition prevails. Same conclusion is made if both test can give S > 1, or if one of the tests converge to the trivial solution and the other gives S > 1.
TEST IT !
Test your program to perform the calculation at the following temperature and pressure conditions:
37
Test number
T (0F)
P (psia)
1
-100
50
2
-100
250
3
-100
500
4
-100
1000
5
-100
1250
6
-100
1500
7
-100
1750
8
-100
2000
9
-100
2250
10
-100
2500
At this
a) b) c) d) e) f) g) h) i)
stage of the project, you will be able to generate and print out:
Detection of Single-Phase/Two-Phase Conditions, number of iterations (for SSM/ASSM convergence), liquid and vapor fractions,
gas and liquid phase composition,
equilibrium ratios,
molecular weight of each phase, fugacities of each component, compressibility factors of each phase, densities and viscosities of both phases.
Some of the results you should be able to match are:
38
Two-Phase (fnl=0.02842)
Two-Phase (fnl=0.15955)
0.9712
Peng-Robinson EOS
P (psia)
50
T (oF)
-100
Phase Detection
Zg
Zl
MWg
MWl
0.0165
17.56
54.86
250
-100
Two-Phase (fnl=0.07448)
0.8572
0.0655
17.00
38.81
500
-100
0.6939
0.1137
16.81
28.17
1000
-100
Single Phase (fnl=1.0)
—
0.2402
—
18.62
1250
1500
-100
—
0.2858
—
18.62
-100
Single Phase (fnl=1.0)
—
18.62
1750
-100
Single Phase (fnl=1.0)
Single Phase (fng=1.0)
—
—
0.3306
0.3743
—
18.62
Two-Phase (fnl=0.02899)
Two-Phase (fnl=0.15458)
0.7183
Soave-Redlich-Kwong EOS
P (psia)
T (oF)
Phase Detection
Zg
Zl
0.1291
MWg
MWl
50
-100
0.9741
0.0186
17.55
54.79
250
500
-100
Two-Phase (fnl=0.07386)
0.8706
0.0743
16.99
39.12
-100
16.80
28.58
1000
1250
-100
-100
Single Phase (fnl=1.0)
Single Phase (fnl=1.0)
—
—
0.2684
—
18.62
0.3194
—
18.62
1500
-100
Single Phase (fnl=1.0)
—
0.3693
—
18.62
1750
-100
Single Phase (fng=1.0)
—
0.4182
—
18.62
References
Risnes, R., Dalen, V., and Jessen, J. I.: Phase Equilibrium Calculations in the Near-Critical Region, Elsevier Sequoia, SA, Proceedings of the European Simposium on EOR, 329-350, Lausanne, 1981.
Michelsen, M.: The Isothermal Flash Problem. Part I. Stability. Fluid Phase Equilibria, v. 9, p. 1-19, 1982.
MODULE 7
Application: Phase Envelope Construction
THE OUTCOME: A: Brute-Force Phase-Envelope Predictor B: Michelsen’s Phase-Envelope Predictor
In order to have more complete information about the phase behavior of a hydrocarbon fluid, it is often useful to generate a phase envelope, which designates the boundary between single-phase fluid (vapor or liquid) and two-phase mixture at equilibrium conditions. Many possible algorithms can be used to construct a phase envelope, but for the purposes of implementing a straight-forward application of the “Phase Stability Subroutine”, we will test the two procedures described below.
PROCEDURE “A”: The “Brute-Force” Approach
The “brute-force” approach outlined here is provided for illustration purposes only. It should be borne in mind that this procedure is numerically inefficient and should be substituted for a more appropriate algorithm (e.g., SPE 11198) for any advanced applications.
1. Select a starting (low) pressure. This could be close to atmospheric, or 20-50 psia.
2. Start at a point found in the two-phase region at this pressure (i.e., select a temperature
at which the stability test predicts the formation of an unstable phase). You could use Wilson Ki’s to estimate the temperature at which fng would be equal to 0.5 at the assumed low pressure.
3. At the given pressure, lower the temperature progressively (dT not higher than 1 F) until a bubble point condition is reached (i.e., the stability test goes from a two-phase prediction to a single-phase prediction). This point could also be a dew-point condition in case that p > pc and T > Tc of the mixture.
4. Starting at the initial point, raise the temperature until a dew point condition is reached (i.e., stability starts predicting a single-phase). This point could also be a bubble-point condition in case that p > pc and T < Tc of the mixture.
5. Increment pressure.
6. Use previous results to get an initial temperature estimate for next set of points. You
could take the average of the previous two saturation temperatures, for example.
7. Close to the critical point, you can encounter either two bubble or two dew points instead
of one each as suggested above.
8. Continue until both predictions (step 3 and step 4) yield approximately the same result
(i.e., you have reached the cricondenbar and the envelope is now complete).
39
For the case of the brute-force approach, please note that you should only call the Phase Stability subroutine to generate the phase envelope indicated here.
PROCEDURE “B”: Pedersen’s approach (+10 bonus points)
Please refer to Pedersen’s book (“Phase Envelope Calculations”). Implement the suggested algorithm for an extra +10 points in your PBM final grade.
TEST IT !
Use the procedure outline above to generate the phase envelope of the fluid described in this handout (using PR OR SRK EOS). For comparison purposes, please refer to the PR- EOS prediction of this envelope that has been included below.
Use the PVTi (ECLIPSE) software in our computer lab for further testing for an extra +5 bonus points in your PBM final grade
From module 7, the PR-EOS prediction of the phase envelope is as follows;
40
From module 7, the SRK-EOS prediction of the phase envelope is as follows;
41