Pure substance phase diagrams

Goal of this notebook

  • Learn how to generate and work with phase diagrams.

  • Learn what PhaseEquilibrium objects are and how to use them.

[3]:
from feos.si import *
from feos.pcsaft import *
from feos.eos import *

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_context('talk')
sns.set_palette('Dark2')
sns.set_style('ticks')
colors = sns.palettes.color_palette('Dark2', 8)

Read parameters from json for a single substance and generate PcSaft object

[4]:
parameters = PcSaftParameters.from_json(
    ['hexane'],
    '../parameters/pcsaft/gross2001.json'
)
parameters
[4]:

component

molarweight

\(m\)

\(\sigma\)

\(\varepsilon\)

\(\mu\)

\(Q\)

\(\kappa_{AB}\)

\(\varepsilon_{AB}\)

\(N_A\)

\(N_B\)

hexane

86.177

3.0576

3.7983

236.77

0

0

0

0

1

1

[5]:
pcsaft = EquationOfState.pcsaft(parameters)

The PhaseDiagram object

A PhaseDiagram object contains multiple thermodyanmic states at phase equilibrium. For a single substance it can be constructed via the PhaseDiagram.pure method by providing

  • an equation of state,

  • the minimum temperature, and

  • the number of points to compute.

The points are evenly distributed between the defined minimum temperature and the critical temperature which is either computed (default) or can be provided via the critical_temperature argument.

Furthermore, we can define options for the numerics: - max_iter to define the maximum number of iterations for the phase equilibrium calculations, and - tol to set the tolerance used for deterimation of phase equilibria.

A Verbosity object can be provided via the verbosity argument to print intermediate calculations to the screen.

Starting from the minimum temperature, phase equilibria are computed in sequence using prior results (i.e. at prior temperature) as input for the next iteration.

[6]:
phase_diagram = PhaseDiagram.pure(pcsaft, min_temperature=200*KELVIN, npoints=501)
[7]:
phase_diagram.vapor[0]
[7]:

temperature

density

200.00000 K

12.21572 mmol/m³

Stored information

A PhaseDiagram object contains the following fields:

  • states: a list (with length of npoints) of PhaseEquilibrium objects at the different temperatures,

  • vapor and liquid: so-called StateVec objects that can be used to compute properties for the vapor and liquid phase, respectively.

The to_dict method can be used to conveniently generate a pandas.DataFrame object (see below).

Building a pandas.DataFrame: the to_dict method

The PhaseDiagramPure object contains some physical properties (such as densities, temperatures and pressures) as well as the PhaseEquilibrium objects at each temperature.

Before we take a look at these objects, a useful tool when working with phase diagrams is the to_dict method. It generates a Python dictionary of some properties (with hard-coded units, see the docstring of to_dict). This dictionary can readily be used to generate a pandas.DataFrame.

[8]:
df = pd.DataFrame(phase_diagram.to_dict())
df.head()
[8]:
molar enthalpy vapor temperature density liquid density vapor molar enthalpy liquid pressure molar entropy vapor molar entropy liquid
0 -16.735400 200.000000 8539.589591 0.012216 -54.057307 20.312763 0.003135 -0.183475
1 -16.639176 200.638669 8532.663964 0.013078 -53.920803 21.816282 0.003021 -0.182794
2 -16.542786 201.277337 8525.749142 0.013994 -53.784188 23.418684 0.002912 -0.182114
3 -16.446230 201.916006 8518.845002 0.014967 -53.647463 25.125608 0.002806 -0.181436
4 -16.349508 202.554674 8511.951422 0.015999 -53.510626 26.942961 0.002703 -0.180759
[9]:
fig, ax = plt.subplots(1, 2, figsize=(15, 6))
ax[0].set_title(f"saturation pressure of {parameters.pure_records[0].identifier.name}")
sns.lineplot(y=df.pressure, x=1.0/df.temperature, ax=ax[0])

# axis and styling
ax[0].set_yscale('log')
ax[0].set_xlabel(r'$\frac{1}{T}$ / K$^{-1}$');
ax[0].set_ylabel(r'$p$ / Pa');
ax[0].set_xlim(0.002, 0.005)
ax[0].set_ylim(1e1, 1e7)

