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 ofnpoints
) ofPhaseEquilibrium
objects at the different temperatures,vapor
andliquid
: so-calledStateVec
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)
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
andPhaseEquilibrium.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]:
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);