Implementing an equation of state in python

In FeOs, you can implement your equation of state in python, register it to the Rust backend, and compute properties and phase equilbria as if you implemented it in Rust. In this tutorial, we will implement the Peng-Robinson equation of state.

Table of contents

• Implementation

• Computing properties

• Critical point

• Phase equilibria and phase diagrams

• Mixtures

• Comparison to Rust implementation

from feos.eos import *
from feos.si import *
import numpy as np

optional = True

import warnings
warnings.filterwarnings('ignore')

if optional:
    import matplotlib.pyplot as plt
    import pandas as pd
    import seaborn as sns

    sns.set_style("ticks")
    sns.set_palette("Dark2")
    sns.set_context("talk")

Implementation

↑ Back to top

To implement an equation of state in python, we have to define a class which has to have the following methods:

class MyEquationOfState:
    def helmholtz_energy(self, state: StateHD) -> D

    def components(self) -> int

    def subset(self, indices: List[int]) -> Self

    def molar_weight(self) -> SIArray1

    def max_density(self, moles: SIArray1) -> f64

• components(self) -> int: Returns the number of components (usually inferred from the shape of the input parameters).

• molar_weight(self) -> SIArray1: Returns an SIArray1 with size equal to the number of components containing the molar mass of each component.

• max_density(self, moles: np.ndarray[float]) -> float: Returns the maximum allowed number density in units of molecules/Angstrom³.

• subset(self, indices: List[int]) -> self: Returns a new equation of state with parameters defined in indices.

• helmholtz_energy(self, state: StateHD) -> Dual: Returns the helmholtz energy as (hyper)-dual number given a StateHD.

SQRT2 = np.sqrt(2)

class PyPengRobinson:
    def __init__(
        self, critical_temperature, critical_pressure,
        acentric_factor, molar_weight, delta_ij=None
    ):
        """Peng-Robinson Equation of State

        Parameters
        ----------
        critical_temperature : SIArray1
            critical temperature of each component.
        critical_pressure : SIArray1
            critical pressure of each component.
        acentric_factor : np.array[float]
            acentric factor of each component (dimensionless).
        molar_weight: SIArray1
            molar weight of each component.
        delta_ij : np.array[[float]], optional
            binary parameters. Shape=[n, n], n = number of components.
            defaults to zero for all binary interactions.

        Raises
        ------
        ValueError: if the input values have incompatible sizes.
        """
        self.n = len(critical_temperature)
        if len(set((
            len(critical_temperature),
            len(critical_pressure),
            len(acentric_factor)
        ))) != 1:
            raise ValueError("Input parameters must all have the same lenght.")

        self.tc = critical_temperature / KELVIN
        self.pc = critical_pressure / PASCAL
        self.omega = acentric_factor
        self.mw = molar_weight / GRAM * MOL

        self.a_r = (0.45724 * critical_temperature**2 * RGAS
                    / critical_pressure / ANGSTROM**3 / NAV / KELVIN)
        self.b = (0.07780 * critical_temperature * RGAS
                  / critical_pressure / ANGSTROM**3 / NAV)
        self.kappa = (0.37464
                      + (1.54226 - 0.26992 * acentric_factor) * acentric_factor)
        self.delta_ij = (np.zeros((self.n, self.n))
                         if delta_ij is None else delta_ij)

    def helmholtz_energy(self, state):
        """Return helmholtz energy.

        Parameters
        ----------
        state : StateHD
            The thermodynamic state.

        Returns
        -------
        helmholtz_energy: float | any dual number
            The return type depends on the input types.
        """
        n = np.sum(state.moles)
        x = state.molefracs
        tr = 1.0 / self.tc * state.temperature
        ak = ((1.0 - np.sqrt(tr)) * self.kappa + 1.0)**2 * self.a_r
        ak_mix = 0.0
        if self.n > 1:
            for i in range(self.n):
                for j in range(self.n):
                    ak_mix += (np.sqrt(ak[i] * ak[j])
                               * (x[i] * x[j] * (1.0 - self.delta_ij[i, j])))
        else:
            ak_mix = ak[0]
        b = np.sum(x * self.b)
        v = state.volume
        rho = np.sum(state.partial_density)
        a = n * (-np.log(1.0 - b * rho)
                 - ak_mix / (b * SQRT2 * 2.0 * state.temperature)
                 * np.log((1.0 + (1.0 + SQRT2) * b * rho) / (1.0 + (1.0 - SQRT2) * b * rho)))
        return a

    def components(self) -> int:
        """Number of components."""
        return self.n

    def subset(self, i: [int]):
        """Return new equation of state containing a subset of all components."""
        if self.n > 1:
            tc = self.tc[i]
            pc = self.pc[i]
            mw = self.mw[i]
            omega = self.omega[i]
            return PyPengRobinson(
                SIArray1(tc * KELVIN),
                SIArray1(pc * PASCAL),
                omega,
                SIArray1(mw * GRAM / MOL)
            )
        else:
            return self

    def molar_weight(self) -> SIArray1:
        return SIArray1(self.mw * GRAM / MOL)

    def max_density(self, moles:[float]) -> float:
        b = np.sum(moles * self.b) / np.sum(moles);
        return 0.9 / b

Computing properties

↑ Back to top

Let’s compute some properties. First, we have to instanciate the class and register it to Rust. This is done using the EquationOfState.python method.

tc = SIArray1(369.96 * KELVIN)
pc = SIArray1(4250000.0 * PASCAL)
omega = np.array([0.153])
molar_weight = SIArray1(44.0962 * GRAM / MOL)

# create an instance of our python class and hand it over to rust
pr = PyPengRobinson(tc, pc, omega, molar_weight)
eos = EquationOfState.python(pr)

Thermodynamic state: the State object

Before we can compute a property, we create a State object. This can be done in several ways depending on what control variables we need. If no total amount of substance is defined, it is set to \(n = \frac{1}{N_{AV}}\). For possible input combinations, you can inspect the signature of the constructor using State?.

State?

Init signature: State(self, /, *args, **kwargs)
Docstring:
A thermodynamic state at given conditions.

Parameters
----------
eos : Eos
    The equation of state to use.
temperature : SINumber, optional
    Temperature.
volume : SINumber, optional
    Volume.
density : SINumber, optional
    Molar density.
partial_density : SIArray1, optional
    Partial molar densities.
total_moles : SINumber, optional
    Total amount of substance (of a mixture).
moles : SIArray1, optional
    Amount of substance for each component.
molefracs : numpy.ndarray[float]
    Molar fraction of each component.
pressure : SINumber, optional
    Pressure.
molar_enthalpy : SINumber, optional
    Molar enthalpy.
molar_entropy : SINumber, optional
    Molar entropy.
molar_internal_energy: SINumber, optional
    Molar internal energy
density_initialization : {'vapor', 'liquid', SINumber, None}, optional
    Method used to initialize density for density iteration.
    'vapor' and 'liquid' are inferred from the maximum density of the equation of state.
    If no density or keyword is provided, the vapor and liquid phase is tested and, if
    different, the result with the lower free energy is returned.
initial_temperature : SINumber, optional
    Initial temperature for temperature iteration. Can improve convergence
    when the state is specified with pressure and molar entropy or enthalpy.

Returns
-------
State : state at given conditions

Raises
------
Error
    When the state cannot be created using the combination of input.
Type:           type
Subclasses:

If we use input variables other than \(\mathbf{N}, V, T\) (the natural variables of the Helmholtz energy), creating a state is an iterative procedure. For example, we can create a state for a give \(T, p\), which will result in a iteration of the volume (density).

# If no amount of substance is given, it is set to 1/NAV.
s = State(eos, temperature=300*KELVIN, pressure=1*BAR)
s.total_moles

$1.6605\times10^{-24}\,\mathrm{ mol}$

s_pt = State(
    eos,
    temperature=300*KELVIN,
    pressure=1*BAR,
    total_moles=1*MOL
)
s_pt.total_moles

$1\,\mathrm{ mol}$

We can use other variables as well. For example, we can create a state at given \(h, p\) (using the enthalpy from the prior computation as input):

s_ph = State(
    eos,
    pressure=1*BAR,
    molar_enthalpy=s_pt.molar_enthalpy()
)

# check if states are equal
print("rel. dev.")
print(
    "entropy    : ",
    (s_ph.molar_entropy() - s_pt.molar_entropy()) / s_pt.molar_entropy()
)
print(
    "density    : ",
    (s_ph.mass_density() - s_pt.mass_density()) / s_pt.mass_density()
)
print(
    "temperature: ",
    (s_ph.temperature - s_pt.temperature) / s_pt.temperature
)

rel. dev.
entropy    :  1.2028878239185388e-16
density    :  -1.8532974044002563e-15
temperature:  1.5158245029548805e-15

Critical point

↑ Back to top

To generate a state at critical conditions, we can use the critical_point constructor.

s_cp = State.critical_point(eos)
print("Critical point")
print("temperature: ", s_cp.temperature)
print("density    : ", s_cp.mass_density())
print("pressure   : ", s_cp.pressure())

Critical point
temperature:  369.9506174234607 K
density    :  198.1862458057177 kg/m³
pressure   :  4.249677749116942 MPa

Phase equilibria and phase diagrams

↑ Back to top

We can also create an object, PhaseEquilibrium, that contains states that are in equilibrium.

vle = PhaseEquilibrium.pure(eos, 300.0*KELVIN)
vle

	temperature	density
phase 1	300.00000 K	488.99014 mol/m³
phase 2	300.00000 K	11.53399 kmol/m³

Each phase is a State object. We can simply access these states and compute properties, just like before.

vle.liquid # the high density phase `State`

temperature	density
300.00000 K	11.53399 kmol/m³

vle.vapor # the low density phase `State`

temperature	density
300.00000 K	488.99014 mol/m³

# we can now easily compute any property:
print("Heat of vaporization: ", vle.vapor.molar_enthalpy() - vle.liquid.molar_enthalpy())
print("for T = {}".format(vle.liquid.temperature))
print("and p = {:.2f} bar".format(vle.liquid.pressure() / BAR))

Heat of vaporization:  14.782343503305114 kJ/mol
for T = 300 K
and p = 9.95 bar

We can also easily compute vapor pressures and boiling temperatures:

# This also works for mixtures, in which case the pure component properties are computed.
# Hence, the result is a list - that is why we use an index [0] here.
print("vapor pressure      (T = 300 K):", PhaseEquilibrium.vapor_pressure(eos, 300*KELVIN)[0])
print("boiling temperature (p = 3 bar):", PhaseEquilibrium.boiling_temperature(eos, 2*BAR)[0])

vapor pressure      (T = 300 K): 994.7761635610093 kPa
boiling temperature (p = 3 bar): 247.84035574956746 K

Phase Diagram

We could repeatedly compute PhaseEquilibrium states for different temperatures / pressures to generate a phase diagram. Because this a common task, there is a object for that as well.

The PhaseDiagram object creates multiple PhaseEquilibrium objects (npoints of those) between a given lower temperature and the critical point.

dia = PhaseDiagram.pure(eos, 230.0 * KELVIN, 500)

We can access each PhaseEquilbrium and then conveniently comput any property we like:

enthalpy_of_vaporization = [
    (vle.vapor.molar_enthalpy() - vle.liquid.molar_enthalpy()) / (KILO * JOULE) * MOL
    for vle in dia.states
]

fig, ax = plt.subplots(figsize=(7, 4))
sns.lineplot(x=dia.vapor.temperature / KELVIN, y=enthalpy_of_vaporization, ax=ax);
ax.set_ylabel(r"$\Delta^{LV}h$ / kJ / mol")
ax.set_xlabel(r"$T$ / K");

