Add VLLE finding for isobar tracing

Closes #46
usnistgov · Sep 6, 2023 · b58f416 · b58f416
1 parent 515d995
commit b58f416
Show file tree

Hide file tree

Showing 6 changed files with 374 additions and 34 deletions.
diff --git a/doc/source/algorithms/VLLE-p.ipynb b/doc/source/algorithms/VLLE-p.ipynb
@@ -0,0 +1,102 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "1ec37f01",
+   "metadata": {},
+   "source": [
+    "# VLLE\n",
+    "\n",
+    "Following the approach described in Bell et al.: https://doi.org/10.1021/acs.iecr.1c04703\n",
+    "\n",
+    "for the mixture of nitrogen + ethane, with the default thermodynamic model in teqp, which is the GERG-2008 mixing parameters (no departure function).\n",
+    "\n",
+    "Two traces are made, and the intersection is obtained, this gives you the VLLE solution."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "29a2031a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import teqp, numpy as np, matplotlib.pyplot as plt, pandas\n",
+    "import CoolProp.CoolProp as CP \n",
+    "\n",
+    "names = ['Nitrogen', 'Ethane']\n",
+    "model = teqp.build_multifluid_model(names, teqp.get_datapath())\n",
+    "pures = [teqp.build_multifluid_model([name], teqp.get_datapath()) for name in names]\n",
+    "p = 29e5 # Pa\n",
+    "\n",
+    "# Trace from both pure fluid endpoints\n",
+    "traces = []\n",
+    "for ipure in [1,0]:\n",
+    "    # Init at the pure fluid endpoint\n",
+    "    anc = pures[ipure].build_ancillaries()\n",
+    "    rhoLpure, rhoVpure = [CP.PropsSI('Dmolar','P',p,'Q',Q,names[ipure]) for Q in [0,1]]\n",
+    "    T = CP.PropsSI('T','P',p,'Q',0,names[ipure])\n",
+    "\n",
+    "    rhovecL = np.array([0.0, 0.0])\n",
+    "    rhovecV = np.array([0.0, 0.0])\n",
+    "    rhovecL[ipure] = rhoLpure\n",
+    "    rhovecV[ipure] = rhoVpure\n",
+    "    opt = teqp.TVLEOptions(); \n",
+    "    opt.p_termination = 1e8; \n",
+    "    opt.crit_termination=1e-4; \n",
+    "    opt.calc_criticality=True\n",
+    "    j = model.trace_VLE_isobar_binary(p, T, rhovecL, rhovecV)\n",
+    "    df = pandas.DataFrame(j)\n",
+    "    plt.plot(df['xL_0 / mole frac.'], df['T / K'])\n",
+    "    plt.plot(df['xV_0 / mole frac.'], df['T / K'])\n",
+    "    traces.append(j)\n",
+    "    \n",
+    "# Do the VLLE solving\n",
+    "for soln in model.find_VLLE_p_binary(traces):\n",
+    "    print(soln)\n",
+    "    T = soln['polished'][-1]\n",
+    "    print('rhovec / mol/m^3 | T / K')\n",
+    "    for rhovec in soln['polished'][0:3]:\n",
+    "        rhovec = np.array(rhovec)\n",
+    "        rhotot = sum(rhovec)\n",
+    "        x = rhovec/rhotot\n",
+    "        p = rhotot*model.get_R(x)*T*(1+model.get_Ar01(T, rhotot, x))\n",
+    "        \n",
+    "        plt.plot(x[0], T, 'X')\n",
+    "        print(rhovec, T)\n",
+    "        \n",
+    "#     # And also carry out the LLE trace for the two liquid phases\n",
+    "#     j = model.trace_VLE_isotherm_binary(T, np.array(soln['polished'][1]), np.array(soln['polished'][2]), opt)\n",
+    "#     df = pandas.DataFrame(j)\n",
+    "#     plt.plot(df['xL_0 / mole frac.'], df['pL / Pa'], 'k')\n",
+    "#     plt.plot(df['xV_0 / mole frac.'], df['pV / Pa'], 'k')\n",
+    "\n",
+    "# Plotting niceties\n",
+    "plt.ylim(top=280, bottom=100)\n",
+    "plt.gca().set(xlabel='$x_1$ / mole frac.', ylabel='$T$ / K', title='nitrogen(1) + ethane(2)')\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "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
+}