{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Comparison of equations of state for model fluids\n",
"\n",
"In this notebook, we take a look at various model fluid equations of state and compare them.\n",
"These EoS are:\n",
"\n",
"- The PeTS equation of state (from `feos.pets`),\n",
"- the uv-theory equation of state (from `feos.uvtheory`), \n",
"- the uv-B3-theory equation of state (from `feos.uvtheory` using `Perturbation.`), and\n",
"- the Thol equation of state (implemented in this notebook as Python class).\n",
"\n",
"The PeTS equation of state is for a Lennard-Jones fluid with cut off distance of $r_c = 2.5 \\sigma$ and a potential shift while the uv-theory and Thol equation of state model the full Lennard-Jones potential. Note that the uv-theory was developed to model Mie potentials and not primarily adjusted to the Lennard-Jones fluid. Thus, we don't expect similar results when comparing uv-theory and Thol to PeTS."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from feos.si import *\n",
"from feos.uvtheory import UVTheoryParameters, Perturbation\n",
"from feos.pets import PetsParameters\n",
"from feos.eos import State, PhaseEquilibrium, PhaseDiagram, EquationOfState, Contributions\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"import pandas as pd\n",
"\n",
"import warnings\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"sns.set_context('talk')\n",
"sns.set_palette(\"hls\", 8)\n",
"sns.set_style('ticks')\n",
"colors = sns.palettes.color_palette('Dark2', 8)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Implement the Thol equation of state as Python class"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Parameters for Thol EoS:\n",
"A = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 2.0,\n",
" 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0])\n",
"N = np.array([0.005208073, 2.186252000, -2.161016000, 1.452700000, -2.041792000, 0.186952860, -0.090988445, \n",
" -0.497456100, 0.109014310, -0.800559220, -0.568839000, -0.620862500, -1.466717700, 1.891469000, \n",
" -0.138370100, -0.386964500, 0.126570200, 0.605781000, 1.179189000, -0.477326790, -9.921857500, -0.574793200, 0.003772923])\n",
"T = np.array([1.000, 0.320, 0.505, 0.672, 0.843, 0.898, 1.294, 2.590, 1.786, 2.770, 1.786,\n",
" 1.205, 2.830, 2.548, 4.650, 1.385, 1.460, 1.351, 0.660, 1.496, 1.830, 1.616, 4.970])\n",
"D = np.array([4.0, 1.0, 1.0, 2.0, 2.0, 3.0, 5.0, 2.0, 2.0, 3.0, 1.0,\n",
" 1.0, 1.0, 1.0, 2.0, 3.0, 3.0, 2.0, 1.0, 2.0, 3.0, 1.0, 1.0])\n",
"L = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 1.0, 2.0, 2.0,\n",
" 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])\n",
"ETA = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.067, 1.522,\n",
" 8.82, 1.722, 0.679, 1.883, 3.925, 2.461, 28.2, 0.753, 0.82])\n",
"BETA = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.625,\n",
" 0.638, 3.91, 0.156, 0.157, 0.153, 1.16, 1.73, 383, 0.112, 0.119])\n",
"GAMMA = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.71,\n",
" 0.86, 1.94, 1.48, 1.49, 1.945, 3.02, 1.11, 1.17, 1.33, 0.24])\n",
"EPSILON = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2053,\n",
" 0.409, 0.6, 1.203, 1.829, 1.397, 1.39, 0.539, 0.934, 2.369, 2.43])\n",
"\n",
"class Thol:\n",
" def __init__(self):\n",
" self.sigma = 3.7039 * ANGSTROM\n",
" self.eps_k = 150.03 * KELVIN\n",
" self.tc = 1.32 * self.eps_k / KELVIN\n",
" self.rhoc = 0.31 / self.sigma**3 * ANGSTROM**3\n",
" \n",
" def components(self): \n",
" return 1\n",
" \n",
" def subset(self, components):\n",
" return self\n",
" \n",
" def molar_weight(self):\n",
" return np.array([1.0])\n",
" \n",
" def max_density(self, moles):\n",
" return 0.04\n",
" \n",
" def helmholtz_energy(self, state):\n",
" \"\"\"\n",
" state (StateHD):\n",
" temperature in Kelvin als Float, Dual oder HD oder HD3, HDD, HDD3,\n",
" partial_density in # / Angstrom^3\n",
" volume in Angstrom^3\n",
" moles in mol\n",
" \"\"\"\n",
" tau = self.tc / state.temperature\n",
" delta = np.sum(state.partial_density) / self.rhoc\n",
" a = 0.0\n",
" \n",
" for i in range(6):\n",
" a = a + N[i] * delta**D[i] * tau**T[i]\n",
" for i in range(6, 12):\n",
" a = a + N[i] * delta**D[i] * tau**T[i] * np.exp(-1.0 * delta**L[i])\n",
" for i in range(12, 23):\n",
" a = a + N[i] * delta**D[i] * tau**T[i] * np.exp(- 1.0 * ETA[i] * (delta - EPSILON[i])**2.0 - 1.0 * BETA[i] * (tau - 1.0 * GAMMA[i])**2.0)\n",
" return a * np.sum(state.moles)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Instantiate all equations of state"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"sigma = 3.7039\n",
"sigma_a = sigma * ANGSTROM\n",
"eps_k = 150.03\n",
"eps_k_k = eps_k * KELVIN\n",
"\n",
"parameters = UVTheoryParameters.new_simple(12.0, 6.0, sigma, eps_k)\n",
"uvtheory_wca = EquationOfState.uvtheory(parameters)\n",
"uvtheory_bh = EquationOfState.uvtheory(parameters, perturbation=Perturbation.BarkerHenderson)\n",
"uvtheory_b3 = EquationOfState.uvtheory(parameters, perturbation=Perturbation.WeeksChandlerAndersenB3)\n",
"\n",
"parameters = PetsParameters.from_values(sigma, eps_k)\n",
"pets = EquationOfState.pets(parameters)\n",
"\n",
"thol = EquationOfState.python_residual(Thol())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compare Helmholtz energies at single state point"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"s1 = State(uvtheory_wca, temperature=300*KELVIN, pressure=1*BAR)\n",
"s2 = State(uvtheory_b3, temperature=300*KELVIN, pressure=1*BAR)\n",
"s3 = State(thol, temperature=300*KELVIN, pressure=1*BAR)\n",
"s4 = State(pets, temperature=300*KELVIN, pressure=1*BAR)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.0016213224956115704"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s1.molar_helmholtz_energy(Contributions.Residual) / (RGAS * 300 * KELVIN)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.0016138652202999685"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s2.molar_helmholtz_energy(Contributions.Residual) / (RGAS * 300 * KELVIN)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.001614415096260159"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s3.molar_helmholtz_energy(Contributions.Residual) / (RGAS * 300 * KELVIN)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-0.0009407040108813862"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s4.molar_helmholtz_energy(Contributions.Residual) / (RGAS * 300 * KELVIN)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compare critical points\n",
"- Reference data for the Lennard-Jones from Potoff and Panagiotopoulos (MC): $T_c^*=1.3120$, $\\rho_c^*=0.316$, $p_c^*=0.1279$"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"uvtheory_b3:\n",
"Tc*=1.31329, Error: 0.099%\n",
"pc*=0.1286, Error: 0.548%\n",
"rhoc*=0.31425, Error: 0.553%\n",
"\n",
"uvtheory_wca:\n",
"Tc*=1.30977, Error: 0.17%\n",
"pc*=0.12448, Error: 2.674%\n",
"rhoc*=0.30233, Error: 4.327%\n",
"\n",
"uvtheory_bh:\n",
"Tc*=1.32083, Error: 0.673%\n",
"pc*=0.13486, Error: 5.442%\n",
"rhoc*=0.31865, Error: 0.838%\n",
"\n",
"thol:\n",
"Tc*=1.32, Error: 0.61%\n",
"pc*=0.13006, Error: 1.689%\n",
"rhoc*=0.3132, Error: 0.886%\n",
"\n"
]
}
],
"source": [
"def evaluate_critical_points():\n",
" names = ['uvtheory_b3', 'uvtheory_wca', 'uvtheory_bh', 'thol']\n",
" print('')\n",
" i = 0\n",
" for eos in [uvtheory_b3, uvtheory_wca, uvtheory_bh, thol]:\n",
" rep = 12\n",
" att = 6.0\n",
" sigma = 3.7039 # in Angstrom\n",
" eps_k = 150.03 # eps / kB in K\n",
" print('{}:'.format(names[i]))\n",
" Tc = State.critical_point(eos).temperature / eps_k_k\n",
" print('Tc*={}, Error: {}%'.format(np.round(Tc, 5), np.round(np.abs(Tc - 1.3120 ) / 1.3120 * 100 ,3)))\n",
" pc = State.critical_point(eos).pressure() * (sigma * ANGSTROM)**3 / KB / eps_k_k\n",
" print('pc*={}, Error: {}%'.format(np.round(pc, 5), np.round(np.abs(pc - 0.1279 ) / 0.1279 * 100 ,3)))\n",
" rhoc = State.critical_point(eos).density * (sigma * ANGSTROM)**3 * NAV\n",
" print('rhoc*={}, Error: {}%'.format(np.round(rhoc, 5), np.round(np.abs(rhoc - 0.316) / 0.316 * 100 ,3)))\n",
" i += 1\n",
" print('')\n",
" return\n",
"\n",
"evaluate_critical_points()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compare Phase diagrams"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 17.9 ms, sys: 1.2 ms, total: 19.1 ms\n",
"Wall time: 23 ms\n"
]
}
],
"source": [
"%%time\n",
"vle_uv = PhaseDiagram.pure(uvtheory_wca, 150*KELVIN, 500)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 16.5 ms, sys: 852 µs, total: 17.4 ms\n",
"Wall time: 19.4 ms\n"
]
}
],
"source": [
"%%time\n",
"vle_uv_bh = PhaseDiagram.pure(uvtheory_bh, 150*KELVIN, 500)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 26.5 ms, sys: 1.29 ms, total: 27.8 ms\n",
"Wall time: 29.7 ms\n"
]
}
],
"source": [
"%%time\n",
"vle_uv_b3 = PhaseDiagram.pure(uvtheory_b3, 150*KELVIN, 500)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 971 ms, sys: 5.32 ms, total: 976 ms\n",
"Wall time: 976 ms\n"
]
}
],
"source": [
"%%time\n",
"vle_thol = PhaseDiagram.pure(thol, 150*KELVIN, 500)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 16.3 ms, sys: 1.41 ms, total: 17.7 ms\n",
"Wall time: 18.3 ms\n"
]
}
],
"source": [
"%%time\n",
"vle_pets = PhaseDiagram.pure(pets, 150*KELVIN, 500)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" t | \n",
" rhov | \n",
" rhol | \n",
" p | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 0.70 | \n",
" 0.002012 | \n",
" 0.843237 | \n",
" 0.001381 | \n",
"
\n",
" \n",
" | 10 | \n",
" 0.71 | \n",
" 0.002284 | \n",
" 0.839106 | \n",
" 0.001586 | \n",
"
\n",
" \n",
" | 20 | \n",
" 0.72 | \n",
" 0.002584 | \n",
" 0.834868 | \n",
" 0.001815 | \n",
"
\n",
" \n",
" | 30 | \n",
" 0.73 | \n",
" 0.002912 | \n",
" 0.830543 | \n",
" 0.002069 | \n",
"
\n",
" \n",
" | 40 | \n",
" 0.74 | \n",
" 0.003271 | \n",
" 0.826147 | \n",
" 0.002350 | \n",
"
\n",
" \n",
"
\n",
"