← Volver al blog

Compressibility Factor Calculation (Peng-Robinson)

03 de junio de 2023 2 min de lectura cubic equation of state peng-robinson compressibility factor code python thermodynamics

Calcular el factor de compresibilidad \(Z\) de un gas real mediante la ecuación cúbica de estado de Peng-Robinson [1]. Las propiedades críticas (\(T_c\), \(P_c\), \(\omega\)) se extraen automáticamente de CoolProp — no hace falta hardcodearlas por compuesto.

La EoS de Peng-Robinson se escribe en forma reducida como:

\[Z = 1 + \beta - q\,\beta \cdot \frac{Z - \beta}{(Z + \varepsilon\,\beta)(Z + \sigma\,\beta)}\]

con:

\[\beta = \frac{\Omega \, P_r}{T_r}, \qquad q = \frac{\Psi \, \alpha(T_r, \omega)}{\Omega \, T_r}\]
\[\alpha(T_r, \omega) = \bigl[1 + \kappa(\omega)(1 - \sqrt{T_r})\bigr]^2, \qquad \kappa(\omega) = 0.37464 + 1.5422\,\omega - 0.26992\,\omega^2\]

donde \(\Omega = 0.07780\), \(\Psi = 0.45724\), \(\sigma = 1 + \sqrt{2}\) y \(\varepsilon = 1 - \sqrt{2}\).

Como la ecuación es implícita en \(Z\), se resuelve por iteración de punto fijo partiendo de \(Z_0 \approx 1\) (gas ideal) y recalculando hasta que \(|Z_{n+1} - Z_n| < \text{tol}\).

import numpy as np
import CoolProp.CoolProp as CP

class Zfactor:
    def __init__(self, component_name, T, P):
        self.component_name = component_name
        self.T = T
        self.P = P

    def calculation_z(self, tol = 1E-4, i_max = 100):
        Tc = CP.PropsSI('Tcrit', self.component_name)
        Pc = CP.PropsSI('Pcrit', self.component_name)
        w  = CP.PropsSI('acentric', self.component_name)

        #Parameters of the Peng-Robinson Equation
        Ω = 0.07780
        Ψ = 0.45724
        σ = 1 + np.sqrt(2)
        ε = 1 - np.sqrt(2)

        #==============Functions==============
        Tr = lambda T: T/Tc
        Pr = lambda P: P/Pc
        acent_w = lambda: 0.37464 + 1.5422*w - 0.26992*w**2
        α = lambda T: (1 + (acent_w())*(1 - np.sqrt(Tr(T))))**2
        β = lambda T, P: (Ω*Pr(P))/Tr(T)
        q = lambda T: (Ψ*α(T))/(Ω*Tr(T))

        #==============Algorithm for calculating compressibility factor==============
        z_previous = 0.9999
        i = 0
        while 1:
            z = 1 + β(self.T, self.P) - q(self.T)*β(self.T, self.P)*((z_previous - β(self.T, self.P))/((z_previous + ε*β(self.T, self.P))*(z_previous + σ*β(self.T, self.P))))
            if abs(z - z_previous) < tol or i >= i_max:
                break
            z_previous = z
            i = i + 1
        return round(z, 4)

Ejemplo con metano a \(T = 300\,\text{K}\) y \(P = 5\,\text{MPa}\):

z = Zfactor("Methane", 300, 5e6).calculation_z()
print(z)   # -> 0.9018

Barrido en presión para distintas temperaturas:

Methane — Z vs P curves

Cómo leer el gráfico:

  • A \(T\) alta (\(400\,\text{K}\), verde) el metano se comporta casi como gas ideal: \(Z\) se mantiene cerca de \(1\) incluso a \(30\,\text{MPa}\).
  • A \(T\) baja (\(200\,\text{K}\), roja, por debajo de \(T_c \approx 190.6\,\text{K}\)) las moléculas se atraen: \(Z\) cae muy por debajo de \(1\) antes de recuperarse a alta presión por el volumen excluido.
  • La curva cruza \(Z = 1\) cerca de la temperatura de Boyle, donde las fuerzas atractivas y repulsivas se cancelan exactamente.

Por qué Peng-Robinson y no otra EoS:

  • Precisión mucho mejor que Van der Waals o Redlich-Kwong para hidrocarburos y mezclas no-polares.
  • La corrección de Soave en \(\alpha(T_r, \omega)\) incorpora la acentricidad \(\omega\), capturando bien la dependencia térmica cerca del punto crítico.
  • Es la EoS estándar en la industria del petróleo y gas — la mayoría de simuladores de procesos (Aspen HYSYS, Pro/II, DWSIM) la usan por defecto.

Referencia

[1] J. M. Smith, H. C. Van Ness, M. M. Abbott, M. T. Swihart. Introduction to Chemical Engineering Thermodynamics (8th ed.). McGraw-Hill, 2018.

Instalación:

pip install coolprop numpy
← Volver al blog