A more convenient way is to create a dictionary. The dictionary can be used to build pandas DataFrame objects. This is a bit less flexible, because the units of the properties are rigid. You can inspect the method signature to check what units are used.

dia.to_dict?

Signature: dia.to_dict()
Docstring:
Returns the phase diagram as dictionary.

Units
-----
temperature : K
pressure : Pa
densities : mol / m³
molar enthalpies : kJ / mol
molar entropies : kJ / mol / K

Returns
-------
dict[str, list[float]]
    Keys: property names. Values: property for each state.

Notes
-----
xi: liquid molefraction of component i
yi: vapor molefraction of component i
i: component index according to order in parameters.
Type:      builtin_function_or_method

data_dia = pd.DataFrame(dia.to_dict())
data_dia.head()

	density liquid	density vapor	molar enthalpy vapor	temperature	molar enthalpy liquid	pressure	molar entropy liquid	molar entropy vapor
0	14125.988947	52.208491	22.140400	230.000000	3.376293	96625.278174	0.039106	0.120689
1	14118.006852	52.811929	22.135738	230.280462	3.383021	97830.133956	0.039135	0.120569
2	14110.010220	53.420767	22.131064	230.560924	3.389761	99046.729400	0.039164	0.120449
3	14101.999011	54.035036	22.126380	230.841386	3.396514	100275.143120	0.039193	0.120330
4	14093.973182	54.654773	22.121683	231.121849	3.403278	101515.453964	0.039221	0.120211

Once we have a dataframe, we can store our results or create a nicely looking plot:

def phase_plot(data, x, y):
    fig, ax = plt.subplots(figsize=(12, 6))
    if x != "pressure" and x != "temperature":
        xl = f"{x} liquid"
        xv = f"{x} vapor"
    else:
        xl = x
        xv = x
    if y != "pressure" and y != "temperature":
        yl = f"{y} liquid"
        yv = f"{y} vapor"
    else:
        yv = y
        yl = y
    sns.lineplot(data=data, x=xv, y=yv, ax=ax, label="vapor")
    sns.lineplot(data=data, x=xl, y=yl, ax=ax, label="liquid")
    ax.set_xlabel(x)
    ax.set_ylabel(y)
    ax.legend(frameon=False)
    sns.despine();

phase_plot(data_dia, "density", "temperature")

phase_plot(data_dia, "molar entropy", "temperature")

Mixtures

↑ Back to top

Fox mixtures, we have to add information about the composition, either as molar fraction, amount of substance per component, or as partial densities.

# propane, butane mixture
tc = np.array([369.96, 425.2]) * KELVIN
pc = np.array([4250000.0, 3800000.0]) * PASCAL
omega = np.array([0.153, 0.199])
molar_weight = np.array([44.0962, 58.123]) * GRAM / MOL

pr = PyPengRobinson(tc, pc, omega, molar_weight)
eos = EquationOfState.python(pr)

s = State(
    eos,
    temperature=300*KELVIN,
    pressure=1*BAR,
    molefracs=np.array([0.5, 0.5]),
    total_moles=MOL
)
s

temperature	density	molefracs
300.00000 K	40.96869 mol/m³	[0.50000, 0.50000]

As before, we can compute properties by calling methods on the State object. Some return vectors or matrices - for example the chemical potential and its derivative w.r.t amount of substance:

s.chemical_potential()

[-15625.3474516824, -12435.866602695127] J/mol

s.dmu_dni() / (KILO * JOULE / MOL**2)

array([[ 4.90827975, -0.10593968],
       [-0.10593968,  4.85467746]])

Phase equilibria can be built from different constructors. E.g. at critical conditions given composition:

s_cp = State.critical_point(
    eos,
    moles=np.array([0.5, 0.5])*MOL
)
s_cp

temperature	density	molefracs
401.65486 K	3.99952 kmol/m³	[0.50000, 0.50000]

Or at given temperature (or pressure) and composition for bubble and dew points.

