Reservoir simulation example (7.12) from AAbou-Kassem, J., Farouq, S. M. & Rafiq, M. (2006) Reservoir simulation. A basic approach. Gulf Publising Company

11-07-2021

A 0.5-ft-diameter wáter well is located in 20-acre spacing. The reservoir thickness, horizontal permeability, and porosity are 30 ft, 150 md, and 0.23, respectively. The flowing fluid has FVF, compressibility and viscosity of 1 RB/B, 1 x 10^-5, and 0.5 cP, respectively. The external boundaries are no-Flow boundaries. The well has open-well completion and is placed on production at a rate of 2000 B/D. Initial reseroir pressure is 4000 psia. The reservoir can be simulated using five gridblocks in the radial direction.Find the pressure distribution in the reservoir after 1 day and 3 days.

#Reservoir simulation example. 1D radial coordinates, one phase only and homogeneous
library(plotly)
library(reshape2)

#input data

nx <- 5   #number of cells
beta <- 0.001127  #convertion factor
alpha <- 5.614  #convertion factor

re <- 526.6040 #Reservoir external radius, ft
h <- 30 #ft, dimension in z direction

#Rock and fluid data
poro <- 0.23 #porosity, fraction
Cr <- 0 #rock Compressibility psi^-1
Cf <- 0.00001 #fluid Compressibility  psi^-1
k <- 150 #permeability, md
Bo <- 1 #Formation volume factor, RB/STB
vis <- 0.5 #viscosity,cp

#Well information
rw <- 3 #in, well radius
s <- 0 #skin
cellP <- 1 #well cell
qo <- -2000 #well rate, STB/D

dt <- 1 #time step, days
TT <- 3 #total simulation time

Po <- 4000 #psi, Initial condition


Grid node and its data are calculated using the following equation:

Geometrical factor $\alpha_{lg} = \left(\frac{re}{rw}\right)^{1/n_r}$ First node radius $r_1 = [\alpha_{lg}ln(\alpha_{lg})/(\alpha_{lg}-1)]r_w$ Nodes radius $r_{i+1} = \alpha_{lg}r_i$

$r_{i+1}^L = \frac{r_{i+1}-r_{i}}{ln(r_{i+1}/r_{i})}$

#location gridblocks
al2 <-  (re / (rw / 12))^(1 / nx) #alpha, geometrical factor
rn <- 0
rn[1]  <-  (al2 / (al2 - 1)) * log(al2) * (rw / 12)

for (i in 2:nx){
rn[i] <-  al2 * rn[i - 1]
}

rf <- 0
for (i in 1:nx-1){
rf[i] <-  (rn[i + 1] - rn[i]) / log(rn[i + 1] / rn[i])
}

rf[nx]  <-  re

#Cells volume
vol <- 0
vol[1] <-  3.1416 * (rf[1] ** 2 - (rw / 12) ** 2) * h
for(i in 2:nx){
vol[i] = 3.1416 * (((rf[i]) ** 2) - ((rf[i-1]) ** 2)) * h
}

#Faces area
Ax = 2 * 3.14 * rf[1] * h


The past data is used to estimate transmissibility and acumulation terms.

Geometrical factor

$G_{ri} = \frac {2 \pi \beta kh}{ln(\alpha_{lg})}$

Transmissibility in the $$r$$ direction

$T_{ri} = G_{ri} \left( \frac{1}{\mu B} \right)$

#Geometrical factor
Gr <- (2*pi*beta*k*h)/(log(al2))
TE <-  Gr*(1/(vis * Bo)) #  East transmisibility
TW <-  Gr*(1/(vis * Bo)) #  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 well factor

FG <- (2 * 3.1416 * beta * k * h) / (log(rn[1] / (rw / 12)) + s)

#Results data.frame
results_cells <- data.frame(0, t(Pt))
colnames(results_cells) <- c("Time", 1:nx)
cells_x <- rn #seq(dx/2, (dx*nx)-(dx/2), length.out = nx)

results_pwf <- data.frame(Time = numeric(), Pwf = numeric())


We have to define thomas algorithm function to solve tridiagonal matrix.

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)
}


Once we have defined the geometry and the equations terms, we can begin the time loop, where we have to define de vector for thomas algoriths and get the pressure in the next time

#Time loop
while( time <= TT ){

aa <- rep(TW, nx)
bb <- -(TW + TE + Acum)
cc <- rep(TE, nx)
dd <- -Acum*Pt

#Boundary condition, No flow
#West
bb[1] <- -(TE + Acum[1])
#East
bb[nx] <- -(TW + Acum[nx])

#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
}


At the end we just generate some plots to visualize the results

#Plot
options(repr.plot.width=16, repr.plot.height=8)
results_cells_time <- reshape2::melt(results_cells,id.vars=c("Time"))
results_cells_time$Time <- as.factor(results_cells_time$Time)
results_cells_time\$Radius <- rep(rn, each = nrow(results_cells))

ggplot(results_cells_time, aes(Radius, Pressure, color = Time)) +
geom_line() +
xlab("Radius (ft)") + ylab("Pressure (psia)") +
theme(text = element_text(size=20))

print(results_cells)

  Time        1        2        3        4        5
1    0 4000.000 4000.000 4000.000 4000.000 4000.000
2    1 3626.282 3674.314 3722.337 3770.211 3815.447
3    2 3438.930 3486.961 3534.988 3582.915 3628.691
4    3 3252.142 3300.173 3348.201 3396.127 3441.910
print(results_cells_time)

   Time      Radius Pressure
1     0   0.4883173 4000.000
2     1   0.4883173 3626.282
3     2   0.4883173 3438.930
4     3   0.4883173 3252.142
5     0   2.2563731 4000.000
6     1   2.2563731 3674.314
7     2   2.2563731 3486.961
8     3   2.2563731 3300.173
9     0  10.4260485 4000.000
10    1  10.4260485 3722.337
11    2  10.4260485 3534.988
12    3  10.4260485 3348.201
13    0  48.1757596 4000.000
14    1  48.1757596 3770.211
15    2  48.1757596 3582.915
16    3  48.1757596 3396.127
17    0 222.6062736 4000.000
18    1 222.6062736 3815.447
19    2 222.6062736 3628.691
20    3 222.6062736 3441.910

Reference

Abou-Kassem, J., Farouq, S. M. & Rafiq, M. (2006) Reservoir simulation. A basic approach. Gulf Publising Compana.

### Citation

Vazquez (2021, Nov. 7). Chato Solutions: 1D radial Reservoir Simulation. Retrieved from https://www.chatosolutions.com/posts/2021-11-07-radial1d/
@misc{vazquez20211d,
}