Del balance de materiales a las curvas de recuperación, GOR y saturación
Al disminuir la presión por debajo del punto de burbuja, parte del gas disuelto se libera. Ese gas modifica las saturaciones y las movilidades, por lo que la relación gas–aceite producida también cambia. El método de Tracy conecta estos cambios con el balance de materiales para estimar cuánto aceite y gas se han producido a cada presión.
Seguiremos el ejemplo 5.4 de Advanced Reservoir Engineering, de Ahmed y McKinney, agregando gráficas para interpretar los resultados. El flujo será: cargar PVT y kr, resolver cada descenso de presión, actualizar los acumulados y graficar.
Este es un pronóstico contra presión, no contra tiempo. Para obtener fechas o caudales harían falta un modelo de productividad, condiciones de operación y una estrategia de explotación. Tampoco debe confundirse con el modelo de invasión de agua de Carter–Tracy.
Consideramos un yacimiento volumétrico con empuje por gas en solución, inicialmente en su presión de burbuja, sin casquete inicial de gas, sin invasión ni producción de agua. Despreciamos la compresibilidad de roca y agua y mantenemos el agua connata inmóvil.
| Dato | Valor |
|---|---|
| Aceite original, N | 15 millones de STB |
| Presión inicial y de burbuja | 4350 psia |
| Presión final del ejercicio | 3350 psia |
| Saturación de agua connata, Swi | 0.30 |
| Viscosidad del aceite, μo | 1.7 cP |
| Viscosidad del gas, μg | 0.023 cP |
| Producción acumulada inicial | 0 |
| GOR inicial | 840 scf/STB |
La tabla PVT procede del ejemplo. Las viscosidades corresponden al cálculo mostrado a 4150 psia; aquí las mantenemos constantes en todo el intervalo por no disponer de una tabla de viscosidades contra presión.
La figura 5.7 contiene las curvas de permeabilidad relativa, pero no una tabla numérica. Usaremos una aproximación tabulada de krg/kro, obtenida por lectura visual de esa figura, con el punto Sg = 0.007, krg/kro = 0.00008 indicado en la solución y un origen de gas inmóvil. Interpolaremos linealmente entre los puntos. Esta entrada es aproximada, especialmente cerca del inicio de movilidad del gas: los resultados son una adaptación reproducible del ejemplo, no una reproducción exacta de su tabla.
Descarga la tabla PVT y la relación krg/kro aproximada, y guarda ambos archivos junto al notebook. Se requieren pandas, NumPy, Matplotlib y SciPy: pip install pandas numpy matplotlib scipy.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import brentq
pvt = pd.read_csv("pvt.csv")
kr = pd.read_csv("kr_aproximada.csv")
N = 15e6
Swi = 0.30
muo, mug = 1.7, 0.023
Boi, Rsi = pvt.loc[0, ["Bo", "Rs"]]
pvt| p | Bo | Bg | Rs | |
|---|---|---|---|---|
| 0 | 4350 | 1.430 | 0.00069 | 840 |
| 1 | 4150 | 1.420 | 0.00071 | 820 |
| 2 | 3950 | 1.395 | 0.00074 | 770 |
| 3 | 3750 | 1.380 | 0.00078 | 730 |
| 4 | 3550 | 1.360 | 0.00081 | 680 |
| 5 | 3350 | 1.345 | 0.00085 | 640 |
Bg está en bbl/scf, mientras Bo está en bbl/STB. Los dos primeros Bg del archivo son 6.9e-4 y 7.1e-4; conservamos la escala coherente con las demás filas y el cálculo del denominador del ejemplo.
fig, ax = plt.subplots(2, 2, figsize=(9, 6))
for a, columna, unidad in zip(ax.flat[:3], ["Bo", "Bg", "Rs"],
["bbl/STB", "bbl/scf", "scf/STB"]):
a.plot(pvt.p, pvt[columna], "o-")
a.set(xlabel="Presión (psia)", ylabel=f"{columna} ({unidad})")
a.invert_xaxis()
ax[1, 1].plot(kr.Sg, kr.krg_kro, "o-")
ax[1, 1].set(xlabel="Saturación de gas, Sg", ylabel="krg/kro")
plt.tight_layout()
plt.show()
El aceite se contrae al liberar gas: Bo y Rs disminuyen. Bg aumenta al caer la presión. Al crecer Sg, aumenta la relación krg/kro y el gas puede adquirir una movilidad importante frente al aceite.
Para trabajar con números pequeños definimos:
\[n=\frac{N_p}{N},\qquad g=\frac{G_p}{N}\]
n es la fracción recuperada; g tiene unidades de scf/STB, no es una fracción de recuperación de gas. El balance de Tracy se expresa como:
\[1=n\Phi_o+g\Phi_g\]
\[D=(B_o-B_{oi})+(R_{si}-R_s)B_g\]
\[\Phi_o=\frac{B_o-R_sB_g}{D},\qquad \Phi_g=\frac{B_g}{D}\]
En la presión inicial, D es cero. Por eso inicializamos el estado a 4350 psia y comenzamos los cálculos a la siguiente presión, sin dividir entre cero.
Para pasar de un estado anterior a uno nuevo aproximamos el GOR medio del intervalo mediante:
\[\overline{R}=\frac{R_{anterior}+R}{2}\]
\[\Delta n=\frac{1-n_{anterior}\Phi_o-g_{anterior}\Phi_g} {\Phi_o+\overline{R}\Phi_g}\]
Luego actualizamos \(n=n_{anterior}+\Delta n\) y \(g=g_{anterior}+\overline{R}\Delta n\).
El GOR de llegada debe ser consistente con las saturaciones y las movilidades:
\[S_o=(1-S_{wi})(1-n)\frac{B_o}{B_{oi}},\qquad S_g=1-S_{wi}-S_o\]
\[R=R_s+\frac{k_{rg}}{k_{ro}}\frac{\mu_o B_o}{\mu_g B_g}\]
Estas ecuaciones están acopladas: R afecta la recuperación y esta afecta Sg y el propio R. En lugar de probar valores manualmente, resolveremos esa igualdad con brentq.
La siguiente función reúne el cálculo repetitivo. En cada presión buscamos el GOR entre Rs y un límite obtenido con la mayor relación krg/kro disponible. Revisamos además que la saturación calculada permanezca dentro de la tabla: no extrapolamos kr.
def pronosticar(tabla):
n, g, R = 0.0, 0.0, Rsi
filas = [[tabla.p.iloc[0], n, g, R, 1-Swi, 0.0, 0.0]]
for v in tabla.iloc[1:].itertuples():
D = (v.Bo-Boi) + (Rsi-v.Rs)*v.Bg
phi_o = (v.Bo-v.Rs*v.Bg)/D
phi_g = v.Bg/D
def estado(Rnuevo):
Rmedio = (R+Rnuevo)/2
dn = (1-n*phi_o-g*phi_g)/(phi_o+Rmedio*phi_g)
nn = n+dn
gg = g+Rmedio*dn
So = (1-Swi)*(1-nn)*v.Bo/Boi
Sg = 1-Swi-So
if not kr.Sg.min() <= Sg <= kr.Sg.max():
raise ValueError("Sg está fuera de la tabla de kr")
relacion = np.interp(Sg, kr.Sg, kr.krg_kro)
Rcalculado = v.Rs+relacion*muo*v.Bo/(mug*v.Bg)
return nn, gg, So, Sg, Rcalculado
limite = v.Rs+kr.krg_kro.max()*muo*v.Bo/(mug*v.Bg)
Rnuevo = brentq(lambda r: estado(r)[4]-r, v.Rs, limite)
n, g, So, Sg, _ = estado(Rnuevo)
R = Rnuevo
error = n*phi_o+g*phi_g-1
filas.append([v.p, n, g, R, So, Sg, error])
return pd.DataFrame(filas, columns=[
"p", "FR", "Gp_N", "GOR", "So", "Sg", "error_MBE"
])
resultado = pronosticar(pvt)
resultado["Np_MMSTB"] = resultado.FR*N/1e6
resultado["Gp_MMscf"] = resultado.Gp_N*N/1e6
resultado[["p", "Np_MMSTB", "Gp_MMscf", "GOR", "So", "Sg"]].round(4)| p | Np_MMSTB | Gp_MMscf | GOR | So | Sg | |
|---|---|---|---|---|---|---|
| 0 | 4350 | 0.0000 | 0.0000 | 840.0000 | 0.7000 | 0.0000 |
| 1 | 4150 | 0.0440 | 36.7921 | 831.7161 | 0.6931 | 0.0069 |
| 2 | 3950 | 0.1675 | 153.7492 | 1062.5052 | 0.6752 | 0.0248 |
| 3 | 3750 | 0.3263 | 349.3291 | 1400.3332 | 0.6608 | 0.0392 |
| 4 | 3550 | 0.4785 | 625.6911 | 2232.3075 | 0.6445 | 0.0555 |
| 5 | 3350 | 0.5953 | 939.0557 | 3134.7170 | 0.6323 | 0.0677 |
A 4150 psia, el cálculo PVT da D = 0.0042, Φo = 199.4762 y Φg = 0.1690476. Conservamos estos valores sin redondearlos antes de resolver. El cálculo avanza en intervalos de 200 psi hasta 3350 psia.
fig, ax = plt.subplots(1, 2, figsize=(9, 4))
ax[0].plot(resultado.p, resultado.Np_MMSTB, "o-")
ax[0].set(ylabel="Aceite acumulado (MMSTB)")
ax[1].plot(resultado.p, resultado.Gp_MMscf, "o-")
ax[1].set(ylabel="Gas acumulado (MMscf)")
for a in ax:
a.set_xlabel("Presión (psia)")
a.invert_xaxis()
plt.tight_layout()
plt.show()
Con pasos de 200 psi, a 3350 psia se obtienen aproximadamente 0.595 MMSTB de aceite, 939.1 MMscf de gas y 3.97 % de recuperación de aceite. El GOR instantáneo es aproximadamente 3135 scf/STB. Estos valores corresponden a la aproximación de kr y viscosidades declarada al inicio.
Ambos acumulados crecen durante el agotamiento. La forma de las curvas permite ver cuánto cambia la producción por cada descenso de presión; su pendiente no es un caudal, porque el eje horizontal no es tiempo.
El GOR instantáneo describe la relación entre las tasas de producción de gas y aceite. Rs representa solamente el gas disuelto por unidad de aceite. Por su parte, \(R_p=G_p/N_p\) es la relación acumulada y contiene la historia de producción.
Rp = resultado.Gp_N / resultado.FR.replace(0, np.nan)
plt.plot(resultado.p, resultado.GOR, "o-", label="GOR instantáneo")
plt.plot(pvt.p, pvt.Rs, "o-", label="Rs")
plt.plot(resultado.p, Rp, "o-", label="Rp = Gp/Np")
plt.xlabel("Presión (psia)")
plt.ylabel("Relación gas–aceite (scf/STB)")
plt.gca().invert_xaxis()
plt.legend()
plt.show()
La separación entre GOR y Rs corresponde al aporte de gas libre móvil. Aunque Rs disminuye, el GOR puede aumentar porque crece la movilidad relativa del gas. Rp suele responder más lentamente, pues promedia toda la producción previa.
fig, ax = plt.subplots(1, 2, figsize=(9, 4))
ax[0].plot(resultado.p, resultado.So, "o-", label="So")
ax[0].plot(resultado.p, resultado.Sg, "o-", label="Sg")
ax[0].axhline(Swi, linestyle="--", label="Sw")
ax[0].set_ylabel("Saturación (fracción)")
ax[0].legend()
ax[1].plot(resultado.p, 100*resultado.FR, "o-")
ax[1].set_ylabel("Recuperación de aceite (%)")
for a in ax:
a.set_xlabel("Presión (psia)")
a.invert_xaxis()
plt.tight_layout()
plt.show()
La saturación de aceite disminuye por la producción y la contracción del aceite remanente. El gas ocupa una fracción creciente del volumen poroso, mientras Sw permanece constante bajo los supuestos de este caso. La recuperación al final del gráfico corresponde a 3350 psia, no a una recuperación final económica del yacimiento.
Cerrar el balance es necesario, pero no demuestra que el paso de presión sea suficientemente pequeño. El GOR medio se aproxima con los valores de los extremos; revisaremos también intervalos de 100 y 50 psi, interpolando linealmente la misma tabla PVT.
assert np.allclose(resultado.So+resultado.Sg+Swi, 1)
assert resultado.error_MBE.abs().max() < 1e-8
base = pvt.sort_values("p")
corridas = {}
for paso in [200, 100, 50]:
presiones = np.arange(4350, 3350-1, -paso)
tabla = pd.DataFrame({"p": presiones})
for columna in ["Bo", "Bg", "Rs"]:
tabla[columna] = np.interp(presiones, base.p, base[columna])
corridas[paso] = pronosticar(tabla)
resumen = pd.DataFrame([
{"Paso_psi": paso, "FR_final_pct": 100*r.FR.iloc[-1],
"GOR_final": r.GOR.iloc[-1],
"Error_MBE_max": r.error_MBE.abs().max()}
for paso, r in corridas.items()
])
resumen| Paso_psi | FR_final_pct | GOR_final | Error_MBE_max | |
|---|---|---|---|---|
| 0 | 200 | 3.968395 | 3134.716988 | 1.110223e-16 |
| 1 | 100 | 3.984968 | 3145.819643 | 1.110223e-16 |
| 2 | 50 | 3.990813 | 3149.735384 | 1.110223e-16 |
for paso, r in corridas.items():
plt.plot(r.p, 100*r.FR, label=f"Paso {paso} psi")
plt.xlabel("Presión (psia)")
plt.ylabel("Recuperación de aceite (%)")
plt.gca().invert_xaxis()
plt.legend()
plt.show()
La recuperación a 3350 psia cambia de 3.9684 % con 200 psi a 3.9908 % con 50 psi: una diferencia de 0.0224 puntos porcentuales. Entre 100 y 50 psi cambia solo 0.0058 puntos. El error del balance queda al nivel del redondeo de máquina.
Esta comparación evalúa la discretización del cálculo, no la incertidumbre del PVT ni de las kr. Reducir el paso no reemplaza una caracterización de fluidos y permeabilidades relativas adecuada.
Pasamos de las propiedades del fluido y la movilidad relativa a un conjunto de curvas de agotamiento. Los gráficos permiten seguir simultáneamente el aceite recuperado, el gas producido y los cambios de saturación, y distinguir el gas disuelto del gas libre que participa en la producción.
El método se aplica aquí con propiedades uniformes y equilibrio de presión representativo del tanque. La segregación gravitacional, heterogeneidad, histéresis, gas atrapado y variaciones de viscosidad pueden modificar el comportamiento real y requieren una descripción adicional.
Nota sobre la fuente. La solución impresa presenta inconsistencias entre algunas fracciones recuperadas, acumulados y conversiones al volumen total. Por ello calculamos los acumulados con las ecuaciones sin redondeo intermedio y verificamos el cierre del balance. Las diferencias respecto a la tabla también dependen de la lectura aproximada de kr y de la hipótesis de viscosidades constantes.
Ahmed, T. y McKinney, P. D. (2005). Advanced Reservoir Engineering. Gulf Professional Publishing. Capítulo 4, sección 4.5: Tracy’s Form of the MBE; capítulo 5, Tracy method, páginas impresas 335–337, ecuaciones 5.1.34–5.1.44, ejemplo 5.4 y figura 5.7. La figura de kr se atribuye en el libro a Economides et al. (1994).