Skip to content

Commit

Permalink
Merge pull request #74 from usnistgov/polarizable
Browse files Browse the repository at this point in the history 
Polarizable model of Gray et al.
  • Loading branch information
@ianhbell
ianhbell authored Oct 20, 2023
2 parents 561c7ec + 3acf828 commit 18256d8
Show file tree
Hide file tree
Showing 4 changed files with 737 additions and 51 deletions.
214 changes: 214 additions & 0 deletions doc/source/models/SAFT-VR-Mie-Dipolar-Polarizable.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,214 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "8482721f-0615-4b73-83ee-187b2f890521",
"metadata": {},
"source": [
"# SAFT-VR-Mie with polar contributions"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d1842386",
"metadata": {},
"outputs": [],
"source": [
"import teqp\n",
"teqp.__version__"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "16f05ec8",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np, io\n",
"import matplotlib.pyplot as plt, pandas\n",
"import math, json"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "be8268a7",
"metadata": {},
"outputs": [],
"source": [
"# These values are not important, get something on the right order of magnitude\n",
"ek = 100 # [K]\n",
"sigma_m = 1e-10\n",
" \n",
"N_A = 6.022e23\n",
"fig, (ax1, ax2) = plt.subplots(2, 1)\n",
"\n",
"kB = 1.380649e-23 # Boltzmann's constant, J/K\n",
"epsilon_0 = 8.8541878128e-12 # Vacuum permittivity\n",
"k_e = 1.0/(4.0*np.pi*epsilon_0*sigma_m**3)\n",
"\n",
"polar_model = 'GrayGubbins+GubbinsTwu'\n",
"\n",
"for mustar in [1, 2]:\n",
" x,TT,DD = [],[],[]\n",
" for alphastar in [0.0, 0.03, 0.06]:\n",
"\n",
" alpha_m3 = alphastar*sigma_m**3\n",
"\n",
" rhostar_guess = 0.27\n",
" Tstar_guess = 1.5\n",
" mu_Cm = (ek*kB/k_e)**0.5*mustar\n",
" j = {\n",
" \"kind\": 'SAFT-VR-Mie',\n",
" \"model\": {\n",
" \"polar_model\": polar_model,\n",
" \"polar_flags\": {\n",
" \"polarizable\": {\n",
" \"alpha_symm / m^3\": [alpha_m3], \n",
" \"alpha_asymm / m^3\": [0.0]\n",
" }\n",
" },\n",
" \"coeffs\": [{\n",
" \"name\": \"PolarizableStockmayer\",\n",
" \"BibTeXKey\": \"me\",\n",
" \"m\": 1.0,\n",
" \"epsilon_over_k\": ek, # [K]\n",
" \"sigma_m\": sigma_m,\n",
" \"lambda_r\": 12.0,\n",
" \"lambda_a\": 6.0,\n",
" \"mu_Cm\": mu_Cm,\n",
" \"nmu\": 1.0\n",
" }]\n",
" }\n",
" }\n",
" model = teqp.make_model(j)\n",
"\n",
" T, rho = model.solve_pure_critical(Tstar_guess*ek, rhostar_guess/(N_A*sigma_m**3))\n",
" # Store the values\n",
" x.append(alphastar)\n",
" TT.append(T/ek)\n",
" DD.append(rho*N_A*sigma_m**3)\n",
" # Update the guess for the next calculation\n",
" Tstar_guess = TT[-1]\n",
" rhostar_guess = DD[-1]\n",
"# print(TT[-1], DD[-1])\n",
"\n",
" ax1.plot(x, TT, label=f'$\\mu^*={mustar}$')\n",
" ax2.plot(x, DD)\n",
" \n",
"ax1.legend(loc='best')\n",
"ax1.set(ylabel=r'$T^*$')\n",
"ax2.set(xlabel=r'$\\alpha^*$', ylabel=r'$\\rho^*$')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5d013035",
"metadata": {},
"outputs": [],
"source": [
"s = io.StringIO(\"\"\"# Kiyohara JCP 1999; doi: 10.1063/1.473082\n",
"alpha* T* mu* p* rhog* rhol* ug* ul* deltah*\n",
"0.00 1.60 -7.177 0.0224 0.0222 0.726 -1.16 -10.16 9.97\n",
"0.00 1.65 -7.078 0.0300 0.0273 0.706 -1.28 -9.89 9.67\n",
"0.00 1.70 -6.986 0.0388 0.0335 0.682 -1.40 -9.58 9.28\n",
"0.00 1.75 -6.900 0.0490 0.0411 0.654 -1.53 -9.22 8.82\n",
"0.00 1.80 -6.822 0.0607 0.0507 0.626 -1.67 -8.87 8.30\n",
"0.00 1.85 -6.750 0.0741 0.0634 0.599 -1.86 -8.53 7.72\n",
"0.00 1.90 -6.683 0.0896 0.0811 0.569 -2.13 -8.17 6.98\n",
"0.03 1.70 -7.893 0.0233 0.0183 0.758 -1.11 -11.48 11.62\n",
"0.03 1.75 -7.783 0.0301 0.0222 0.741 -1.22 -11.21 11.30\n",
"0.03 1.80 -7.679 0.0379 0.0270 0.720 -1.34 -10.91 10.91\n",
"0.03 1.85 -7.582 0.0469 0.0327 0.697 -1.46 -10.58 10.48\n",
"0.03 1.90 -7.492 0.0572 0.0396 0.674 -1.58 -10.25 10.02\n",
"0.03 1.95 -7.407 0.0690 0.0480 0.650 -1.71 -9.92 9.53\n",
"0.03 2.00 -7.329 0.0823 0.0587 0.624 -1.87 -9.56 8.95\n",
"0.03 2.05 -7.255 0.0974 0.0730 0.593 -2.10 -9.16 8.21\n",
"0.03 2.10 -7.187 0.1146 0.0927 0.556 -2.44 -8.69 7.27\n",
"0.06 2.00 -8.695 0.0357 0.0232 0.761 -1.19 -13.00 13.30\n",
"0.06 2.05 -8.581 0.0440 0.0275 0.749 -1.31 -12.75 12.98\n",
"0.06 2.10 -8.471 0.0535 0.0325 0.732 -1.44 -12.44 12.57\n",
"0.06 2.15 -8.367 0.0641 0.0385 0.709 -1.58 -12.07 12.07\n",
"0.06 2.20 -8.270 0.0761 0.0455 0.686 -1.71 -11.69 11.54\n",
"0.06 2.25 -8.178 0.0895 0.0539 0.663 -1.86 -11.33 11.00\n",
"0.06 2.30 -8.092 0.1044 0.0644 0.639 -2.05 -10.95 10.36\n",
"0.06 2.35 -8.010 0.1211 0.0784 0.609 -2.30 -10.51 9.55\n",
"0.06 2.40 -7.934 0.1397 0.0969 0.573 -2.66 -9.99 8.53\"\"\")\n",
"\n",
"df = pandas.read_csv(s, sep='\\s+', engine='python', comment='#')\n",
"\n",
"mustar = 2.0\n",
"for alphastar, gp in df.groupby('alpha*'):\n",
" \n",
" alpha_m3 = alphastar*sigma_m**3 \n",
" \n",
" j = {\n",
" \"kind\": 'SAFT-VR-Mie',\n",
" \"model\": {\n",
" \"polar_model\": polar_model,\n",
" \"polar_flags\": {\n",
" \"polarizable\": {\n",
" \"alpha_symm / m^3\": [alpha_m3], \n",
" \"alpha_asymm / m^3\": [0.0]\n",
" }\n",
" },\n",
" \"coeffs\": [{\n",
" \"name\": \"PolarizableStockmayer\",\n",
" \"BibTeXKey\": \"me\",\n",
" \"m\": 1.0,\n",
" \"epsilon_over_k\": ek, # [K]\n",
" \"sigma_m\": sigma_m,\n",
" \"lambda_r\": 12.0,\n",
" \"lambda_a\": 6.0,\n",
" \"mu_Cm\": mu_Cm,\n",
" \"nmu\": 1.0\n",
" }]\n",
" }\n",
" }\n",
" model = teqp.make_model(j)\n",
" Tc, rhoc = model.solve_pure_critical(Tstar_guess*ek, rhostar_guess/(N_A*sigma_m**3))\n",
" anc = teqp.build_ancillaries(model, Tc, rhoc, Tc/2.0, {})\n",
" Tvec = np.linspace(Tc/2.0, Tc, 1000)\n",
" \n",
" line, = plt.plot(gp['rhol*'], gp['T*'], 'o')\n",
" plt.plot(gp['rhog*'], gp['T*'], 'o', color=line.get_color())\n",
" \n",
" RHOL = np.array([anc.rhoL(T) for T in Tvec])\n",
" RHOV = np.array([anc.rhoV(T) for T in Tvec])\n",
" \n",
" plt.plot(RHOL*N_A*sigma_m**3, Tvec/ek, '-', color=line.get_color(), label=rf'$\\alpha^*$: {alphastar}')\n",
" plt.plot(RHOV*N_A*sigma_m**3, Tvec/ek, '-', color=line.get_color())\n",
"\n",
"plt.title('Comparison with the MC data of Kiyohara')\n",
"plt.xlabel(r'$\\rho^*$')\n",
"plt.ylabel(r'$T^*$')\n",
"plt.legend();"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
16 changes: 15 additions & 1 deletion include/teqp/constants.hpp
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,21 @@
#pragma once

namespace teqp {
const double N_A = 6.02214076e23; ///< Avogadro's number

namespace constants{

const double R_CODATA2002 = 8.314472; // CODATA 2002: https://doi.org/10.1103/RevModPhys.77.1
const double R_CODATA2006 = 8.314472; // CODATA 2006: https://doi.org/10.1103/RevModPhys.80.633
const double R_CODATA2010 = 8.3144621; // CODATA 2010: https://doi.org/10.1103/RevModPhys.84.1527
const double R_CODATA2017 = 8.31446261815324; // CODATA 2017: https://doi.org/10.1103/RevModPhys.84.1527

const double N_A = 6.02214076e23; ///< Avogadro's number
const double k_B = 1.380649e-23; ///< Boltzmann constant
const double epsilon_0 = 8.8541878128e-12; // https://en.wikipedia.org/wiki/Vacuum_permittivity, in F/m, or C^2⋅N^−1⋅m^−2
const double PI = 3.141592653589793238462643383279502884197;
const double k_e = 1.0/(4.0*PI*epsilon_0); // coulomb constant, with units of N m^2 / C^2
}
using constants::N_A;

///< Gas constant, according to CODATA 2019, in the given number type
template<typename NumType>
Expand Down
Loading

0 comments on commit 18256d8

Please sign in to comment.