Simulación de Yacimientos 1D en Coordenadas Radiales

Ejemplo 7.12 de Abou-Kassem et al. (2006) resuelto en R

Yacimientos
Simulación
R

Simulación numérica de un yacimiento en coordenadas radiales con 5 bloques, un pozo productor y fronteras de no flujo. Resolvemos el sistema con el algoritmo de Thomas y obtenemos la distribución de presión a 1 y 3 días.

Author

Rigoberto Chandomi

Published

June 10, 2026

Descripción del Problema

Un pozo de agua de 0.5 pies de diámetro está ubicado en un espaciamiento de 20 acres. El espesor del yacimiento, la permeabilidad horizontal y la porosidad son 30 ft, 150 md y 0.23, respectivamente. El fluido tiene un factor volumétrico de 1 RB/B, compresibilidad de 1×10⁻⁵ psi⁻¹ y viscosidad de 0.5 cp. Las fronteras externas son de no flujo. El pozo tiene terminación abierta y produce a un gasto de 2,000 B/D. La presión inicial del yacimiento es 4,000 psia.

El yacimiento se simula usando 5 bloques en la dirección radial. Encontrar la distribución de presión después de 1 día y 3 días.

Esquema del yacimiento radial con 5 bloques

Video

Referencia

Este ejemplo corresponde al Ejemplo 7.12 del libro Reservoir Simulation: A Basic Approach de Abou-Kassem, Farouq Ali & Rafiq (2006), Gulf Publishing Company.

Datos de Entrada

Ver código 
# Simulación 1D radial, una fase, yacimiento homogéneo
library(ggplot2)
library(reshape2)

# Datos del grid
nx <- 5            # Número de celdas
beta <- 0.001127   # Factor de conversión de Darcy
alpha <- 5.614     # Factor de conversión volumétrico

re <- 526.6040     # Radio externo del yacimiento (ft)
h <- 30            # Espesor (ft)

# Datos de roca y fluidos
poro <- 0.23       # Porosidad (fracción)
Cr <- 0            # Compresibilidad de la roca (psi⁻¹)
Cf <- 0.00001      # Compresibilidad del fluido (psi⁻¹)
k <- 150           # Permeabilidad (md)
Bo <- 1            # Factor volumétrico (RB/STB)
vis <- 0.5         # Viscosidad (cp)

# Información del pozo
rw <- 3            # Radio del pozo (in)
s <- 0             # Factor de daño (skin)
cellP <- 1         # Celda del pozo
qo <- -2000        # Gasto del pozo (STB/D), negativo = producción

# Control de tiempo
dt <- 1            # Paso de tiempo (días)
TT <- 3            # Tiempo total de simulación (días)

# Condición inicial
Po <- 4000         # Presión inicial (psia)

Cálculo del Grid Radial

En coordenadas radiales, los nodos del grid no están igualmente espaciados como en el caso cartesiano. Se utiliza una distribución logarítmica para que las celdas cerca del pozo sean más pequeñas (mayor resolución donde los gradientes de presión son mayores).

Factor geométrico:

\[\alpha_{lg} = \left(\frac{r_e}{r_w}\right)^{1/n_r}\]

Radio del primer nodo:

\[r_1 = \frac{\alpha_{lg} \ln(\alpha_{lg})}{\alpha_{lg} - 1} \cdot r_w\]

Radio de los nodos siguientes:

\[r_{i+1} = \alpha_{lg} \cdot r_i\]

Radio de las caras de los bloques:

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

Ver código 
# Factor geométrico
al2 <- (re / (rw / 12))^(1 / nx)

# Radio de los nodos
rn <- numeric(nx)
rn[1] <- (al2 / (al2 - 1)) * log(al2) * (rw / 12)

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

# Radio de las caras de los bloques
rf <- numeric(nx)
for (i in 1:(nx - 1)) {
  rf[i] <- (rn[i + 1] - rn[i]) / log(rn[i + 1] / rn[i])
}
rf[nx] <- re

# Volumen de las celdas
vol <- numeric(nx)
vol[1] <- pi * (rf[1]^2 - (rw / 12)^2) * h
for (i in 2:nx) {
  vol[i] <- pi * (rf[i]^2 - rf[i - 1]^2) * h
}

cat("Radios de nodos (ft):", round(rn, 2), "\n")
Radios de nodos (ft): 0.49 2.26 10.43 48.18 222.61
Ver código 
cat("Radios de caras (ft):", round(rf, 2), "\n")
Radios de caras (ft): 1.16 5.34 24.66 113.97 526.6
Ver código 
cat("Volúmenes (ft³):", round(vol, 1), "\n")
Volúmenes (ft³): 119.9 2559.5 54647.7 1166781 24911905

Transmisibilidad y Acumulación

Con la geometría definida, calculamos los términos de transmisibilidad y acumulación que forman el sistema de ecuaciones lineales.

Factor geométrico de transmisibilidad:

\[G_r = \frac{2\pi \beta k h}{\ln(\alpha_{lg})}\]

Transmisibilidad en dirección radial:

\[T_r = G_r \cdot \frac{1}{\mu B_o}\]

Término de acumulación:

\[A_i = \frac{V_i \cdot \phi \cdot (c_f + c_r)}{\alpha \cdot B_o \cdot \Delta t}\]

Ver código 
# Factor geométrico de transmisibilidad
Gr <- (2 * pi * beta * k * h) / log(al2)

# Transmisibilidades Este y Oeste
TE <- Gr * (1 / (vis * Bo))
TW <- Gr * (1 / (vis * Bo))

# Término de acumulación
Acum <- (vol * poro * (Cf + Cr)) / (alpha * Bo * dt)

# Inicialización
time <- dt
Pt <- rep(Po, nx)       # Presión en tiempo n
Ptdt <- rep(0, nx)      # Presión en tiempo n+1

# Factor geométrico del pozo
FG <- (2 * pi * beta * k * h) / (log(rn[1] / (rw / 12)) + s)

# Dataframes de resultados
results_cells <- data.frame(0, t(Pt))
colnames(results_cells) <- c("Time", 1:nx)

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

cat("Transmisibilidad:", round(TE, 4), "STB/(D·psi)\n")
Transmisibilidad: 41.6389 STB/(D·psi)
Ver código 
cat("Factor geométrico del pozo:", round(FG, 4), "\n")
Factor geométrico del pozo: 47.5952

Algoritmo de Thomas

Para resolver el sistema tridiagonal de ecuaciones lineales en cada paso de tiempo, utilizamos el algoritmo de Thomas (también conocido como TDMA — TriDiagonal Matrix Algorithm).

Ver código 
thomas <- function(a, b, c, d, x, n) {
  # Resuelve un sistema tridiagonal Ax = d
  # a = vector subdiagonal
  # b = vector diagonal principal
  # c = vector superdiagonal
  # d = vector del lado derecho
  # x = vector solución
  # n = número de ecuaciones
  
  # Eliminación hacia adelante
  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]
  }
  
  # Sustitución hacia atrás
  x[n] <- d[n] / b[n]
  for (i in (n - 1):1) {
    x[i] <- (d[i] - c[i] * x[i + 1]) / b[i]
  }
  
  return(x)
}

Ciclo de Tiempo y Solución

Con la geometría, las transmisibilidades y el algoritmo de Thomas definidos, ejecutamos el ciclo de tiempo. En cada paso:

  1. Armamos los vectores del sistema tridiagonal (a, b, c, d)
  2. Aplicamos las condiciones de frontera (no flujo en ambos extremos)
  3. Incluimos el pozo como término fuente
  4. Resolvemos con Thomas
  5. Actualizamos la presión
Ver código 
# Ciclo de tiempo
while (time <= TT) {
  
  # Vectores del sistema tridiagonal
  aa <- rep(TW, nx)                    # Subdiagonal
  bb <- -(TW + TE + Acum)             # Diagonal principal
  cc <- rep(TE, nx)                    # Superdiagonal
  dd <- -Acum * Pt                     # Lado derecho
  
  # Condiciones de frontera: No flujo
  bb[1] <- -(TE + Acum[1])            # Frontera oeste (pozo)
  bb[nx] <- -(TW + Acum[nx])          # Frontera este (externa)
  
  # Pozo: término fuente
  dd[cellP] <- dd[cellP] - qo
  
  # Resolver sistema tridiagonal
  Ptdt <- thomas(aa, bb, cc, dd, Ptdt, nx)
  
  # Actualizar presión
  Pt <- Ptdt
  
  # Calcular presión de fondo fluyente (Pwf)
  pwf <- qo / (FG / (Bo * vis)) + Pt[cellP]
  
  # Guardar resultados
  results_cells <- rbind(results_cells, c(time, Ptdt))
  results_pwf <- rbind(results_pwf, c(time, pwf))
  
  time <- time + dt
}

cat("Simulación completada.\n")
Simulación completada.

Resultados

Distribución de Presión

Ver código 
# Preparar datos para graficar
results_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$Radius <- rep(rn, each = nrow(results_cells))

# Gráfico de presión vs radio
ggplot(results_cells_time[results_cells_time$Time != 0, ], 
       aes(Radius, Pressure, color = Time)) + 
  geom_line(linewidth = 1.2) +
  geom_point(size = 3) +
  xlab("Radio (ft)") + 
  ylab("Presión (psia)") +
  labs(color = "Tiempo (días)") +
  theme_minimal() +
  theme(text = element_text(size = 14))
Figure 1: Distribución de presión radial a 1, 2 y 3 días. La caída de presión es mayor cerca del pozo y se atenúa hacia el exterior.

Tabla de Resultados

Ver código 
# Tabla de presión por celda y tiempo
colnames(results_cells) <- c("Tiempo (días)", paste0("Celda ", 1:nx))
knitr::kable(results_cells, digits = 2, 
             caption = "Presión (psia) en cada celda para cada paso de tiempo")
Presión (psia) en cada celda para cada paso de tiempo
Tiempo (días) Celda 1 Celda 2 Celda 3 Celda 4 Celda 5
0 4000.00 4000.00 4000.00 4000.00 4000.00
1 3626.28 3674.31 3722.34 3770.21 3815.45
2 3438.93 3486.96 3534.99 3582.91 3628.69
3 3252.14 3300.17 3348.20 3396.13 3441.91

Observaciones

  • La caída de presión es mayor cerca del pozo (celda 1) debido a la convergencia del flujo radial — el área de flujo es menor cerca del pozo.
  • La frontera externa (celda 5) apenas siente el efecto del pozo después de 3 días — la señal de presión no ha llegado completamente.
  • La distribución logarítmica del grid captura adecuadamente el gradiente de presión cerca del pozo sin necesidad de muchas celdas.

Referencias

  • Abou-Kassem, J., Farouq Ali, S.M. & Islam, M.R. (2006). Reservoir Simulation: A Basic Approach. Gulf Publishing Company. Ejemplo 7.12.
  • Aziz, K. & Settari, A. (1979). Petroleum Reservoir Simulation. Applied Science Publishers.
🛒
Propiedades De Los Fluidos Del Yacimiento

Disponible en Mercado Libre

Ver oferta →