Para la malla 1D con bloques centradas mostrada en la figura inferior, determinar la distribución de presión durante el primer año de producción. La presión inicial es de 6000 psi. Las propiedades de la roca y del fluido son Δ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. Emplea la formulación implícita y un paso de tiempo de Δt = 15 days.
En la siguiente sección definimos los datos de entrada. Importamos las librerias numpy y pandas para la creación de arreglos y dataframe. Para la visualización de los resultados importamos matplotlib.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import plotly.express as px
# Definición de parámetros iniciales
nx = 5 # numero de celdas
beta = 0.001127 # factor de conversión
alpha = 5.614 # factor de conversión
# Dimensiones
dx = 1000 # ft, dimension in x direction
dy = 1000 # ft, dimension in y direction
dz = 75 # ft, dimension in z direction
# Propiedades del yacimiento
poro = 0.18 # porosidad fracción
Cr = 0 # compresibilidad de la roca psi^-1
Cf = 0.0000035 # compresibilidad del fluido psi^-1
k = 15 # permeabilidad, md
Bo = 1 # factor de volumen de formación del aceite, RB/STB
vis = 10 # viscosidad, cp
# Información del pozo
rw = 3.5 # in, radio de pozo
s = 0 # skin
cellP = 4 # celda de pozo
qo = -150 # gasto de producción del pozo, STB/D
# Tiempo
dt = 15 # paso de tiempo, días
TT = 360 # tiempo total de simulación
Po = 6000 # psi, Icondición inicial Definimos los términos de transmisibilidad en la dirección x
\[T_{x_{i\pm1/2}} = \beta_c\frac{k_x A_x}{\mu B \Delta x}\]
Definimos el termino de acumulación
\[\frac{V_b \phi (c_o + c_{\phi})}{\alpha B_o \Delta t}\]
# Cálculos adicionales
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 = np.full(nx, Po) # Presión al tiempo n
Ptdt = np.full(nx, Po) # Presión al tiempo n + 1
# Results dataframes initialization
# For simplicity in Python, we'll use pandas DataFrame
results_cells = pd.DataFrame(data={'Time': [0], **{i+1: Pt[i] for i in range(nx)}})
cells_x = np.linspace(dx/2, (dx*nx)-(dx/2), num=nx)
results_pwf = pd.DataFrame(columns=['Time', 'Pwf'])Thomas’ algorithm
Este algoritmo es aplicable a problemas de flujo en una dimensión
def thomas(a,b,c,d,n):
cp = [0] * n
dp = [0] * n
cp[0] = c[0]/b[0]
dp[0] = d[0]/b[0]
#Forward elimination
for j in range(1,n):
m = b[j]-cp[j-1]*a[j]
cp[j] = c[j]/m
dp[j] = (d[j]-dp[j-1]*a[j])/m
x = [0] * n
x[n-1] = dp[n-1]
for j in reversed(range(0,n-1)):
x[j] = dp[j]-cp[j]*x[j+1]
return(x)Definimos un ciclo de tiempo usando while(). Como se definió anteriormente, el vector Pt representan los valores de presión al tiempo n, valores conocidos, y el vector Ptdt, los valores de presión al tiempo n + 1, que son los valores incógnita del problema. En este caso la función thomas() nos estará regresando como resultados los valores del vector Ptdt. En cada iteración reasignamos los valores de ptdt a pt para calcular al siguiente paso de tiempo.
Para este ejemplo definimos los vectores a,b,c,d empleando de la ceda 2 a 5, ya que la presión de la celda 1 se mantiene constante por el tipo de frontera.
Debemos de modificar los valores del vector b de las celdas 2 y 5 por las condiciones de frontera. De manera similar modificamos el valor de la celda 2 del vector d, para el cual agregamos el valor conocido de la celda 1.
Con estas modificaciones corremos el ciclo while y obtenemos los resultados, observamos que los valores de la celda 1 se mantiene con un valor de 6000 para todos los pasos de tiempo.
# Variables para el control del ciclo
current_time = 0
steps = 0
# Preparar vectores para el sistema tridiagonal
a = np.zeros(nx-1)
b = np.zeros(nx-1)
c = np.zeros(nx-1)
d = np.zeros(nx-1)
while current_time < TT:
# Actualizar vectores para el sistema tridiagonal basado en condiciones de simulación
for i in range(nx-1):
a[i] = TW
b[i] = -(Acum + TE + TW)
c[i] = TE
d[i] = -Acum * Pt[i+1]
# Condiciones de frontera
b[0] = -(Acum + TE + TW)
b[nx-2] = -(Acum + TW)
d[cellP-2] = d[cellP-2] - qo
d[0] = -Acum * Pt[1] - Po * TW
# Resolver el sistema para encontrar presiones en el siguiente paso de tiempo
x = thomas(a, b, c, d,nx-1)
# Actualizar presiones y tiempo
#Pt = Ptdt
for j in range(1,nx):
Ptdt[j] = x[j-1]
Pt[j] = Ptdt[j]
current_time += dt
steps += 1
# Almacenar resultados
results_cells = results_cells.append(pd.Series([current_time] + list(Ptdt), index=results_cells.columns), ignore_index=True)
"Simulación completada con {} pasos.".format(steps)'Simulación completada con 24 pasos.'print(results_cells) Time 1 2 3 4 5
0 0 6000 6000 6000 6000 6000
1 15 6000 5995 5968 5805 5964
2 30 6000 5984 5921 5654 5907
3 45 6000 5968 5867 5531 5837
4 60 6000 5948 5811 5426 5761
5 75 6000 5925 5754 5332 5682
6 90 6000 5901 5697 5246 5601
7 105 6000 5876 5642 5166 5520
8 120 6000 5850 5588 5090 5440
9 135 6000 5824 5536 5017 5362
10 150 6000 5798 5485 4947 5285
11 165 6000 5773 5435 4879 5210
12 180 6000 5748 5387 4814 5137
13 195 6000 5723 5340 4751 5065
14 210 6000 5699 5294 4689 4995
15 225 6000 5675 5250 4629 4927
16 240 6000 5652 5207 4570 4861
17 255 6000 5630 5165 4513 4797
18 270 6000 5608 5124 4458 4734
19 285 6000 5587 5084 4404 4673
20 300 6000 5566 5045 4352 4613
21 315 6000 5546 5007 4301 4555
22 330 6000 5527 4970 4251 4499
23 345 6000 5508 4934 4203 4444
24 360 6000 5489 4899 4156 4391Por último, graficamos la presión de la celda 4, donde esta ubicado el pozo, contra tiempo, y generamos los perfiles de presión a lo largo del yacimiento para cada paso de tiempo.
#Suponiendo que 'results_cells' es tu DataFrame y 'cellP' es el número de celda del pozo
plt.figure(figsize=(10, 6));
plt.plot(results_cells['Time'], results_cells[cellP], label=f'Presión en Celda {cellP}');
plt.title('Evolución de la Presión en el Tiempo');
plt.xlabel('Tiempo (días)');
plt.ylabel('Presión (psi)');
plt.legend();
plt.grid(True);
plt.show();
#Grafio perfiles de presión
plt.figure(figsize=(10, 6));
# Iterar sobre cada fila de la tabla para graficarfor index, row in results_cells.iterrows():
plt.plot(row[1:], label=f"Time: {row['Time']}");
# Configuración del gráfico[<matplotlib.lines.Line2D object at 0x00000000647B8128>]
[<matplotlib.lines.Line2D object at 0x00000000647B8400>]
[<matplotlib.lines.Line2D object at 0x00000000647B86A0>]
[<matplotlib.lines.Line2D object at 0x00000000647B8940>]
[<matplotlib.lines.Line2D object at 0x00000000647B8BE0>]
[<matplotlib.lines.Line2D object at 0x00000000647B8E80>]
[<matplotlib.lines.Line2D object at 0x0000000064C39160>]
[<matplotlib.lines.Line2D object at 0x0000000064C39400>]
[<matplotlib.lines.Line2D object at 0x0000000064C396A0>]
[<matplotlib.lines.Line2D object at 0x0000000064C39940>]
[<matplotlib.lines.Line2D object at 0x00000000647B8160>]
[<matplotlib.lines.Line2D object at 0x0000000064C39E48>]
[<matplotlib.lines.Line2D object at 0x0000000064C46128>]
[<matplotlib.lines.Line2D object at 0x0000000064C463C8>]
[<matplotlib.lines.Line2D object at 0x0000000064C46668>]
[<matplotlib.lines.Line2D object at 0x0000000064C46908>]
[<matplotlib.lines.Line2D object at 0x0000000064C46BA8>]
[<matplotlib.lines.Line2D object at 0x0000000064C46E48>]
[<matplotlib.lines.Line2D object at 0x0000000064C55128>]
[<matplotlib.lines.Line2D object at 0x0000000064C553C8>]
[<matplotlib.lines.Line2D object at 0x0000000064C55668>]
[<matplotlib.lines.Line2D object at 0x0000000064C55908>]
[<matplotlib.lines.Line2D object at 0x0000000064C55BA8>]
[<matplotlib.lines.Line2D object at 0x0000000064C55E48>]
[<matplotlib.lines.Line2D object at 0x0000000064C64128>]plt.xlabel('Celda');
plt.ylabel('Presión');
plt.title('Perfil de presión');
plt.legend(bbox_to_anchor=(1.12, 1),loc='upper right',prop={'size': 5.8});
plt.grid(True);
# Mostrar el gráfico
plt.show();
Usando la librería pyvista podemos generar una malla estructurada con las dimensiones de las celdas definidas, coloreamos las celdas de acuerdo con el valor de presión y visualizamos cada paso de tiempo añadiendo una slider widget. Por último, añadimos una línea para representar la ubicación del pozo en la celda 4.
import pyvista as pv
ni, nj, nk = 5, 1, 1
si, sj, sk = 1000, 1000, 75
xcorn = np.arange(0, (ni + 1) * si, si)
xcorn = np.repeat(xcorn, 2)
xcorn = xcorn[1:-1]
xcorn = np.tile(xcorn, 4 * nj * nk)
ycorn = np.arange(0, (nj + 1) * sj, sj)
ycorn = np.repeat(ycorn, 2)
ycorn = ycorn[1:-1]
ycorn = np.tile(ycorn, (2 * ni, 2 * nk))
ycorn = np.transpose(ycorn)
ycorn = ycorn.flatten()
zcorn = np.arange(0, (nk + 1) * sk, sk)
zcorn = np.repeat(zcorn, 2)
zcorn = zcorn[1:-1]
zcorn = np.repeat(zcorn, (4 * ni * nj))
corners = np.stack((xcorn, ycorn, zcorn))
corners = corners.transpose()
dims = np.asarray((ni, nj, nk)) + 1
grid = pv.ExplicitStructuredGrid(dims, corners)
grid = grid.compute_connectivity()
#grid.cell_data["colors"] = [5000,4000,3500,3500,2000]
#grid = grid.plot(show_edges=True)
pl = pv.Plotter()
#pl.add_mesh(grid, show_edges=True)
def create_mesh(value):
fila = int(value)
presion = list(results_cells.iloc[fila,1:])
grid.cell_data["colors"] = presion
pl.add_lines(np.array([[3500, 500, 0], [3500, 500, 350]]), color='w', width=5, label=None, name=None)
pl.add_mesh(grid, show_edges=True)
#grid.cell_data["colors"] = presion
#pl.add_mesh(sphere, name="sphere", show_edges=True)
slider = pl.add_slider_widget(
create_mesh,
[1, len(results_cells)],
fmt="%0.0f"
)
pl.show();