{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "df447d16",
   "metadata": {},
   "source": [
    "# General cubics\n",
    "\n",
    "The reduced residual Helmholtz energy for the main cubic EOS (van der Waals, Peng-Robinson, and Soave-Redlich-Kwong) can be written in a common form (see https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7365965/)\n",
    "\n",
    "$$ \\alpha^r = \\psi^{(-)} - \\dfrac{\\tau a_m}{RT_r } \\psi^{(+)} $$\n",
    "\n",
    "$$ \\psi^{(-)} =-\\ln(1-b_m\\rho ) $$\n",
    "\n",
    "$$ \\psi^{(+)} = \\dfrac{\\ln\\left(\\dfrac{\\Delta_1 b_m\\rho+1}{\\Delta_2b_m\\rho+1}\\right)}{b_m(\\Delta_1-\\Delta_2)} $$\n",
    "\n",
    "with the constants given by:\n",
    "\n",
    "* vdW: $\\Delta_1=0$, $\\Delta_2=0$\n",
    "* SRK: $\\Delta_1=1$, $\\Delta_2=0$\n",
    "* PR: $\\Delta_1=1+\\sqrt{2}$, $\\Delta_2=1-\\sqrt{2}$\n",
    "\n",
    "The quantities $a_m$ and $b_m$ are described (with exact solutions for the numerical coefficients) for each of these EOS in https://pubs.acs.org/doi/abs/10.1021/acs.iecr.1c00847.\n",
    "\n",
    "The models in teqp are instantiated based on knowledge of the critical temperature, pressure, and acentric factor.  Thereafter all quantities are obtained from derivatives of $\\alpha^r$."
   ]
  },
  {
   "cell_type": "raw",
   "id": "a65e7316",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "The Python class is here: :py:class:`GeneralizedCubic <teqp.teqp.GeneralizedCubic>`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f2564948",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-15T22:39:58.033948Z",
     "iopub.status.busy": "2024-03-15T22:39:58.033617Z",
     "iopub.status.idle": "2024-03-15T22:39:58.048341Z",
     "shell.execute_reply": "2024-03-15T22:39:58.047763Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.19.1'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import teqp\n",
    "teqp.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9a884dab",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-15T22:39:58.051063Z",
     "iopub.status.busy": "2024-03-15T22:39:58.050684Z",
     "iopub.status.idle": "2024-03-15T22:39:58.148895Z",
     "shell.execute_reply": "2024-03-15T22:39:58.148347Z"
    }
   },
   "outputs": [],
   "source": [
    "import json\n",
    "import CoolProp.CoolProp as CP\n",
    "\n",
    "# Values taken from http://dx.doi.org/10.6028/jres.121.011\n",
    "Tc_K = [ 190.564, 154.581, 150.687 ]\n",
    "pc_Pa = [ 4599200, 5042800, 4863000 ]\n",
    "acentric = [ 0.011, 0.022, -0.002 ]\n",
    "\n",
    "# Instantiate Peng-Robinson model\n",
    "modelPR = teqp.canonical_PR(Tc_K, pc_Pa, acentric)\n",
    "\n",
    "# Instantiate Soave-Redlich-Kwong model\n",
    "modelSRK = teqp.canonical_SRK(Tc_K, pc_Pa, acentric)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "75961e5b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-15T22:39:58.152150Z",
     "iopub.status.busy": "2024-03-15T22:39:58.151219Z",
     "iopub.status.idle": "2024-03-15T22:39:58.156906Z",
     "shell.execute_reply": "2024-03-15T22:39:58.156442Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Delta1': 2.414213562373095,\n",
       " 'Delta2': -0.41421356237309515,\n",
       " 'OmegaA': 0.4572355289213822,\n",
       " 'OmegaB': 0.07779607390388846,\n",
       " 'kind': 'Peng-Robinson'}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# And you can get information about the model in JSON format \n",
    "# from the get_meta function\n",
    "modelPR.get_meta()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "615c7830",
   "metadata": {},
   "source": [
    "## Adjusting k_ij\n",
    "\n",
    "Fine-tuned values of $k_{ij}$ can be provided when instantiating the model, for Peng-Robinson and SRK.  A complete matrix of all the $k_{ij}$ values must be provided. This allows for asymmetric mixing models in which $k_{ij}\\neq k_{ji}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a8c87890",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-15T22:39:58.160114Z",
     "iopub.status.busy": "2024-03-15T22:39:58.159271Z",
     "iopub.status.idle": "2024-03-15T22:39:58.165513Z",
     "shell.execute_reply": "2024-03-15T22:39:58.165081Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<teqp.teqp.AbstractModel at 0x7894f86dd850>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k_12 = 0.01\n",
    "kmat = [[0,k_12,0],[k_12,0,0],[0,0,0]]\n",
    "teqp.canonical_PR(Tc_K, pc_Pa, acentric, kmat)\n",
    "teqp.canonical_SRK(Tc_K, pc_Pa, acentric, kmat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3a62795",
   "metadata": {},
   "source": [
    "## Superancillary\n",
    "\n",
    "The superancillary equation gives the co-existing liquid and vapor (orthobaric) densities as a function of temperature. The set of Chebyshev expansions was developed in https://pubs.acs.org/doi/abs/10.1021/acs.iecr.1c00847 .  These superancillary equations are more accurate than iterative calculations in double precision arithmetic and also at least 10 times faster to calculate, and cannot fail in iterative routines, even extremely close to the critical point. \n",
    "\n",
    "The superancillary equation is only exposed for pure fluids to remove ambiguity when considering mixtures.  The returned tuple is the liquid and vapor densities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1b02ef72",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-15T22:39:58.168628Z",
     "iopub.status.busy": "2024-03-15T22:39:58.167794Z",
     "iopub.status.idle": "2024-03-15T22:39:58.173489Z",
     "shell.execute_reply": "2024-03-15T22:39:58.173060Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(30846.392909514052, 42.480231719002326)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "teqp.canonical_PR([Tc_K[0]], [pc_Pa[0]], [acentric[0]]).superanc_rhoLV(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da00ff20",
   "metadata": {},
   "source": [
    "## a and b\n",
    "\n",
    "For the cubic EOS, it can be useful to obtain the a and b parameters directly. The b parameter is particularly useful because 1/b is the maximum allowed density in the EOS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b5f6aaea",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-15T22:39:58.176645Z",
     "iopub.status.busy": "2024-03-15T22:39:58.175807Z",
     "iopub.status.idle": "2024-03-15T22:39:58.181003Z",
     "shell.execute_reply": "2024-03-15T22:39:58.180474Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.1874177858906821, 2.1984349667726406e-05)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "z = np.array([0.3, 0.4, 0.3])\n",
    "modelPR.get_a(140, z), modelPR.get_b(140, z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6de8071a",
   "metadata": {},
   "source": [
    "## alpha functions\n",
    "\n",
    "It can be advantageous to modify the alpha function to allow for more accurate handling of the attractive interactions. Coefficients are tabulated for many species in https://pubs.acs.org/doi/10.1021/acs.jced.7b00967 for the Peng-Robinson EOS with Twu alpha function and the values from the SI of that paper are in the csv file next to this file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0cfc2bf3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-15T22:39:58.183127Z",
     "iopub.status.busy": "2024-03-15T22:39:58.182814Z",
     "iopub.status.idle": "2024-03-15T22:39:58.347779Z",
     "shell.execute_reply": "2024-03-15T22:39:58.347334Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas\n",
    "\n",
    "dfTwu = pandas.read_csv('fitted_Twu_coeffs.csv')\n",
    "def get_model(INCHIKey):\n",
    "    row = dfTwu.loc[dfTwu['inchikey']==INCHIKey]\n",
    "    if len(row) == 1:\n",
    "        row = row.iloc[0]\n",
    "    Tc_K = row['Tc_K']\n",
    "    pc_MPa = row['pc_MPa']\n",
    "    c = [row['c0'], row['c1'], row['c2']]\n",
    "\n",
    "    # The JSON definition of the EOS,\n",
    "    # here a generic cubic EOS to allow for \n",
    "    # specification of the alpha function(s)\n",
    "    j = {\n",
    "        'kind': 'cubic',\n",
    "        'model': {\n",
    "            'type': 'PR',\n",
    "            'Tcrit / K': [Tc_K],\n",
    "            'pcrit / Pa': [pc_MPa*1e6],\n",
    "            'acentric': [0.1],\n",
    "            'alpha': [{'type': 'Twu', 'c': c}]\n",
    "        }\n",
    "    }\n",
    "    model = teqp.make_model(j)\n",
    "    return model, j \n",
    "\n",
    "# Hexane\n",
    "model, j = get_model(INCHIKey='VLKZOEOYAKHREP-UHFFFAOYSA-N')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2a5b8f8e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-15T22:39:58.350140Z",
     "iopub.status.busy": "2024-03-15T22:39:58.349840Z",
     "iopub.status.idle": "2024-03-15T22:39:58.354611Z",
     "shell.execute_reply": "2024-03-15T22:39:58.354195Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "102739.27983424198 6.03697343297877\n"
     ]
    }
   ],
   "source": [
    "# And how about we calculate the pressure and s^+ = -sr/R at NBP of water\n",
    "model, j = get_model(INCHIKey='XLYOFNOQVPJJNP-UHFFFAOYSA-N') # WATER\n",
    "\n",
    "T = 373.1242958476844 # K, NBP of water\n",
    "rhoL, rhoV = model.superanc_rhoLV(T)\n",
    "z = np.array([1.0])\n",
    "pL = rhoL*model.get_R(z)*T*(1.0 + model.get_Ar01(T, rhoL, z))\n",
    "splusL = model.get_splus(T, rhoL*z)\n",
    "print(pL, splusL)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b42d2266",
   "metadata": {},
   "source": [
    "Also implemented in version 0.17 are the alpha functions of Mathias-Copeman.\n",
    "\n",
    "$$\n",
    "\\alpha = (1+c_0x + c_1x^2 + c_2x^3)^2\n",
    "$$\n",
    "with\n",
    "$$\n",
    "x = 1 + \\sqrt{\\frac{T}{T_ci}}\n",
    "$$\n",
    "Parameters are tabulated for many fluids in the paper of Horstmann (https://doi.org/10.1016/j.fluid.2004.11.002) for the SRK EOS (only)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "7e2df55e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-15T22:39:58.356475Z",
     "iopub.status.busy": "2024-03-15T22:39:58.356172Z",
     "iopub.status.idle": "2024-03-15T22:39:58.360408Z",
     "shell.execute_reply": "2024-03-15T22:39:58.359952Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "And with SRK and Mathias-Copeman parameters: 101639.22259842217 Pa\n"
     ]
    }
   ],
   "source": [
    "# Here is an example from Horstmann\n",
    "j = {\n",
    "    \"kind\": \"cubic\",\n",
    "    \"model\": {\n",
    "        \"type\": \"SRK\",\n",
    "        \"Tcrit / K\": [647.30],\n",
    "        \"pcrit / Pa\": [22.048321e6],\n",
    "        \"acentric\": [0.3440],\n",
    "        \"alpha\": [\n",
    "            {\"type\": \"Mathias-Copeman\", \"c\": [1.07830, -0.58321, 0.54619]}\n",
    "        ]\n",
    "    }\n",
    "}\n",
    "\n",
    "model = teqp.make_model(j)\n",
    "T = 373.1242958476844 # K\n",
    "rhoL, rhoV = model.superanc_rhoLV(T)\n",
    "z = np.array([1.0])\n",
    "pL = rhoL*model.get_R(z)*T*(1.0 + model.get_Ar01(T, rhoL, z))\n",
    "print('And with SRK and Mathias-Copeman parameters:', pL, 'Pa')"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}