-
Notifications
You must be signed in to change notification settings - Fork 13
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Implement the RK-PR method of Cismondi
- Loading branch information
Showing
4 changed files
with
332 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,202 @@ | ||
| { | ||
| "cells": [ | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "e86a8de0", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "# RK-PR\n", | ||
| "\n", | ||
| "The EOS can be given as \n", | ||
| "\n", | ||
| "$$ \\alpha^{\\rm r} = \\psi^{(-)} - \\dfrac{a_m}{RT } \\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 EOS fixed constants of\n", | ||
| "$$\n", | ||
| "\\Delta_1 = \\sum_i x_i \\delta_{1,i}\n", | ||
| "$$\n", | ||
| "$$\n", | ||
| "\\Delta_2 = \\frac{1-\\Delta_1}{1+\\Delta_1}\n", | ||
| "$$\n", | ||
| "\n", | ||
| "The attractive term goes like\n", | ||
| "$$\n", | ||
| "a_{i} = a_{c,i}\\left(\\frac{2}{3+T/T_{c,i}}\\right)^{k_i}\n", | ||
| "$$\n", | ||
| "with quadratic mixing rules\n", | ||
| "$$\n", | ||
| "a_m = \\sum_i\\sum_jx_ix_j(1-k_{ij})\\sqrt{a_{i}(T)a_{j}(T)}\n", | ||
| "$$\n", | ||
| "And the covolume also gets quadratic mixing rules\n", | ||
| "$$\n", | ||
| "b_m = \\sum_i\\sum_jx_ix_j(1-l_{ij})(b_{i}+b_{j})/2\n", | ||
| "$$\n", | ||
| "\n", | ||
| "Thus, to implement the RK-PR model in predictive mode, the following steps are required:\n", | ||
| "\n", | ||
| "1. Obtain the critical parameters Tc, pc\n", | ||
| "2. Solve for delta_1 from the experimental critical compressibility factor, begin with the values from the correlation\n", | ||
| "3. Solve for k by fixing the pressure at the T=0.7*Tc. In the case (e.g, CO$_2$) that Tt < 0.7*Tc, use instead Tr=Tt/Tc \n", | ||
| "\n", | ||
| "It may be necessary to adjust the values of $\\delta_{1,i}$ and $k_i$ for an individual component to better match the behavior of more polar components." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "e1aa4062", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "import numpy as np \n", | ||
| "import scipy.optimize \n", | ||
| "import matplotlib.pyplot as plt \n", | ||
| "import teqp, numpy as np\n", | ||
| "import CoolProp.CoolProp as CP \n", | ||
| "import pandas\n", | ||
| "\n", | ||
| "def delta1_correlation(Zc):\n", | ||
| " # Eq. B.4 of Cismondi FPE 2005\n", | ||
| " d1 = 0.428363\n", | ||
| " d2 = 18.496215\n", | ||
| " d3 = 0.338426\n", | ||
| " d4 = 0.660000\n", | ||
| " d5 = 789.723105\n", | ||
| " d6 = 2.512392\n", | ||
| " return d1 + d2*(d3-Zc)**d4 + d5*(d3-Zc)**d6 \n", | ||
| "\n", | ||
| "def Zc_delta1(delta1):\n", | ||
| " # Eqs. B.1 to B.3 of Cismondi FPE 2005\n", | ||
| " d1 = (1+delta1**2)/(1+delta1)\n", | ||
| " y = 1 + (2*(1+delta1))**(1/3) + (4/(1+delta1))**(1/3)\n", | ||
| " return y/(3*y + d1 - 1)\n", | ||
| "\n", | ||
| "DELTA1 = np.linspace(np.sqrt(2)-1, 25, 1000)\n", | ||
| "ZZ = Zc_delta1(DELTA1)\n", | ||
| "plt.plot(DELTA1, ZZ, label='values')\n", | ||
| "DELTA1back = delta1_correlation(ZZ)\n", | ||
| "plt.axvline(np.sqrt(2)-1)\n", | ||
| "plt.plot(DELTA1back, ZZ, dashes=[2,2], label='correlation')\n", | ||
| "plt.gca().set(ylabel='$Z_c$', xlabel='$\\delta_1$')\n", | ||
| "plt.legend(loc='best')\n", | ||
| "plt.show()\n", | ||
| "\n", | ||
| "# for Zc in np.linspace(0.2, 0.3383, 1000):\n", | ||
| "# resid = lambda x: Zc_delta1(x)-Zc\n", | ||
| "# # print(resid(delta1_correlation(Zc)))\n", | ||
| "# print(Zc, scipy.optimize.newton(resid, delta1_correlation(Zc)), delta1_correlation(Zc))" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "8fb5d777", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "names = ['CO2', 'n-Decane']\n", | ||
| "\n", | ||
| "R = 8.31446261815324\n", | ||
| "Tc = np.array([CP.PropsSI(k,\"Tcrit\") for k in names])\n", | ||
| "pc = np.array([CP.PropsSI(k,\"pcrit\") for k in names])\n", | ||
| "rhoc = np.array([CP.PropsSI(k,\"rhomolar_critical\") for k in names])\n", | ||
| "Zcexp = pc/(rhoc*R*Tc)\n", | ||
| "\n", | ||
| "# Use a rescaled Zc to obtain delta_1\n", | ||
| "Zc = 1.168*Zcexp\n", | ||
| "\n", | ||
| "delta_1 = [scipy.optimize.newton(lambda x: Zc_delta1(x)-Zc_, delta1_correlation(Zc_)) for Zc_ in Zc]\n", | ||
| "\n", | ||
| "def solve_for_k(i, p_target, Tr):\n", | ||
| " \"\"\" \n", | ||
| " The value of k for the i-th component is based on getting \n", | ||
| " the right vapor pressure, so a rootfinding routing is\n", | ||
| " used to obtain these values\n", | ||
| " \"\"\"\n", | ||
| " def objective(k):\n", | ||
| " j = {\n", | ||
| " \"kind\": \"RKPRCismondi2005\",\n", | ||
| " \"model\": {\n", | ||
| " \"delta_1\": [delta_1[i]],\n", | ||
| " \"Tcrit / K\": [Tc[i]],\n", | ||
| " \"pcrit / Pa\": [pc[i]],\n", | ||
| " \"k\": [k],\n", | ||
| " \"kmat\": [[0.0]],\n", | ||
| " \"lmat\": [[0.0]],\n", | ||
| " }\n", | ||
| " }\n", | ||
| " model = teqp.make_model(j)\n", | ||
| " T = Tr*Tc[i]\n", | ||
| " z = np.array([1.0])\n", | ||
| " a, b = model.get_ab(T, z)\n", | ||
| "\n", | ||
| " anc = teqp.build_ancillaries(model, Tc[i], rhoc[i], 150)\n", | ||
| " rhoL, rhoV = model.pure_VLE_T(T, anc.rhoL(T), anc.rhoV(T), 10)\n", | ||
| " p = T*R*rhoL*(1+model.get_Ar01(T, rhoL, z))\n", | ||
| " \n", | ||
| " return p-p_target\n", | ||
| " return scipy.optimize.newton(objective, 2.1)\n", | ||
| "\n", | ||
| "Tr = 0.7\n", | ||
| "i = 1\n", | ||
| "k_C10 = solve_for_k(i, CP.PropsSI('P','T',Tr*Tc[i],'Q',0,names[i]), Tr)\n", | ||
| "\n", | ||
| "model = teqp.make_model({\n", | ||
| " \"kind\": \"RKPRCismondi2005\",\n", | ||
| " \"model\": {\n", | ||
| " \"delta_1\": delta_1,\n", | ||
| " \"Tcrit / K\": Tc.tolist(),\n", | ||
| " \"pcrit / Pa\": pc.tolist(),\n", | ||
| " \"k\": [2.23854, k_C10],\n", | ||
| " \"kmat\": [[0,0],[0,0]],\n", | ||
| " \"lmat\": [[0,0],[0,0]],\n", | ||
| " }\n", | ||
| "})\n", | ||
| "\n", | ||
| "# Start at both pures\n", | ||
| "for ipure in [0, 1]:\n", | ||
| " Tc, rhoc = model.solve_pure_critical(300, 5000, {\"alternative_pure_index\":ipure, \"alternative_length\": 2})\n", | ||
| " z = np.array([0.0, 0.0]); z[ipure] = 1.0\n", | ||
| " pc = Tc*R*rhoc*(1+model.get_Ar01(Tc, rhoc, z))\n", | ||
| " plt.plot(Tc, pc/1e5, 'o')\n", | ||
| " \n", | ||
| " opt = teqp.TCABOptions(); opt.polish=True; opt.verbosity=100; opt.integration_order=5; opt.rel_err=1e-10; opt.abs_err=1e-10\n", | ||
| " trace = model.trace_critical_arclength_binary(Tc, z*rhoc, options=opt)\n", | ||
| " df = pandas.DataFrame(trace)\n", | ||
| " plt.plot(df['T / K'], df['p / Pa']/1e5)\n", | ||
| " \n", | ||
| "# Overlay the data from Reamer and Sage, Cismondi additional data points not present in Reamer and Sage\n", | ||
| "Tc_K = [310.928, 344.261, 377.594, 410.928, 444.261, 477.594, 510.928]\n", | ||
| "pc_kPa = np.array([7997.92, 12824.25, 16492.26, 18560.69, 18836.48, 17836.74, 15333.94])\n", | ||
| "plt.plot(Tc_K, pc_kPa/1e2, 'o')\n", | ||
| "\n", | ||
| "plt.gca().set(xlabel='$T$ / K', ylabel='$p$ / bar');" | ||
| ] | ||
| } | ||
| ], | ||
| "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 | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters