diff --git a/doc/source/models/RK-PR.ipynb b/doc/source/models/RK-PR.ipynb new file mode 100644 index 00000000..a3fc4e05 --- /dev/null +++ b/doc/source/models/RK-PR.ipynb @@ -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 +} diff --git a/include/teqp/models/cubics.hpp b/include/teqp/models/cubics.hpp index a8970f42..ccc3c24b 100644 --- a/include/teqp/models/cubics.hpp +++ b/include/teqp/models/cubics.hpp @@ -881,4 +881,120 @@ class QuantumCorrectedPR{ } }; + +/** + The RK-PR method as defined in: + + Martín Cismondi, Jørgen Mollerup, + Development and application of a three-parameter RK–PR equation of state, + Fluid Phase Equilibria, + Volume 232, Issues 1–2, + 2005, + Pages 74-89, + https://doi.org/10.1016/j.fluid.2005.03.020. + */ +class RKPRCismondi2005{ +public: + const double Ru = get_R_gas(); /// Universal gas constant, exact number + +private: + const std::vector delta_1, Tc_K, pc_Pa, k; + const Eigen::ArrayXXd kmat, lmat; + const std::vector a_c, b_c; + + /// A convenience function to save some typing + std::vector get_(const nlohmann::json &j, const std::string& k) const { return j.at(k).get>(); } + + /// Calculate the parameters \f$y\f$ and \f$d_1\f$ + auto get_yd1(double delta_1){ + return std::make_tuple(1 + cbrt(2*(1+delta_1)) + cbrt(4/(1+delta_1)), (1+delta_1*delta_1)/(1+delta_1)); + } + + auto build_ac(){ + std::vector a_c(delta_1.size()); + for (auto i = 0; i < delta_1.size(); ++i){ + auto [y, d_1] = get_yd1(delta_1[i]); + auto Omega_a = (3*y*y + 3*y*d_1 + d_1*d_1 + d_1 - 1.0)/pow(3.0*y + d_1 - 1.0, 2); + a_c[i] = Omega_a*pow(Ru*Tc_K[i], 2)/pc_Pa[i]; + } + return a_c; + } + auto build_bc(){ + std::vector b_c(delta_1.size()); + for (auto i = 0; i < delta_1.size(); ++i){ + auto [y, d_1] = get_yd1(delta_1[i]); + auto Omega_b = 1.0/(3.0*y + d_1 - 1.0); + b_c[i] = Omega_b*Ru*Tc_K[i]/pc_Pa[i]; + } + return b_c; + } +public: + + RKPRCismondi2005(const nlohmann::json &j) : delta_1(get_(j, "delta_1")), Tc_K(get_(j, "Tcrit / K")), pc_Pa(get_(j, "pcrit / Pa")), k(get_(j, "k")), kmat(build_square_matrix(j.at("kmat"))), lmat(build_square_matrix(j.at("lmat"))), a_c(build_ac()), b_c(build_bc()) {} + + template + auto R(const VecType& /*molefrac*/) const { + return Ru; + } + auto get_delta_1() const { return delta_1; } + auto get_Tc_K() const { return Tc_K; } + auto get_pc_Pa() const { return pc_Pa; } + auto get_k() const { return k; } + auto get_kmat() const { return kmat; } + auto get_lmat() const { return lmat; } + + template + auto get_bi(std::size_t i, const TType& T) const { + return b_c[i]; + } + + template + auto get_ai(std::size_t i, const TType& T) const { + return forceeval(a_c[i]*pow(3.0/(2.0+T/Tc_K[i]), k[i])); + } + + template + auto get_ab(const TType& T, const FractionsType& z) const{ + using numtype = std::common_type_t; + numtype b = 0.0; + numtype a = 0.0; + std::size_t N = delta_1.size(); + for (auto i = 0; i < N; ++i){ + auto bi = get_bi(i, T); + auto ai = get_ai(i, T); + for (auto j = 0; j < N; ++j){ + auto aj = get_ai(j, T); + auto bj = get_bi(j, T); + + a += z[i]*z[j]*sqrt(ai*aj)*(1.0 - kmat(i,j)); + b += z[i]*z[j]*(bi + bj)/2.0*(1.0 - lmat(i,j)); + + } + } + return std::make_tuple(a, b); + } + + template + auto alphar(const TType& T, const RhoType& rho, const FractionsType& molefrac) const { + if (molefrac.size() != delta_1.size()) { + throw std::invalid_argument("Sizes do not match"); + } + + // The mixture Delta_1 is a mole-fraction-weighted average of the pure values of delta_1 + const auto delta_1_view = Eigen::Map(&(delta_1[0]), delta_1.size()); + const decltype(molefrac[0]) Delta1 = (molefrac*delta_1_view).eval().sum(); + + auto Delta2 = (1.0-Delta1)/(1.0+Delta1); + auto [a, b] = get_ab(T, molefrac); + + auto Psiminus = -log(1.0 - b*rho); + auto Psiplus = log((Delta1 * b * rho + 1.0) / (Delta2 * b * rho + 1.0)) / (b * (Delta1 - Delta2)); + auto val = Psiminus - a/(Ru*T)*Psiplus; + return forceeval(val); + } +}; +using RKPRCismondi2005_t = decltype(RKPRCismondi2005({})); + + + }; // namespace teqp diff --git a/interface/CPP/teqp_impl_factory.cpp b/interface/CPP/teqp_impl_factory.cpp index 860f8e89..7c001acf 100644 --- a/interface/CPP/teqp_impl_factory.cpp +++ b/interface/CPP/teqp_impl_factory.cpp @@ -24,6 +24,7 @@ namespace teqp { {"cubic", [](const nlohmann::json& spec){ return make_owned(make_generalizedcubic(spec));}}, {"QCPRAasen", [](const nlohmann::json& spec){ return make_owned(QuantumCorrectedPR(spec));}}, {"advancedPRaEres", [](const nlohmann::json& spec){ return make_owned(make_AdvancedPRaEres(spec));}}, + {"RKPRCismondi2005", [](const nlohmann::json& spec){ return make_owned(RKPRCismondi2005(spec));}}, {"CPA", [](const nlohmann::json& spec){ return make_owned(CPA::CPAfactory(spec));}}, {"PCSAFT", [](const nlohmann::json& spec){ return make_owned(PCSAFT::PCSAFTfactory(spec));}}, diff --git a/interface/pybind11_wrapper.cpp b/interface/pybind11_wrapper.cpp index fcd0b907..8ed00880 100644 --- a/interface/pybind11_wrapper.cpp +++ b/interface/pybind11_wrapper.cpp @@ -110,6 +110,7 @@ const std::type_index EXP6_Kataoka1992_i{std::type_index(typeid(EXP6_Kataoka1992 const std::type_index twocenterLJF_i{std::type_index(typeid(twocenterLJF_t))}; const std::type_index QuantumPR_i{std::type_index(typeid(QuantumPR_t))}; const std::type_index advancedPRaEres_i{std::type_index(typeid(advancedPRaEres_t))}; +const std::type_index RKPRCismondi2005_i{std::type_index(typeid(RKPRCismondi2005_t))}; /** At runtime we can add additional model-specific methods that only apply for a particular model. We take in a Python-wrapped @@ -186,6 +187,18 @@ void attach_model_specific_methods(py::object& obj){ setattr("get_pc_Pa", MethodType(py::cpp_function([](py::object& o){ return get_typed(o).get_pc_Pa(); }), obj)); // setattr("get_meta", MethodType(py::cpp_function([](py::object& o){ return get_typed(o).get_meta(); }), obj)); // setattr("set_meta", MethodType(py::cpp_function([](py::object& o, const std::string& s){ return get_mutable_typed(o).set_meta(s); }, "self"_a, "s"_a), obj)); + } + else if (index == RKPRCismondi2005_i){ + setattr("get_ab", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed(o).get_ab(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj)); + setattr("get_delta_1", MethodType(py::cpp_function([](py::object& o){ return get_typed(o).get_delta_1(); }), obj)); + setattr("get_Tc_K", MethodType(py::cpp_function([](py::object& o){ return get_typed(o).get_Tc_K(); }), obj)); + setattr("get_pc_Pa", MethodType(py::cpp_function([](py::object& o){ return get_typed(o).get_pc_Pa(); }), obj)); + setattr("get_k", MethodType(py::cpp_function([](py::object& o){ return get_typed(o).get_k(); }), obj)); + + setattr("get_kmat", MethodType(py::cpp_function([](py::object& o){ return get_typed(o).get_kmat(); }), obj)); + setattr("get_lmat", MethodType(py::cpp_function([](py::object& o){ return get_typed(o).get_lmat(); }), obj)); +// setattr("get_meta", MethodType(py::cpp_function([](py::object& o){ return get_typed(o).get_meta(); }), obj)); +// setattr("set_meta", MethodType(py::cpp_function([](py::object& o, const std::string& s){ return get_mutable_typed(o).set_meta(s); }, "self"_a, "s"_a), obj)); } else if (index == AmmoniaWaterTillnerRoth_i){ setattr("TcNH3", get_typed(obj).TcNH3);