vle = PhaseEquilibrium.bubble_point(
    eos,
    350*KELVIN,
    liquid_molefracs=np.array([0.5, 0.5])
)
vle

	temperature	density	molefracs
phase 1	350.00000 K	879.47505 mol/m³	[0.67625, 0.32375]
phase 2	350.00000 K	8.96382 kmol/m³	[0.50000, 0.50000]

Similar to pure substance phase diagrams, there is a constructor for binary systems.

vle = PhaseDiagram.binary_vle(eos, 350*KELVIN, npoints=50)

fig, ax = plt.subplots(1, 2, figsize=(18, 6))
# fig.title("T = 350 K, Propane (1), Butane (2)")
sns.lineplot(x=vle.liquid.molefracs[:,0], y=vle.liquid.pressure / BAR, ax=ax[0])
sns.lineplot(x=vle.vapor.molefracs[:,0], y=vle.vapor.pressure / BAR, ax=ax[0])
ax[0].set_xlabel(r"$x_1$, $y_1$")
ax[0].set_ylabel(r"$p$ / bar")
ax[0].set_xlim(0, 1)
ax[0].set_ylim(5, 35)
# ax[0].legend(frameon=False);

sns.lineplot(x=vle.liquid.molefracs[:,0], y=vle.vapor.molefracs[:,0], ax=ax[1])
sns.lineplot(x=np.linspace(0, 1, 10), y=np.linspace(0, 1, 10), color="black", alpha=0.3, ax=ax[1])
ax[1].set_xlabel(r"$x_1$")
ax[1].set_ylabel(r"$y_1$")
ax[1].set_xlim(0, 1)
ax[1].set_ylim(0, 1);

Comparison to Rust implementation

↑ Back to top

Implementing an equation of state in Python is nice for quick prototyping and development but when it comes to performance, implementing the equation of state in Rust is the way to go. For each non-cached call to the Helmholtz energy, we have to transition between Rust and Python with our Python implementation which generates quite some overhead.

Here are some comparisons between the Rust and our Pyhton implemenation:

# rust
from feos.cubic import PengRobinsonParameters
eos_rust = EquationOfState.peng_robinson(PengRobinsonParameters.from_json(["propane"], "peng-robinson.json"))

# python
tc = SIArray1(369.96 * KELVIN)
pc = SIArray1(4250000.0 * PASCAL)
omega = np.array([0.153])
molar_weight = SIArray1(44.0962 * GRAM / MOL)
eos_python = EquationOfState.python(PyPengRobinson(tc, pc, omega, molar_weight))

# let's first test if both actually yield the same results ;)
assert abs(State.critical_point(eos_python).pressure() / BAR - State.critical_point(eos_rust).pressure() / BAR) < 1e-13
assert abs(State.critical_point(eos_python).temperature / KELVIN - State.critical_point(eos_rust).temperature / KELVIN) < 1e-13

import timeit

time_python = timeit.timeit(lambda: State.critical_point(eos_python), number=2_500)
time_rust = timeit.timeit(lambda: State.critical_point(eos_rust), number=2_500)

rel_dev = (time_rust - time_python) / time_rust
print(f"Critical point for pure substance")
print(f"Python implementation is {'slower' if rel_dev < 0 else 'faster'} by a factor of {abs(time_python / time_rust):.0f}.")

Critical point for pure substance
Python implementation is slower by a factor of 48.

time_python = timeit.timeit(lambda: PhaseDiagram.pure(eos_python, 300*KELVIN, 100), number=100)
time_rust = timeit.timeit(lambda: PhaseDiagram.pure(eos_rust, 300*KELVIN, 100), number=100)

rel_dev = (time_rust - time_python) / time_rust
print(f"Phase diagram for pure substance")
print(f"Python implementation is {'slower' if rel_dev < 0 else 'faster'} by a factor of {abs(time_python / time_rust):.0f}.")

Phase diagram for pure substance
Python implementation is slower by a factor of 31.