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.
Ver código
import numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport plotly.express as px# Definición de parámetros inicialesnx =5# numero de celdasbeta =0.001127# factor de conversiónalpha =5.614# factor de conversión# Dimensionesdx =1000# ft, dimension in x directiondy =1000# ft, dimension in y directiondz =75# ft, dimension in z direction# Propiedades del yacimientoporo =0.18# porosidad fracciónCr =0# compresibilidad de la roca psi^-1Cf =0.0000035# compresibilidad del fluido psi^-1k =15# permeabilidad, mdBo =1# factor de volumen de formación del aceite, RB/STBvis =10# viscosidad, cp # Información del pozorw =3.5# in, radio de pozos =0# skincellP =4# celda de pozoqo =-150# gasto de producción del pozo, STB/D# Tiempodt =15# paso de tiempo, díasTT =360# tiempo total de simulaciónPo =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}\]
# Cálculos adicionalesAx = dy * dz # Faces areavol = dx * dy * dz # Cells volumeTE = ((beta * Ax * k) / (vis * Bo * dx)) # East transmisibilityTW = ((beta * Ax * k) / (vis * Bo * dx)) # West transmisibilityAcum = ((vol * poro * (Cf + Cr)) / (alpha * Bo * dt))time = dtPt = np.full(nx, Po) # Presión al tiempo nPtdt = np.zeros(nx) # Presión al tiempo n + 1# Results dataframes initialization# For simplicity in Python, we'll use pandas DataFrameresults_cells = pd.DataFrame(data={'Time': [0], **{i+1: Pt[i] for i inrange(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
Ver código
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 eliminationfor j inrange(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 inreversed(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.
Ver código
# Variables para el control del ciclocurrent_time =0steps =0# Preparar vectores para el sistema tridiagonala = 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ónfor i inrange(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 = Ptdtfor j inrange(nx): Ptdt[j] = x[j] Pt[j] = Ptdt[j] current_time += dt steps +=1# Almacenar resultados results_cells = pd.concat([results_cells, pd.DataFrame([[current_time] +list(Ptdt)], columns=results_cells.columns)], ignore_index=True)"Simulación completada con {} pasos.".format(steps)
Por ú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.
Ver código
#Suponiendo que 'results_cells' es tu DataFrame y 'cellP' es el número de celda del pozoplt.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();
Ver código
#Grafio perfiles de presiónplt.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áficoplt.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áficoplt.show();
