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.zeros(nx) # 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.
# Variables para el control del ciclo
current_time = 0
steps = 0
# Preparar vectores para el sistema tridiagonal
a = np.zeros(nx)
b = np.zeros(nx)
c = np.zeros(nx)
d = np.zeros(nx)
while current_time < TT:
# Actualizar vectores para el sistema tridiagonal basado en condiciones de simulación
for i in range(nx):
a[i] = TW
b[i] = -(Acum + TE + TW)
c[i] = TE
d[i] = -Acum * Pt[i]
# Condiciones de frontera
b[0] = -(Acum + TE)
b[nx-1] = -(Acum + TW)
d[cellP-1] = d[cellP-1] - qo
# Resolver el sistema para encontrar presiones en el siguiente paso de tiempo
x = thomas(a, b, c, d,nx)
# Actualizar presiones y tiempo
#Pt = Ptdt
for j in range(nx):
Ptdt[j] = x[j]
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.0 6000.000000 6000.000000 6000.000000 6000.000000 6000.000000
1 15.0 5999.083033 5995.025001 5968.950142 5805.464338 5964.144153
2 30.0 5996.195345 5983.783360 5921.732208 5654.923245 5907.032509
3 45.0 5990.504508 5966.184242 5867.445950 5531.704558 5837.827409
4 60.0 5981.559802 5944.207760 5810.414301 5426.201198 5761.283606
5 75.0 5969.416982 5918.156413 5752.525214 5332.539726 5682.028332
6 90.0 5954.230874 5888.870205 5694.595817 5246.279591 5601.690179
7 105.0 5936.013228 5856.412885 5637.023968 5165.487994 5520.728591
8 120.0 5915.044673 5822.307008 5580.461228 5088.403482 5440.450275
9 135.0 5891.301004 5786.421250 5524.087967 5014.316343 5361.540104
10 150.0 5864.800198 5748.853145 5468.512928 4942.615033 5283.885362
11 165.0 5835.517418 5709.467939 5412.895074 4872.455717 5207.330518
12 180.0 5804.427827 5669.130875 5357.393367 4803.994559 5132.720040
13 195.0 5771.548769 5627.936059 5302.595144 4736.485890 5059.100805
14 210.0 5736.728776 5585.061695 5247.796920 4670.656770 4987.422506
15 225.0 5700.232032 5541.941123 5193.093436 4605.682602 4916.717472
16 240.0 5662.652104 5497.369235 5138.998697 4541.645746 4847.000884
17 255.0 5623.374651 5452.438421 5084.295213 4478.481821 4779.076561
18 270.0 5583.148138 5406.783980 5030.316621 4416.273816 4712.144112
19 285.0 5541.951156 5360.289765 4976.338028 4354.905585 4646.182132
20 300.0 5499.073995 5313.530732 4922.338028 4293.664618 4581.059293
21 315.0 5456.066772 5266.066050 4868.359436 4233.266855 4516.907555
22 330.0 5412.222573 5218.485853 4814.475583 4173.593351 4452.889308
23 345.0 5367.404476 5170.047272 4760.475583 4114.031931 4389.707405
24 360.0 5321.750035 5121.496607 4706.591729 4055.328896 4327.499400Por ú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 0x0000000063282898>]
[<matplotlib.lines.Line2D object at 0x0000000063282BA8>]
[<matplotlib.lines.Line2D object at 0x0000000063282E48>]
[<matplotlib.lines.Line2D object at 0x000000006323AE48>]
[<matplotlib.lines.Line2D object at 0x00000000632923C8>]
[<matplotlib.lines.Line2D object at 0x0000000063292668>]
[<matplotlib.lines.Line2D object at 0x0000000063292908>]
[<matplotlib.lines.Line2D object at 0x0000000063292BA8>]
[<matplotlib.lines.Line2D object at 0x0000000063292E48>]
[<matplotlib.lines.Line2D object at 0x0000000063282E80>]
[<matplotlib.lines.Line2D object at 0x00000000632828D0>]
[<matplotlib.lines.Line2D object at 0x00000000632A1630>]
[<matplotlib.lines.Line2D object at 0x00000000632A18D0>]
[<matplotlib.lines.Line2D object at 0x00000000632A1B70>]
[<matplotlib.lines.Line2D object at 0x00000000632A1E10>]
[<matplotlib.lines.Line2D object at 0x0000000063292E80>]
[<matplotlib.lines.Line2D object at 0x00000000632AF390>]
[<matplotlib.lines.Line2D object at 0x00000000632AF630>]
[<matplotlib.lines.Line2D object at 0x00000000632AF8D0>]
[<matplotlib.lines.Line2D object at 0x00000000632AFB70>]
[<matplotlib.lines.Line2D object at 0x00000000632AFE10>]
[<matplotlib.lines.Line2D object at 0x00000000632A1E48>]
[<matplotlib.lines.Line2D object at 0x00000000632BF390>]
[<matplotlib.lines.Line2D object at 0x00000000632BF630>]
[<matplotlib.lines.Line2D object at 0x00000000632BF8D0>]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();
import pyvista as pv
ni, nj, nk = nx, 1, 1
si, sj, sk = dx, dy, dz
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.plot(show_edges=True)