Hi,
I have a project in class that I'm a bit stuck on. I know the overflow problem is caused by the Z and Pi1 being zero's but don't know why they come out to zero's in the code. any help on getting this program to run?
******CODE*********
Option Explicit
Public Sub Project1()
Dim Qgsc As Double
Dim Pwh As Double
Dim Twh As Double
Dim D As Double
Dim deltaL As Double
Dim gradG As Double
Dim Sg As Double
Dim Ls As Double
Dim i As Integer
Dim n As Integer
Dim Li As Double
Dim Ti1 As Double
Dim Ti As Double
Dim Pi As Double
Dim Pi1 As Double
Dim Tavg As Double
Dim Pavg As Double
Dim Pc As Double
Dim Tc As Double
Dim Z As Double
Dim VisG As Double
Dim Re As Double
Dim Le As Double
Dim f As Double
Dim skin As Double
Dim Tol As Double
Dim Pi1guess As Double
Dim vg As Double
Dim error As Double
Qgsc = ActiveSheet.Cells(1, 2)
Pwh = ActiveSheet.Cells(2, 2)
Twh = ActiveSheet.Cells(3, 2)
D = ActiveSheet.Cells(4, 2)
deltaL = ActiveSheet.Cells(5, 2)
gradG = ActiveSheet.Cells(6, 2)
Sg = ActiveSheet.Cells(7, 2)
n = 250
For i = 2 To n
Ls = deltaL / (n - 1)
Li = Ls * (i - 1)
Ti1 = Twh + (gradG) * (Li)
Pi = Pwh
Ti = Twh
Pi1guess = 2000
Tavg = 0.5 * (Ti + Ti1)
Pavg = 0.5 * (Pi + Pi1guess)
Pc = 709.6 - 58.718 * Sg
Tc = 170.5 + 307.3 * Sg
Do While (error > Tol)
Z = Zfactor(Tavg, Pavg, Pc, Tc)
VisG = Viscosity(Sg, Pavg, Tavg, Z)
Re = 1667 * ((Qgsc * Sg) / (VisG * D))
f = Ffactor(Re, D)
Pi1 = (1.0144 * (10 ^ (-16)) * ((f * Sg * Qgsc ^ (2) * Z * Tavg * Le) / D ^ (5)) + Pi ^ (2) * Exp(skin)) ^ (1 / 2)
skin = 0.0375 * (Sg / Z * Tavg) * Ls
Le = (Ls * Exp(skin) - 1) / skin
Tol = 0.0001
error = Abs(Pi1 - Pi1guess)
If error > Tol Then
Pi1guess = Pi1
Else
Pi1 = Pi1guess
ActiveSheet.Cells(i + 1, 8) = Z
ActiveSheet.Cells(i + 1, 7) = Ti1
ActiveSheet.Cells(i + 1, 6) = Pi1
ActiveSheet.Cells(i + 1, 5) = Li
End If
Loop
vg = 5.98107 * (10 ^ (-5)) * ((Z * Ti1 / Pi1) * (Qgsc / D ^ (2)))
ActiveSheet.Cells(i + 1, 9) = vg
Pi = Pi1
Ti = Ti1
Next i
End Sub
Private Function Zfactor(Tavg As Double, Pavg As Double, Pc As Double, Tc As Double) As Double
Dim Pr As Double
Dim Tr As Double
Dim A As Double
Dim B As Double
Dim Z As Double
Pr = Pavg / Pc
Tr = Tavg / Tc
A = 0.274 * (Pr) ^ 2 / (10) ^ (0.8157 * Tr)
B = (3.52 * Pr) / (10) ^ (0.9813 * Tr)
Z = 1 - B - A
Zfactor = Z
End Function
Private Function Viscosity(Sg As Double, Pavg As Double, Tavg As Double, Z As Double) As Double
Dim DenG As Double
Dim M As Double
Dim X As Double
Dim Y As Double
Dim k As Double
Dim VisG As Double
DenG = 0.0433 * Sg * Pavg / (Z * Tavg)
M = 29 * Sg
X = 3.448 + 986.4 / Tavg + 0.01009 * M
Y = 2.4 - 0.2 * X
k = (9.379 + 0.016 * M) * (Tavg) ^ 1.5
k = k / (209.2 + 19.26 * M + Tavg)
VisG = k * 0.0001 * Exp(X * (DenG) ^ Y)
Viscosity = VisG
End Function
Private Function Ffactor(Re As Double, D As Double)
Dim B As Double
Dim A As Double
Dim f As Double
Dim error As Double
B = (37530 / Re) ^ (16)
A = (2.547 * Log(1 / ((7 / Re) ^ (0.9) + (0.27 * (error / D))))) ^ (16)
f = 8 * ((8 / Re) ^ (12) + (1 / (A + B) ^ (3 / 2))) ^ (1 / 12)
Ffactor = f
End Function
*****Project Instructions***********
(a) The project consists of calculating gas velocity in a production tubing string knowing the
following information: gas flow rate at standard conditions (Qg,sc ), wellhead pressure ( pwh ),
wellhead temperature (Twh ), pipe diameter ( D ), well total depth (∆L ), geothermal gradient (
gradG ), and gas specific gravity ( g
γ ). These are going to be the input for your code.
The suggested pseudo-code is the following
(1) Take the total depth of the well and split it in segments of constant length as shown in the
below figure
Where, Ls is the length of each segment that is calculated as,
−1
∆
=
n
L
Ls
∆L is the well depth and “n”is the total number of segments or intervals. Note that a node
“i” is located at the edge of each interval. It is at those nodes where you will calculate the
pressure, fluid properties, temperature and gas velocity. The depth of each node is
calculated as,
L = L ( ) i −1 i s
(2) Compute the temperature downstream and upstream temperature,
( )( ) Ti+1 = Twh + gradG Li
(3) Make pi = pwh and Ti = Twh (4) The upstream pressure ( pi+1
) is calculated as follows,
a. Guess upstream pressure values i guess p +1
b. Calculate the average temperature and pressure
( )1
5.0 T = Ti +Ti+
( ) i i guess p p p ,1
5.0 = + +
Ti: downstream temperature, Ti+1: upstream temperature
Pi: downstream pressure, pi+1: upstream pressure
c. Compute deviation factor, gas density and gas viscosity using the average
pressure and temperature. Look for your old functions that we have worked in
quizzes and add it your project. If you don’t have it, you will need to create it
again.
d. Calculate Reynolds number
D
Q
g
g sc g
µ
γ , Re =1667
e. Write a function that calculate friction factor
16
Re
37530
B =
16
0 9
0 27
Re
7
1
2 547 ln
+
=
D
ε
.
A .
.
( )
112
3 2
12
Re
8 1
8
+
+
=
N A B
f
f. Compute the upstream pressure
[ ] ( ) ( )
5
2
2 2 16 ,
1
exp .1 0114 10
D
f Q ZT L
p p S
g g sc e
i i
γ −
+ − = ×
S
g
L
ZT
S
γ
= .0 0375
( ) [ ]
S
L S
L
S
e
exp −1
=
Units: Qg,sc = SCFD, D =ft,T = deg R, L Le
, =ft, i i
p , p +1 =psia µ g
= cP,
Pi: downstream pressure, pi+1: upstream pressure calculated
g. Compute the error
= pi+ ,1 actual − pi+1guess ε h. For error>tolerance, make i guess i actual p p +1 = + ,1
. Otherwise, upstream pressure has
been calculated successfully ( i i guess p p +1 = +1
)
(5) Once the upstream pressure is calculated, Estimate the gas velocity at the node
2
,
1
5 1
.5 98107 10
D
Q
p
zT
v
g sc
i
i
g
+
− + = × T R p psia D in Q scfd
s
ft
vg = , = deg , = , = ,
g,sc =
Deviation factor is evaluated at the node pressure (pi+1) and temperature (Ti+1)
(6) Print the results
(7) Make pi = pi+1
and Ti = Ti+1
(8) Repeat step 3 through 6 until reaching i=n
Your subroutine must print out for each node: depth, pressure, temperature, deviation factor, and
gas velocity.
(b) Solve the following cases:
A 4,000 ftvertical gas well requires a wellhead pressure of at least 100 psia to enable flow of the
0.85 gravity gas from the wellhead through the surface facility. The geothermal gradient is 2.875
deg F/100 ft while the wellhead temperature is 82 deg F.
(a) Determine the bottom-hole pressure for a gas flow rate of 10 MMscfd if it is flowing through
a 4.5 in ID casing
(b) Repeat the part (a) assuming that the gas flows through a tubing with an ID of 1.83 in.
(b) Repeat the part (a) assuming that the gas flows through the casing-tubing annulus. The casing
ID is 4.5” while the tubing OD is 2.1”. In this case, you can replace the diameter in the
calculation by an equivalent diameter given by,
2
2
2 Deq = D1 − D where D1 is the outer diameter and D2 is inner diameter
Prepare a short report with the solution for these problems. Add to the reports the following
plots: pressure (x-axis) vs. depth (y-axis), gas velocity (x-axis) vs. depth (y-axis), and Z (x-axis)
vs. depth (y-axis). Add a short manual about how to use your program.
Bookmarks