Cubic Equations of State

There is one more thing to show of the Python API that yaeos supports. At this point of the tutorial we have been using cubic equations of state to perform different calculations. The experience users may have noticed that we never set any interaction parameters for the cubic equations of state.

yaeos support the following mixing rules for cubic equations of state:

• QMR: quadratic mixing rule

• MHV: modified Huron-Vidal mixing rule

You can check the available mixing rules by running the following command:

import yaeos

?yaeos.models

Type:        module
String form: <module 'yaeos.models' from '/home/salvador/code/yaeos/python/yaeos/models/__init__.py'>
File:        ~/code/yaeos/python/yaeos/models/__init__.py
Docstring:
Models module.

Yaeos models module. This module provides the following submodules:

- excess_gibbs: Excess Gibbs energy models
    - NRTL: non-random two-liquid model
    - UNIFACVLE: Original UNIFAC VLE model
    - UNIQUAC: UNIversal QUAsiChemical Excess Gibbs free energy model

- residual_helmholtz: Residual Helmholtz energy models
    - Cubic EoS:
        - PengRobinson76: Peng-Robinson model (1976)
        - PengRobinson78: Peng-Robinson model (1978)
        - SoaveRedlichKwong: Soave-Redlich-Kwong model
        - RKPR: RKPR model
    - Mixing rules: mixing rules for cubic EoS
        - QMR: cuadratic mixing rule
        - MHV: modified Huron-Vidal mixing rule

Quadratic mixing rule (QMR)

The QMR is the default mixing rule for all cubic equations of state. But it is always set with all the interaction parameters equals to zero by default.

To set the \(k_{ij}\) and \(l_{ij}\) interaction matrix you can do the following:

from yaeos import PengRobinson76, QMR
import numpy as np

tcs = np.array([190.6, 369.8])  # critical temperatures [K]
pcs = np.array([45.99, 42.48])  # critical pressures [bar]
w = np.array([0.0115, 0.0866])  # acentric factors [-]

# binary interaction parameters
kij = np.array([[0.0, 0.1], [0.1, 0.0]])

lij = np.array([[0.0, 0.01], [0.01, 0.0]])

# Set QMR
mixing_rule = QMR(kij=kij, lij=lij)

# Set a model with the QMR with parameters
model_kij = PengRobinson76(
    critical_temperatures=tcs,
    critical_pressures=pcs,
    acentric_factors=w,
    mixrule=mixing_rule,
)

# Set model without QMR (QMR with all kij and lij as zero)
model = PengRobinson76(tcs, pcs, w)

There we have instantiated a PengRobinson76 with the QMR mixing rule with fantasy interaction parameters. And then we instantiated a PengRobinson76 with no interaction parameters as we have been doing.

With the model with interaction parameters we can calculate the same things that we have been calculating at this tutorial. But to show that there is a difference, we can calculate a phase envelope and see what happens.

import matplotlib.pyplot as plt

z = np.array([0.5, 0.5])
t0 = 60.0
p0 = 1.0

# Envelope with the model with kij and lij
env_kij = model_kij.phase_envelope_pt(z, t0=t0, p0=p0)

# Envelope with the model without kij and lij
env = model.phase_envelope_pt(z, t0=t0, p0=p0)

# Plot
plt.plot(env_kij["Ts"], env_kij["Ps"], label="With kij and lij")
plt.plot(env["Ts"], env["Ps"], label="Without kij and lij")

plt.legend()
plt.xlabel("Temperature [K]")
plt.ylabel("Pressure [bar]")

Text(0, 0.5, 'Pressure [bar]')

The QMR mix rule can be used with any cubic equation of state that yaeos supports.

Huron-Vidal mixing rule

The Huron-Vidal mixing rule allows linking a Cubic Equation of State with an Excess Gibbs model, this can be helpful to predict behaviour of systems with high asymmetry in polarity.

from yaeos import NRTL, HV, PengRobinson76
import numpy as np

tcs = np.array([190.6, 369.8])  # critical temperatures [K]
pcs = np.array([45.99, 42.48])  # critical pressures [bar]
w = np.array([0.0115, 0.0866])  # acentric factors [-]

# binary interaction parameters for NRTL model
aij = np.array([[0.0, -1.2], [0.5, 0.0]])

bij = np.array([[0.0, 100.0], [-100.0, 0.0]])

cij = np.array([[0.0, 0.3], [0.3, 0.0]])

nrtl = NRTL(aij, bij, cij)

mixing_rule = HV(ge=nrtl)

# Set a model with the MHV with parameters
model_hv_nrtl = PengRobinson76(
    critical_temperatures=tcs,
    critical_pressures=pcs,
    acentric_factors=w,
    mixrule=mixing_rule,
)

Modified Huron-Vidal mixing rule (MHV)

The other mixing rule that yaeos supports is the MHV mixing rule. To use it we need to instantiate an Excess Gibbs free energy model and pass it to the MHV mixing rule. The only Excess Gibbs free energy model that yaeos supports is the NRTL model. Let’s see an example:

from yaeos import NRTL, MHV, PengRobinson76
import numpy as np

mixing_rule = MHV(ge=nrtl, q=-0.53)

# Set a model with the MHV with parameters
model_mhv_nrtl = PengRobinson76(
    critical_temperatures=tcs,
    critical_pressures=pcs,
    acentric_factors=w,
    mixrule=mixing_rule,
)

# Set a model with null kij and lij to compare later
model = PengRobinson76(
    critical_temperatures=tcs, critical_pressures=pcs, acentric_factors=w
)

Let’s make first an envelope to see the difference and then explain what is happening.

import matplotlib.pyplot as plt

z = np.array([0.2, 0.8])
t0 = 100.0
p0 = 1.0

# Envelope with the model with HV and NRTL
env_mhv_nrtl = model_mhv_nrtl.phase_envelope_pt(
    z, t0=t0, p0=p0, max_points=5000
)

# Envelope with the model with MHV and NRTL
env_hv_nrtl = model_hv_nrtl.phase_envelope_pt(z, t0=t0, p0=p0, max_points=5000)

# Envelope with the model without kij and lij
env = model.phase_envelope_pt(z, t0=t0, p0=p0, max_points=5000)

# Plot
plt.plot(env_hv_nrtl["Ts"], env_hv_nrtl["Ps"], label="HV-NRTL")
plt.plot(env_mhv_nrtl["Ts"], env_mhv_nrtl["Ps"], label="MHV-NRTL")
plt.plot(env["Ts"], env["Ps"], label="Without kij and lij")

plt.legend()
plt.xlabel("Temperature [K]")
plt.ylabel("Pressure [bar]")

 WARN: possible bad root solving:                   Infinity   7.2869289762081603E+026
 WARN: possible bad root solving:                   Infinity   7.2869289762081603E+026

Text(0, 0.5, 'Pressure [bar]')

We have instantiated a NRTL model with some interaction parameters and then instantiated a MHV mixing rule with the NRTL model. The q value of the MHV mixing rule is taken from bibliography. The value to use are:

• q = -0.594 for Soave-Redlich-Kwong

• q = -0.53 for Peng-Robinson

• q = -0.85 for Van der Waals