For the 1D, block-centered grid shown below, determine the pressure distribution during the first year of production. The initial reservoir pressure is 6000 psia. The rock and fluid properties for this problem are Δx = 1000 ft, Δy = 1000 ft and Δz = 75 ft, Bo = 1 RB/STB, co = 3.5 x 10^-6 psi-1 , kx = 15 md, ϕ = 0.18, μo = 10 cp. Using the implicit formulation and a timestep size of Δt = 15 days.
Define imput data
library(ggplot2)
library(reshape2)
library(plotly)
library(scales)
library(DT)
nx <- 5 #number of cells
beta <- 0.001127 #convertion factor
alpha <- 5.614 #convertion factor
dx <- 1000 #ft, dimension in x direction
dy <- 1000 #ft, dimension in y direction
dz <- 75 #ft, dimension in z direction
poro <- 0.18 #porosity, fraction
Cr <- 0 #rock Compressibility psi^-1
Cf <- 0.0000035 #fluid Compressibility psi^-1
k <- 15 #permeability, md
Bo <- 1 #Formation volume factor, RB/STB
vis <- 10 #viscosity,cp
#Well information
rw <- 3.5 #in, well radius
s <- 0 #skin
cellP <- 4 #well cell
qo <- -150 #well rate, STB/D
dt <- 15 #time step, days
TT <- 360 #total simulation time
Po <- 6000 #psi, Initial condition
Transmissibilities in the x-direction
\[T_{x_{i\pm1/2}} = \beta_c\frac{k_x A_x}{\mu B \Delta x}\] Fluid accumulation term
\[\frac{V_b \phi (c_o + c_{\phi})}{\alpha B_o \Delta t}\]
Ax <- dy * dz #Faces area
vol <- dx * dy * dz #Cells volume
TE <- ((beta * Ax * k) / (vis * Bo * dx)) # East transmisibility
TW <- ((beta * Ax * k) / (vis * Bo * dx)) # West transmisibility
Acum <- ((vol * poro * (Cf + Cr)) / (alpha * Bo * dt))
time <- dt
Pt <- rep(Po, nx) #Pressure at time n
Ptdt <- rep(0, nx) #Pressure at time n + 1
#geometry factor
req <- 0.14 * ((dx) **2 + (dy) ** 2) ** 0.5
FG <- (2 * 3.1416 * beta * k * dz) / (log(req / (rw / 12)) + s)
#Results data.frame
results_cells <- data.frame(0, t(Pt))
colnames(results_cells) <- c("Time", 1:nx)
cells_x <- seq(dx/2, (dx*nx)-(dx/2), length.out = nx)
results_pwf <- data.frame(Time = numeric(), Pwf = numeric())
Thomas’ algorithm
This algorithm is applicable for a reservoir where flow takes place in the x-direction in rectangular flow problems
thomas <- function(a,b,c,d,x,n){
# Subroutine to solve a tridiagonal system
# a = subdiagonal vector
# B = diagonal vector
# c = superdiagonal vector
# d = right hand side vector
# x = solution vector
# n = number of diagonal vector elements
#Forward elimination
for(i in 2:n){
b[i] <- b[i]-a[i]*c[i-1]/b[i-1]
d[i] <- d[i]-a[i]*d[i-1]/b[i-1]
}
#Back substitution
x[n] <- d[n]/b[n]
for(i in (nx-1):1){
x[i] <- (d[i]-c[i]*x[i+1])/b[i]
}
return(x)
}
Time loop
while( time <= TT ){
aa <- rep(TW, nx)
bb <- rep(-(TW + TE + Acum), nx)
cc <- rep(TE, nx)
dd <- rep(-Acum, nx)*Pt
#Boundary condition, No flow
#West
bb[1] <- -(TE + Acum)
#East
bb[nx] <- -(TW + Acum)
#Well
dd[cellP] = dd[cellP] - qo
Ptdt <- thomas(aa,bb,cc,dd,Ptdt,nx)
Pt <- Ptdt
pwf = qo / (FG / (Bo * vis)) + Pt[cellP]
results_cells <- rbind(results_cells, c(time,Ptdt))
results_pwf <-rbind(results_pwf, c(time, pwf))
time = time + dt
}
print(results_cells)
Time 1 2 3 4 5
1 0 6000.000 6000.000 6000.000 6000.000 6000.000
2 15 5999.083 5995.025 5968.950 5805.464 5964.144
3 30 5996.291 5983.932 5922.483 5655.392 5907.236
4 45 5990.910 5967.100 5868.802 5532.939 5838.248
5 60 5982.518 5945.380 5812.120 5428.013 5762.635
6 75 5970.932 5919.654 5754.523 5334.501 5683.724
7 90 5956.140 5890.678 5696.986 5248.661 5603.535
8 105 5938.245 5859.054 5639.909 5168.168 5523.291
9 120 5917.416 5825.238 5583.411 5091.553 5443.715
10 135 5893.854 5789.577 5527.480 5017.865 5365.224
11 150 5867.769 5752.329 5472.052 4946.474 5288.042
12 165 5839.371 5713.695 5417.043 4876.952 5212.272
13 180 5808.859 5673.830 5362.376 4808.993 5137.942
14 195 5776.419 5632.859 5307.981 4742.375 5065.033
15 210 5742.223 5590.886 5253.798 4676.927 4993.499
16 225 5706.424 5548.000 5199.782 4612.516 4923.278
17 240 5669.166 5504.277 5145.894 4549.031 4854.299
18 255 5630.574 5459.786 5092.105 4486.382 4786.486
19 270 5590.765 5414.589 5038.392 4424.489 4719.765
20 285 5549.842 5368.741 4984.738 4363.286 4654.060
21 300 5507.901 5322.292 4931.127 4302.711 4589.301
22 315 5465.028 5275.291 4877.551 4242.710 4525.419
23 330 5421.300 5227.780 4824.001 4183.236 4462.350
24 345 5376.787 5179.799 4770.471 4124.244 4400.032
25 360 5331.556 5131.385 4716.956 4065.695 4338.409results_cells_time <-reshape2::melt(results_cells,id.vars=c("Time"))
colnames(results_cells_time) <- c("Time", "Cell", "Pressure")
results_cells_time$Time <- as.factor(results_cells_time$Time)
results_cells_time$Cell <- as.numeric(results_cells_time$Cell)
ggplot(results_cells_time, aes(Cell, Pressure, color = Time)) +
geom_line()
colnames(results_pwf) <-c("Time", "Pwf")
ggplot(results_pwf, aes(Time, Pwf)) +
geom_point()
