Correlaciones de Aceite en Python: Flujo de Cálculo de Rs, Bo y Viscosidad

De los datos de entrada a las curvas PVT, por debajo y por encima de la presión de burbuja

Yacimientos
PVT
Python

Calcula la presión de burbuja, el gas en solución, el factor de volumen y la viscosidad del aceite con un ejemplo en Python. Flujo paso a paso con Standing, Beggs–Robinson y Vázquez–Beggs.

Autor/a
Fecha de publicación

27 de septiembre de 2026

Un mismo conjunto de datos, tres propiedades

Cuando calculamos propiedades PVT del aceite, las ecuaciones no trabajan de manera aislada. Primero necesitamos conocer la presión de burbuja, Pb, para decidir cómo calcular el gas en solución, el factor de volumen y la viscosidad a cada presión.

En este ejemplo usamos los siguientes métodos: Standing para Pb, Rs y Bo, Beggs–Robinson para viscosidad de aceite muerto y saturado, y Vázquez–Beggs para las correcciones por encima de Pb. No aplicamos corrección de la gravedad del gas por condiciones de separador.

El objetivo es seguir el cálculo completo y generar las tres curvas con fragmentos sencillos de Python.

Propiedad Qué representa Unidades
Rs Gas disuelto por barril de aceite a condiciones de tanque scf/STB
Bo Volumen de aceite con su gas disuelto a condiciones de yacimiento por barril de aceite de tanque rb/STB
μo Viscosidad dinámica del aceite cP
Visita el canal @rigopetrodata

Tutoriales de R, Python y Excel para ingeniería petrolera

Ver canal

Datos del ejemplo

Entrada Valor Unidad
Gravedad del aceite 35 °API
Gravedad relativa del gas 0.75 Aire = 1
Presión de evaluación 4000 psia
Temperatura 180 °F
Gas en solución a Pb, Rsb 650 scf/STB

La temperatura se mantiene constante. Para visualizar ambos regímenes calcularemos de 500 a 5000 psia e incluiremos explícitamente Pb y la presión de evaluación. Este intervalo es una elección para el ejemplo, no un rango universal de validez de las correlaciones.

Rsb no es el GOR producido

Rsb es el gas disuelto en el aceite a la presión de burbuja. No debe sustituirse automáticamente por la relación gas–aceite de producción, que puede incluir gas libre. Las ecuaciones siguientes requieren presión absoluta en psia y temperatura en °F.

El flujo de cálculo

  1. Con API, gravedad del gas, temperatura y Rsb, calcular Pb.
  2. Calcular las propiedades de referencia: Bob, viscosidad de aceite muerto μod y viscosidad a Pb μob.
  3. Para P ≤ Pb, calcular Rs con Standing; usar ese Rs en Bo y μo.
  4. Para P > Pb, mantener Rs = Rsb; corregir Bob y μob por presión.
  5. Reunir los resultados en una tabla y graficar Rs(P), Bo(P) y μo(P).

Paso 1 — Entradas y presión de burbuja

Necesitamos NumPy, pandas y Matplotlib. Puedes instalarlos con python -m pip install numpy pandas matplotlib. Ejecuta los fragmentos en orden, dentro del mismo notebook o sesión.

La gravedad relativa del aceite se obtiene de \(\gamma_o=141.5/(API+131.5)\). Para Pb usamos Standing:

\[ P_b=18.2\left[\left(\frac{R_{sb}}{\gamma_g}\right)^{0.83} 10^{0.00091T-0.0125API}-1.4\right] \]

Ver código
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

API = 35
gamma_g = 0.75
T = 180          # °F
P_eval = 4000    # psia
Rsb = 650       # scf/STB

gamma_o = 141.5 / (API + 131.5)
Pb = 18.2 * ((Rsb / gamma_g)**0.83
             * 10**(0.00091*T - 0.0125*API) - 1.4)

P = np.unique(np.append(np.linspace(500, 5000, 150), [Pb, P_eval]))
print(f"Presión de burbuja: {Pb:.2f} psia")
print("A 4000 psia:", "aceite subsaturado" if P_eval > Pb else "aceite saturado")
Presión de burbuja: 2633.95 psia
A 4000 psia: aceite subsaturado

Pb es el punto de cambio entre las dos ramas del cálculo. A presiones mayores que Pb, el gas del sistema permanece disuelto en el aceite; al reducir la presión por debajo de Pb comienza a liberarse gas.

Paso 2 — Gas en solución, Rs