ax[1].set_title(r"$T$-$\rho$-diagram of {}".format(parameters.pure_records[0].identifier.name))
sns.lineplot(y=df.temperature, x=df['density vapor'], ax=ax[1], color=colors[0])
sns.lineplot(y=df.temperature, x=df['density liquid'], ax=ax[1], color=colors[0])

# axis and styling
ax[1].set_ylabel(r'$T$ / K');
ax[1].set_xlabel(r'$\rho$ / mol/m³');
ax[1].set_ylim(200, 600)
ax[1].set_xlim(0, 9000)

sns.despine(offset=10)
../../../_images/tutorials_eos_pcsaft_pcsaft_phase_diagram_11_0.png

The PhaseEquilibrium object

A PhaseEquilibrium object contains two thermodynamic states (State objects) that are in thermal, mechanical and chemical equilibrium. The PhaseEquilibrium.pure constructor can be used to compute a phase equilibrium for a single substance. We have to provide the equation of state and either temperature or pressure. Optionally, we can also provide PhaseEquilibrium object which can be used as starting point for the calculation which possibly speeds up the computation.

[10]:
vle = PhaseEquilibrium.pure(
    pcsaft,
    temperature_or_pressure=300.0*KELVIN
)
vle
[10]:

temperature

density

phase 1

300.00000 K

8.86860 mol/m³

phase 2

300.00000 K

7.51850 kmol/m³

The two equilibrium states can be extracted via the liquid and vapor getters, respectively. Returned are State objects, for which we can now compute any property that is available for the State object and the given equation of state.

[11]:
liquid = vle.liquid
vapor = vle.vapor

assert(abs((liquid.pressure() - vapor.pressure()) / BAR) < 1e-10)
print(f'saturation pressure      p_sat(T = {liquid.temperature}) = {liquid.pressure() / BAR:6.2f} bar')
print(f'enthalpy of vaporization h_lv (T = {liquid.temperature}) = {(vapor.specific_enthalpy() - liquid.specific_enthalpy()) / (KILO*JOULE/KILOGRAM):6.2f} kJ/kg')
saturation pressure      p_sat(T = 300 K) =   0.22 bar
enthalpy of vaporization h_lv (T = 300 K) = 365.83 kJ/kg

If you want to compute a boiling temperature or saturation pressure without needing the PhaseEquilibrium object you can use the

  • PhaseEquilibrium.boiling_temperature and

  • PhaseEquilibrium.vapor_pressure

methods. Note that these methods return lists (even for pure substance systems) where each entry contains the pure substance property.

[12]:
PhaseEquilibrium.boiling_temperature(pcsaft, liquid.pressure())[0]
[12]:
$300\,\mathrm{K}$

The states method returns all PhaseEquilibrium objects from PhaseDiagram

Once a PhaseDiagram object is created, we can access all underlying PhaseEquilibrium objects via the states field. In the following cell, we compute the enthalpy of vaporization by iterating through all states and calling the specific_enthalpy method on the vapor and liquid states of the PhaseEquilibrium object, respectively.

Note that this is merely an example to show how to compute any property of the states. The total value of enthalpy of vaporization may not be correct, depending on the ideal gas model used.

[13]:
# Add enthalpy of vaporization to dataframe
df['hlv'] = [
    (vle.vapor.specific_enthalpy() - vle.liquid.specific_enthalpy())
    / (KILO * JOULE / KILOGRAM)
    for vle in phase_diagram.states
]
[14]:
fig, ax = plt.subplots(figsize=(7, 6))
sns.lineplot(y=df.hlv, x=df.temperature, ax=ax)

# axis and styling
ax.set_xlabel(r'$T$ / K');
ax.set_ylabel(r'$\Delta_{LV} h$ / kJ / kg');
ax.set_xlim(200, 600)
ax.set_ylim(0, 500)

sns.despine(offset=10);
../../../_images/tutorials_eos_pcsaft_pcsaft_phase_diagram_20_0.png