Del gasto y las propiedades del yacimiento a los perfiles de presión a 1, 12 y 24 horas
Cuando un pozo comienza a producir a gasto constante, la caída de presión se extiende gradualmente hacia el yacimiento. En este ejemplo calcularemos la presión a diferentes distancias del pozo después de 1, 12 y 24 horas.
Resolveremos el ejemplo siguiendo un flujo sencillo de Python: datos → integral exponencial → presión → tabla → gráficas.
Un pozo produce aceite a un gasto constante de 300 STB/día, con las siguientes propiedades:
| Variable | Valor | Unidad |
|---|---|---|
| Gasto de aceite, q | 300 | STB/día |
| Espesor, h | 15 | ft |
| Permeabilidad, k | 60 | mD |
| Porosidad, φ | 0.15 | fracción |
| Viscosidad, μ | 1.5 | cP |
| Factor de volumen, B | 1.25 | rb/STB |
| Presión inicial, pi | 4000 | psi |
| Radio del pozo, rw | 0.25 | ft |
| Compresibilidad total, ct | 0.000012 | psi⁻¹ |
Evaluaremos los radios 0.25, 5, 10, 50, 100, 500, 1000, 1500, 2000 y 2500 ft. El radio se mide desde el centro del pozo.
El factor de volumen B se expresa en rb/STB y convierte el gasto a condiciones de tanque en gasto a condiciones de yacimiento.
Usamos la solución de fuente lineal para flujo radial transitorio a gasto constante:
\[ u = \frac{948\,\phi\,\mu\,c_t\,r^2}{k\,t} \]
\[ p(r,t)=p_i+\frac{70.6\,q\,\mu\,B}{k\,h}\operatorname{Ei}(-u) \]
Las constantes corresponden a las unidades de la tabla, con r en ft y t en horas, aunque q esté en STB/día. No convertimos el gasto a STB/h.
Para evitar confusiones con el signo emplearemos \(E_1(u)=-\operatorname{Ei}(-u)\):
\[ p(r,t)=p_i-\frac{70.6\,q\,\mu\,B}{k\,h}E_1(u) \]
La función exp1 de SciPy evalúa \(E_1\) directamente; la relación entre ambas funciones está documentada en SciPy.
Suponemos un yacimiento homogéneo, isotrópico, de espesor constante y comportamiento infinito durante los tiempos evaluados, con un fluido monofásico ligeramente compresible y propiedades constantes. No incluimos daño, almacenamiento del pozo ni efectos de fronteras. La evaluación en rw es la aproximación de fuente lineal utilizada en el ejercicio.
Necesitamos NumPy, pandas, Matplotlib y SciPy. Si no están instalados, ejecuta pip install numpy pandas matplotlib scipy en la terminal. Después ejecuta los fragmentos en orden dentro del mismo cuaderno o sesión de Python.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.special import exp1
q, h, k = 300, 15, 60
phi, mu, B = 0.15, 1.5, 1.25
pi, rw, ct = 4000, 0.25, 0.000012
radios = np.array([rw, 5, 10, 50, 100, 500, 1000, 1500, 2000, 2500])
tiempos = [1, 12, 24]Una función pequeña permite repetir el cálculo para cualquier radio y tiempo positivos del ejemplo.
def presion(r, t):
u = 948 * phi * mu * ct * r**2 / (k * t)
caida = 70.6 * q * mu * B / (k * h) * exp1(u)
return pi - caida
print(f"Presión a 0.25 ft y 1 h: {presion(rw, 1):.2f} psi")Presión a 0.25 ft y 1 h: 3459.13 psiEn este primer punto, \(u=2.66625\times10^{-6}\), \(E_1(u)\approx12.257625\) y el multiplicador de la integral es 44.125 psi. Por tanto, la presión es 3459.13 psi.
tabla = pd.DataFrame({"Radio (ft)": radios})
for t in tiempos:
tabla[f"Presión a {t} h (psi)"] = presion(radios, t)
tabla.round(2)| Radio (ft) | Presión a 1 h (psi) | Presión a 12 h (psi) | Presión a 24 h (psi) | |
|---|---|---|---|---|
| 0 | 0.25 | 3459.13 | 3349.49 | 3318.90 |
| 1 | 5.00 | 3723.46 | 3613.86 | 3583.27 |
| 2 | 10.00 | 3784.49 | 3675.01 | 3644.44 |
| 3 | 50.00 | 3922.13 | 3816.67 | 3786.28 |
| 4 | 100.00 | 3970.89 | 3876.68 | 3846.87 |
| 5 | 500.00 | 4000.00 | 3988.29 | 3972.05 |
| 6 | 1000.00 | 4000.00 | 3999.71 | 3997.05 |
| 7 | 1500.00 | 4000.00 | 4000.00 | 3999.83 |
| 8 | 2000.00 | 4000.00 | 4000.00 | 4000.00 |
| 9 | 2500.00 | 4000.00 | 4000.00 | 4000.00 |
La presión en el radio del pozo pasa de 3459.13 psi a 1 hora a 3349.49 psi a 12 horas y 3318.90 psi a 24 horas. Redondeamos solo al mostrar la tabla, sin modificar los valores usados en las gráficas.
Graficamos la presión en los diez radios del ejemplo y unimos los puntos para visualizar el perfil a cada tiempo.
plt.figure(figsize=(8, 4.5))
for t in tiempos:
plt.plot(radios, presion(radios, t), "o-", label=f"{t} h")
plt.xlabel("Radio (ft)")
plt.ylabel("Presión (psi)")
plt.grid(alpha=0.3)
plt.legend(title="Tiempo")
plt.tight_layout()
plt.show()
La escala lineal muestra la extensión espacial de la perturbación, pero concentra los puntos cercanos al pozo en el extremo izquierdo.
La misma información puede visualizarse con semilogx. Introducimos los radios directamente: no calculamos su logaritmo antes de graficar.
plt.figure(figsize=(8, 4.5))
for t in tiempos:
plt.semilogx(radios, presion(radios, t), "o-", label=f"{t} h")
plt.xlabel("Radio (ft), escala logarítmica")
plt.ylabel("Presión (psi)")
plt.grid(alpha=0.3)
plt.legend(title="Tiempo")
plt.tight_layout()
plt.show()
Esta escala permite distinguir mejor las variaciones cerca del pozo. Las líneas conectan puntos de evaluación; no representan un cambio brusco real en el perfil entre radios consecutivos.
Que la tabla muestre 4000.00 psi por redondeo no significa que exista una frontera exacta sin perturbación. Esta solución se aproxima gradualmente a la presión inicial.