Invertimos la ecuación de Standing para calcular Rs por debajo de Pb. Por encima de Pb, fijamos Rs en el valor de entrada Rsb:

\[ R_s(P)=\begin{cases} \gamma_g\left[\left(\frac{P}{18.2}+1.4\right) 10^{0.0125API-0.00091T}\right]^{1/0.83}, & P<P_b\\ R_{sb}, & P\geq P_b \end{cases} \]

Usamos \(1/0.83\) como inverso del exponente de Pb para que ambas expresiones coincidan en el punto de burbuja; suele encontrarse redondeado a 1.2048.

Ver código
Rs_sat = gamma_g * ((P / 18.2 + 1.4)
                    * 10**(0.0125*API - 0.00091*T))**(1 / 0.83)
Rs = np.where(P < Pb, Rs_sat, Rsb)

np.where selecciona la expresión correspondiente a cada presión. No continuamos aumentando Rs después de Pb: para este fluido el límite es 650 scf/STB.

Paso 3 — Factor de volumen del aceite, Bo

Para el aceite saturado, Standing relaciona Bo con Rs, las gravedades relativas y la temperatura:

\[ B_{o,sat}=0.9759+0.00012\left[R_s\sqrt{\frac{\gamma_g}{\gamma_o}}+1.25T\right]^{1.2} \]

Al sustituir Rs por Rsb obtenemos Bob, el valor a Pb. Por encima de Pb necesitamos representar la compresión del aceite. Usamos la compresibilidad de Vázquez–Beggs, sin corrección de separador:

\[ c_o(P)=\frac{C}{P},\qquad C=\frac{-1433+5R_{sb}+17.2T-1180\gamma_g+12.61API}{10^5} \]

Como \(c_o=-d\ln B_o/dP\), integrar esta expresión desde Pb hasta P da:

\[B_o(P)=B_{ob}\left(\frac{P}{P_b}\right)^{-C},\qquad P>P_b\]

Esta forma conserva la dependencia de \(c_o\) con la presión. La expresión \(B_{ob}\exp[-c_o(P-P_b)]\) corresponde a tratar \(c_o\) como constante en el intervalo; no son exactamente la misma aproximación.

Ver código
Bo_sat = 0.9759 + 0.00012 * (Rs * np.sqrt(gamma_g / gamma_o) + 1.25*T)**1.2
Bob = 0.9759 + 0.00012 * (Rsb * np.sqrt(gamma_g / gamma_o) + 1.25*T)**1.2

C = (-1433 + 5*Rsb + 17.2*T - 1180*gamma_g + 12.61*API) / 100000
Bo = np.where(P <= Pb, Bo_sat, Bob * (P / Pb)**(-C))

print(f"Bo a Pb: {Bob:.4f} rb/STB")
Bo a Pb: 1.3610 rb/STB

Paso 4 — Viscosidad del aceite, μo

El cálculo tiene tres niveles. Primero estimamos la viscosidad del aceite muerto, sin gas disuelto, a la temperatura del ejemplo. Después incorporamos el efecto de Rs para el aceite saturado. Finalmente corregimos por presión en la región subsaturada.

Con Beggs–Robinson:

\[ X=10^{3.0324-0.02023API}\,T^{-1.163},\qquad \mu_{od}=10^X-1 \]

\[ \mu_{o,sat}=A\mu_{od}^{B},\qquad A=10.715(R_s+100)^{-0.515},\qquad B=5.44(R_s+150)^{-0.338} \]

Evaluando esta última expresión con Rsb obtenemos μob. Para P mayor que Pb usamos Vázquez–Beggs:

\[ \mu_o=\mu_{ob}\left(\frac{P}{P_b}\right)^m,\qquad m=2.6P^{1.187}\exp(-11.513-8.98\times10^{-5}P) \]

Ver código
X = 10**(3.0324 - 0.02023*API) * T**(-1.163)
mu_dead = 10**X - 1

A = 10.715 * (Rs + 100)**(-0.515)
B = 5.44 * (Rs + 150)**(-0.338)
mu_sat = A * mu_dead**B
mu_b = 10.715 * (Rsb + 100)**(-0.515) * mu_dead**(5.44 * (Rsb + 150)**(-0.338))

m = 2.6 * P**1.187 * np.exp(-11.513 - 8.98e-5*P)
mu_o = np.where(P <= Pb, mu_sat, mu_b * (P / Pb)**m)

print(f"Viscosidad de aceite muerto: {mu_dead:.4f} cP")
print(f"Viscosidad a Pb: {mu_b:.4f} cP")
Viscosidad de aceite muerto: 2.1833 cP
Viscosidad a Pb: 0.5520 cP

μod es una referencia de cálculo, no la viscosidad del aceite vivo a la presión de yacimiento. El gas disuelto reduce la viscosidad respecto al aceite muerto en este ejemplo.

🛒
Propiedades De Los Fluidos Del Yacimiento

Disponible en Mercado Libre

Ver oferta →

Paso 5 — Consultar las propiedades a 4000 psia

Reunimos las propiedades calculadas en una tabla. Como incluimos P_eval y Pb explícitamente en el arreglo, podemos consultar esos valores sin interpolar.

Ver código
pvt = pd.DataFrame({
    "P (psia)": P,
    "Rs (scf/STB)": Rs,
    "Bo (rb/STB)": Bo,
    "Viscosidad (cP)": mu_o
})

print(pvt.loc[np.isin(P, [Pb, P_eval])].round(4).to_string(index=False))
 P (psia)  Rs (scf/STB)  Bo (rb/STB)  Viscosidad (cP)
2633.9504         650.0       1.3610           0.5520
4000.0000         650.0       1.3358           0.6369

La fila de Pb muestra las condiciones de referencia. A 4000 psia, Rs permanece en 650 scf/STB, Bo es menor que Bob y la viscosidad es mayor que μob.

Paso 6 — Graficar Rs, Bo y viscosidad

Usamos tres paneles porque las propiedades tienen unidades y escalas distintas. La línea punteada marca Pb y el punto negro identifica la presión de evaluación.

Ver código
fig, axes = plt.subplots(3, 1, figsize=(9, 10), sharex=True)
propiedades = [(Rs, "Rs (scf/STB)"), (Bo, "Bo (rb/STB)"),
               (mu_o, "Viscosidad (cP)")]

for ax, (valores, nombre) in zip(axes, propiedades):
    ax.plot(P, valores, color="steelblue", linewidth=2)
    ax.axvline(Pb, color="gray", linestyle="--", label=f"Pb = {Pb:.0f} psia")
    ax.scatter(P_eval, valores[P == P_eval][0], color="black", label="4000 psia")
    ax.set_ylabel(nombre)
    ax.grid(alpha=0.25)

axes[0].legend()
axes[0].set_title("Correlaciones PVT de aceite: 35 °API y 180 °F")
axes[-1].set_xlabel("Presión (psia)")
plt.tight_layout()
plt.show()
Rs aumenta hasta Pb y luego permanece constante. Bo alcanza su máximo y la viscosidad su mínimo en Pb.
Figura 1: Propiedades del aceite a 180 °F. Línea punteada: Pb; punto negro: 4000 psia.

¿Cómo leer las curvas?

Si recorremos el gráfico desde baja hacia alta presión:

  • Rs aumenta hasta Pb y después queda constante. El aceite no puede incorporar más gas que el Rsb definido para este sistema.
  • Bo aumenta hasta Pb por el efecto del gas disuelto. Por encima de Pb disminuye debido a la compresión del aceite.
  • La viscosidad disminuye hasta Pb al aumentar el gas en solución. Por encima de Pb aumenta con la presión, manteniendo Rs constante.

En una depleción partiendo de 4000 psia, el recorrido es en sentido contrario: primero el aceite se expande al reducir la presión hasta Pb; después libera gas, se reduce Bo y aumenta la viscosidad del aceite remanente.

Alcance del ejemplo

Estas curvas son estimaciones por correlaciones, no datos de laboratorio. Se mantienen constantes API, gravedad del gas, temperatura y Rsb; no se modelan cambios composicionales, emulsiones ni una separación multietapa. La selección de métodos sigue la pantalla de PVT oil; para Bo subsaturado se ha explicitado la integración de la compresibilidad, por lo que una app que use compresibilidad constante puede dar un valor ligeramente distinto. Antes de aplicar el resultado a otro crudo, revisa el dominio de cada correlación y contrasta con datos PVT medidos.

Para repetir el ejercicio

Modifica únicamente las entradas del primer fragmento y vuelve a ejecutar todas las celdas. Comprueba que el intervalo de presión cubra las condiciones que deseas estudiar y revisa la unión de las curvas en Pb. No se necesita un dataset externo: todos los datos de entrada están al inicio del ejemplo.

Referencias