From e7fbb79f947cb2e3576b510a21aa2db699b15911 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Mon, 16 Oct 2023 09:50:41 -0400 Subject: [PATCH 1/3] Add advanced cubic mixing rules from the work of Lasala Is working! --- include/teqp/models/cubics.hpp | 277 ++++++++++++++++++++++++++++++++ src/tests/catch_test_cubics.cxx | 49 ++++++ 2 files changed, 326 insertions(+) diff --git a/include/teqp/models/cubics.hpp b/include/teqp/models/cubics.hpp index 621ef357..9f666e3d 100644 --- a/include/teqp/models/cubics.hpp +++ b/include/teqp/models/cubics.hpp @@ -400,6 +400,283 @@ inline auto make_generalizedcubic(const nlohmann::json& spec){ return cub; } +enum class AdvancedPRaEMixingRules {knotspecified, kLinear, kQuadratic}; + +struct AdvancedPRaEOptions{ + AdvancedPRaEMixingRules brule = AdvancedPRaEMixingRules::kQuadratic; + double s = 2.0; + double CEoS = -sqrt(2.0)/2.0*log(1.0 + sqrt(2.0)); +}; + + +/** + A residual Helmholtz term that returns nothing (the empty term) + */ +template +class NullResidualHelmholtzOverRT { +public: + template + auto operator () (const TType& T, const MoleFractions& molefracs) const { + std::common_type_t val = 0.0; + return val; + } +}; + +/** + + \f[ + \frac{a^{E,\gamma}_{total}}{RT} = -sum_iz_i\ln\left(\sum_jz_jOmega_{ij}(T)\right) + \f] + + \f[ + \frac{a^{E,\gamma}_{comb}}{RT} = -sum_iz_i\ln\left(\frac{\Omega_i}{z_i}\right) + \f] + + \f[ + \frac{a^{E,\gamma}_{res}}{RT} = \frac{a^{E,\gamma}_{total}}{RT} - \frac{a^{E,\gamma}_{comb}}{RT} + \f] + + Volume fraction of component \f$i\f$ +\f[ + \phi_i = \frac{z_iv_i}{\sum_j z_j v_j} + \f] + with \f$v_i = b_i\f$ + */ +template +class WilsonResidualHelmholtzOverRT { + +public: + const double R = 8.31446261815324; + const std::vector b; + const Eigen::ArrayXXd m, n; + WilsonResidualHelmholtzOverRT(const std::vector& b, const Eigen::ArrayXXd& m, const Eigen::ArrayXXd& n) : b(b), m(m), n(n) {}; + + template + auto combinatorial(const TType& T, const MoleFractions& molefracs) const { + if (b.size() != molefracs.size()){ + throw teqp::InvalidArgument("Bad size of molefracs"); + } + + using TYPE = std::common_type_t; + // The denominator in Phi + TYPE Vtot = 0.0; + for (auto i = 0; i < molefracs.size(); ++i){ + auto v_i = b[i]; + Vtot += molefracs[i]*v_i; + } + + TYPE summer = 0.0; + for (auto i = 0; i < molefracs.size(); ++i){ + auto v_i = b[i]; + // The ratio phi_i/z_i is expressed like this to better handle + // the case of z_i = 0, which would otherwise be a divide by zero + // in the case that the composition of one component is zero + auto phi_i_over_z_i = v_i/Vtot; + summer += molefracs[i]*log(phi_i_over_z_i); + } + return summer; + }; + + template + auto get_Aij(std::size_t i, std::size_t j, const TType& T) const{ + return forceeval(m(i,j)*T + n(i,j)); + } + + template + auto total(const TType& T, const MoleFractions& molefracs) const { + + using TYPE = std::common_type_t; + TYPE summer = 0.0; + for (auto i = 0; i < molefracs.size(); ++i){ + auto v_i = b[i]; + TYPE summerj = 0.0; + for (auto j = 0; j < molefracs.size(); ++j){ + auto v_j = b[j]; + auto Aij = get_Aij(i,j,T); + auto Omega_ji = v_j/v_i*exp(-Aij/T); + summerj += molefracs[j]*Omega_ji; + } + summer += molefracs[i]*log(summerj); + } + return forceeval(-summer); + }; + + // Returns ares/RT + template + auto operator () (const TType& T, const MoleFractions& molefracs) const { + return forceeval(total(T, molefracs) - combinatorial(T, molefracs)); + } +}; + +using ResidualHelmholtzOverRTOptions = std::variant, WilsonResidualHelmholtzOverRT>; + +/** + Cubic EOS with advanced mixing rules, the EoS/aE method of Jaubert and co-workers + + */ +template > +class AdvancedPRaEres { +public: + // Hard-coded values for Peng-Robinson + const NumType Delta1 = 1+sqrt(2.0); + const NumType Delta2 = 1-sqrt(2.0); + // See https://doi.org/10.1021/acs.iecr.1c00847 + const NumType OmegaA = 0.45723552892138218938; + const NumType OmegaB = 0.077796073903888455972; + const int superanc_code = CubicSuperAncillary::PR_CODE; + const double CEoS; + +protected: + + std::valarray ai, bi; + + const AlphaFunctions alphas; + const ResidualHelmholtzOverRTOptions ares; + Eigen::ArrayXXd lmat; + const double s; + const AdvancedPRaEMixingRules brule; + + nlohmann::json meta; + + template + auto get_ai(TType& T, IndexType i) const { + auto alphai = std::visit([&](auto& t) { return t(T); }, alphas[i]); + return forceeval(ai[i]*alphai); + } + + template + auto get_bi(TType& /*T*/, IndexType i) const { return bi[i]; } + + template + void check_lmat(IndexType N) { + if (lmat.cols() != lmat.rows()) { + throw teqp::InvalidArgument("lmat rows [" + std::to_string(lmat.rows()) + "] and columns [" + std::to_string(lmat.cols()) + "] are not identical"); + } + if (lmat.cols() == 0) { + lmat.resize(N, N); lmat.setZero(); + } + else if (lmat.cols() != N) { + throw teqp::InvalidArgument("lmat needs to be a square matrix the same size as the number of components [" + std::to_string(N) + "]"); + } + }; + +public: + AdvancedPRaEres(const std::valarray& Tc_K, const std::valarray& pc_Pa, const AlphaFunctions& alphas, const ResidualHelmholtzOverRTOptions& ares, const Eigen::ArrayXXd& lmat, const AdvancedPRaEOptions& options = {}) + : alphas(alphas), ares(ares), lmat(lmat), s(options.s), brule(options.brule), CEoS(options.CEoS) + { + ai.resize(Tc_K.size()); + bi.resize(Tc_K.size()); + for (auto i = 0; i < Tc_K.size(); ++i) { + ai[i] = OmegaA * pow2(Ru * Tc_K[i]) / pc_Pa[i]; + bi[i] = OmegaB * Ru * Tc_K[i] / pc_Pa[i]; + } + check_lmat(ai.size()); + }; + + void set_meta(const nlohmann::json& j) { meta = j; } + auto get_meta() const { return meta; } + auto get_lmat() const { return lmat; } + static double get_bi(double Tc_K, double pc_Pa){ + const NumType OmegaB = 0.077796073903888455972; + const NumType R = 8.31446261815324; + return OmegaB*R*Tc_K/pc_Pa; + } + + /// Return a tuple of saturated liquid and vapor densities for the EOS given the temperature + /// Uses the superancillary equations from Bell and Deiters: + /// \param T Temperature + /// \param ifluid Must be provided in the case of mixtures + auto superanc_rhoLV(double T, std::optional ifluid = std::nullopt) const { + + std::valarray molefracs(ai.size()); molefracs = 1.0; + + // If more than one component, must provide the ifluid argument + if(ai.size() > 1){ + if (!ifluid){ + throw teqp::InvalidArgument("For mixtures, the argument ifluid must be provided"); + } + if (ifluid.value() > ai.size()-1){ + throw teqp::InvalidArgument("ifluid must be less than "+std::to_string(ai.size())); + } + molefracs = 0.0; + molefracs[ifluid.value()] = 1.0; + } + + auto b = get_b(T, molefracs); + auto a = get_am_over_bm(T, molefracs)*b; + auto Ttilde = R(molefracs)*T*b/a; + return std::make_tuple( + CubicSuperAncillary::supercubic(superanc_code, CubicSuperAncillary::RHOL_CODE, Ttilde)/b, + CubicSuperAncillary::supercubic(superanc_code, CubicSuperAncillary::RHOV_CODE, Ttilde)/b + ); + } + + const NumType Ru = get_R_gas(); /// Universal gas constant, exact number + + template + auto R(const VecType& /*molefrac*/) const { + return Ru; + } + + template + auto get_a(TType T, const CompType& molefracs) const { + return forceeval(get_am_over_bm(T, molefracs)*get_b(T, molefracs)); + } + + template + auto get_am_over_bm(TType T, const CompType& molefracs) const { + auto aEresRT = std::visit([&](auto& aresRTfunc) { return aresRTfunc(T, molefracs); }, ares); // aEres/RT, so a non-dimensional quantity + std::common_type_t summer = aEresRT*Ru*T/CEoS; + for (auto i = 0; i < molefracs.size(); ++i) { + summer += molefracs[i]*get_ai(T,i)/get_bi(T,i); + } + return forceeval(summer); + } + + template + auto get_b(TType T, const CompType& molefracs) const { + std::common_type_t b_ = 0.0; + + switch (brule){ + case AdvancedPRaEMixingRules::kQuadratic: + for (auto i = 0; i < molefracs.size(); ++i) { + auto bi_ = get_bi(T, i); + for (auto j = 0; j < molefracs.size(); ++j) { + auto bj_ = get_bi(T, j); + + auto bij = (1 - lmat(i,j)) * pow((pow(bi_, 1.0/s) + pow(bj_, 1.0/s))/2.0, s); + b_ += molefracs[i] * molefracs[j] * bij; + } + } + break; + case AdvancedPRaEMixingRules::kLinear: + for (auto i = 0; i < molefracs.size(); ++i) { + b_ += molefracs[i] * get_bi(T, i); + } + break; + default: + throw teqp::InvalidArgument("Mixing rule for b is invalid"); + } + return forceeval(b_); + } + + template + auto alphar(const TType& T, + const RhoType& rho, + const MoleFracType& molefrac) const + { + if (molefrac.size() != alphas.size()) { + throw std::invalid_argument("Sizes do not match"); + } + auto b = get_b(T, molefrac); + auto a = get_am_over_bm(T, molefrac)*b; + 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); + } +}; + /** The quantum corrected Peng-Robinson model as developed in diff --git a/src/tests/catch_test_cubics.cxx b/src/tests/catch_test_cubics.cxx index 71a40005..34c5fdf6 100644 --- a/src/tests/catch_test_cubics.cxx +++ b/src/tests/catch_test_cubics.cxx @@ -628,3 +628,52 @@ TEST_CASE("QCPR", "[QCPR]"){ auto z = (Eigen::ArrayXd(2) << 0.3, 0.7).finished(); CHECK(std::isfinite(model->get_B12vir(T, z))); } + +TEST_CASE("Advanced cubic EOS", "[AdvancedPR]"){ + + // Values for CO2 + N2 from Table 6.3 of Lasala dissertation which come from DIPPR + std::valarray Tc_K = {304.21, 126.19}, + pc_Pa = {7.383e6, 3395800.0}, + acentric = {0.22394, 0.0372}; + + // Classic Peng-Robinson alpha function and covolume + std::vector alphas; + std::vector b; + for (auto i = 0; i < acentric.size(); ++i){ + auto mi = 0.37464 + 1.54226*acentric[i] - 0.26992*acentric[i]*acentric[i]; + alphas.push_back(BasicAlphaFunction(Tc_K[i], mi)); + b.push_back(teqp::AdvancedPRaEres::get_bi(Tc_K[i], pc_Pa[i])); + } + + SECTION("Check pcrit"){ + + // Matrices for putting the coefficients in directly +// Eigen::ArrayXXd mWilson = (Eigen::ArrayXXd(2,2) << 0.0, -3.4768, 3.5332, 0.0).finished().transpose(); +// Eigen::ArrayXXd nWilson = (Eigen::ArrayXXd(2,2) << 0.0, 825, -585, 0.0).finished().transpose(); + + Eigen::ArrayXXd mWilson = (Eigen::ArrayXXd(2,2) << 0.0, 0.0, 0.0, 0.0).finished(); + + for (double T: {223.1, 253.05, 273.1}){ + std::size_t ipure = 0; + + double A12 = -3.4768*T + 825; + double A21 = 3.5332*T - 585; + Eigen::ArrayXXd nWilson = (Eigen::ArrayXXd(2,2) << 0.0, A12, A21, 0.0).finished(); + auto aresmodel = WilsonResidualHelmholtzOverRT(b, mWilson, nWilson); + AdvancedPRaEOptions options; options.CEoS = -0.52398; + auto model = teqp::AdvancedPRaEres(Tc_K, pc_Pa, alphas, aresmodel, Eigen::ArrayXXd::Zero(2, 2), options); + + // Solve for starting point with superancillary function + auto [rhoL, rhoV] = model.superanc_rhoLV(T, ipure); + + Eigen::ArrayXd rhovecL0 = Eigen::ArrayXd::Zero(2); rhovecL0(ipure) = rhoL; + Eigen::ArrayXd rhovecV0 = Eigen::ArrayXd::Zero(2); rhovecV0(ipure) = rhoV; + + TVLEOptions opt; opt.revision = 2; + auto J = trace_VLE_isotherm_binary(model, T, rhovecL0, rhovecV0, opt); + auto pcrit = J["data"].back()["pL / Pa"]; + std::cout << T << ", " << pcrit << std::endl; + std::ofstream file("isoT_advcub"+std::to_string(T)+".json"); file << J; + } + } +} From 08d751cd800371cf22a0b9507421abdbf9d3b58f Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Mon, 16 Oct 2023 09:52:06 -0400 Subject: [PATCH 2/3] Also update the superancillary function for normal cubics --- include/teqp/models/cubics.hpp | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/include/teqp/models/cubics.hpp b/include/teqp/models/cubics.hpp index 9f666e3d..31940341 100644 --- a/include/teqp/models/cubics.hpp +++ b/include/teqp/models/cubics.hpp @@ -142,14 +142,28 @@ class GenericCubic { auto get_kmat() const { return kmat; } /// Return a tuple of saturated liquid and vapor densities for the EOS given the temperature - /// Uses the superancillary equations from Bell and Deiters: - auto superanc_rhoLV(double T) const { - if (ai.size() != 1) { - throw std::invalid_argument("function only available for pure species"); + /// Uses the superancillary equations from Bell and Deiters: + /// \param T Temperature + /// \param ifluid Must be provided in the case of mixtures + auto superanc_rhoLV(double T, std::optional ifluid = std::nullopt) const { + + std::valarray molefracs(ai.size()); molefracs = 1.0; + + // If more than one component, must provide the ifluid argument + if(ai.size() > 1){ + if (!ifluid){ + throw teqp::InvalidArgument("For mixtures, the argument ifluid must be provided"); + } + if (ifluid.value() > ai.size()-1){ + throw teqp::InvalidArgument("ifluid must be less than "+std::to_string(ai.size())); + } + molefracs = 0.0; + molefracs[ifluid.value()] = 1.0; } - const std::valarray z = { 1.0 }; - auto b = get_b(T, z); - auto Ttilde = R(z)*T*b/get_a(T,z); + + auto b = get_b(T, molefracs); + auto a = get_a(T, molefracs); + auto Ttilde = R(molefracs)*T*b/a; return std::make_tuple( CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOL_CODE, Ttilde)/b, CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOV_CODE, Ttilde)/b From 1dfaa274ef12c707e032a426df6f29b7d8a1275e Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Mon, 16 Oct 2023 11:56:37 -0400 Subject: [PATCH 3/3] Finalized python interfacing and example --- doc/source/models/AdvancedCubicMixing.ipynb | 179 ++++++++++++++++++++ doc/source/models/index.rst | 1 + include/teqp/models/cubics.hpp | 126 ++++++++++---- interface/CPP/teqp_impl_factory.cpp | 1 + interface/pybind11_wrapper.cpp | 10 ++ src/tests/catch_test_cubics.cxx | 14 ++ 6 files changed, 301 insertions(+), 30 deletions(-) create mode 100644 doc/source/models/AdvancedCubicMixing.ipynb diff --git a/doc/source/models/AdvancedCubicMixing.ipynb b/doc/source/models/AdvancedCubicMixing.ipynb new file mode 100644 index 00000000..d282ac9e --- /dev/null +++ b/doc/source/models/AdvancedCubicMixing.ipynb @@ -0,0 +1,179 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "92b3c17a", + "metadata": {}, + "source": [ + "# Advanced cubic mixing rules\n", + "\n", + "\n", + "In the advanced cubic mixing rules, the conventional cubic EOS is taken as the basis for the method (usually Peng-Robinson), but different rules are used for the attractive term $a_m$. The formulation reads:\n", + "\n", + "$$\n", + "\\frac{a_m}{b_m} = \\sum_i z_i\\frac{a_i}{b_i} + \\frac{a^{E,\\gamma}_{\\rm res}}{CEoS}\n", + "$$\n", + "\n", + "where $CEoS$ is a scaling parameter that is in principle linked with the EOS coefficients, but can also be allowed to be an adjustable parameter. The $a_i$ and $b_i$ are the pure fluid values of component $i$. The $z_i$ are mole fractions. The mixture covolume is given by \n", + "$$\n", + "b_m = \\sum_i\\sum_j z_iz_jb_{ij}\n", + "$$\n", + "with \n", + "$$\n", + "b_{ij} = \\left(\\frac{b_i^{1/s} + b_j^{1/s}}{2}\\right)^s\n", + "$$\n", + "\n", + "The heart of the method is the definition of $a^{E,\\gamma}_{\\rm res}$, the residual contribution (not in the conventional thermodynamic sense) to the excess Helmholtz energy. There are many possible models here, but one that seems to work well is that of Wilson, for which the expression reads:\n", + "\n", + "$$\n", + "\\frac{a^{E,\\gamma}_{\\rm res}}{RT} = -\\sum_i z_i\\ln\\left(\\sum_j z_j \\Omega_{ji}(T)\\right) - \\sum_iz_i\\ln\\left(\\frac{\\phi_i}{z_i}\\right)\n", + "$$\n", + "with \n", + "$$\n", + "\\Omega_{ji}=\\frac{v_j}{v_i}\\exp(-A_{ij}/T)\n", + "$$\n", + "and \n", + "$$\n", + "\\frac{\\phi_i}{z_i} = \\frac{v_i}{\\sum_kz_kv_k}\n", + "$$\n", + "with the $v_i=b_i$. The parameter $A_{ij}\\neq A_{ji}$ in general, and is also given temperature dependence, which is also not supposed to be present according to the derivation. Thus, the models for $A_{ij}$ read something like this here:\n", + "$$\n", + "A_{ij} = m_{ij}T + n_{ij}\n", + "$$\n", + "so $m$ is non-dimensional and $n$ has units of temperature." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "de8cc3c3", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy, matplotlib.pyplot as plt, numpy as np, pandas\n", + "import teqp\n", + "teqp.__version__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7e7e2f96", + "metadata": {}, + "outputs": [], + "source": [ + "# Four isotherms of experimental data from doi: 10.1016/j.fluid.2016.05.015\n", + "import io, pandas\n", + "dat = pandas.read_csv(io.StringIO(\"\"\"PointID y1 uy1 x1 ux1 p/bar up T/K\n", + "1 0.0274 0.0007 0.0068 0.0002 59.830 0.053 293.10 \n", + "2 0.0664 0.0014 0.0183 0.0004 64.864 0.080 293.10 \n", + "3 0.0978 0.0020 0.0298 0.0007 69.772 0.080 293.10 \n", + "4 0.1199 0.0024 0.0424 0.0009 74.737 0.080 293.10 \n", + "5 0.1219 0.0028 0.1132 0.0023 89.869 0.080 293.10 \n", + "6 0.1339 0.0024 0.0995 0.0022 89.198 0.080 293.10 \n", + "7 0.1399 0.0026 0.0943 0.0020 88.853 0.080 293.10 \n", + "8 0.1461 0.0027 0.0823 0.0019 86.962 0.080 293.10 \n", + "9 0.1466 0.0028 0.0778 0.0017 85.942 0.080 293.10 \n", + "10 0.1466 0.0028 0.0772 0.0016 85.868 0.080 293.10 \n", + "1 0.1378 0.0027 0.0159 0.0004 42.667 0.051 273.08 \n", + "2 0.2143 0.0038 0.0297 0.0007 49.547 0.051 273.08 \n", + "3 0.2612 0.0043 0.0411 0.0009 55.238 0.051 273.08 \n", + "4 0.3209 0.0049 0.0609 0.0013 65.069 0.088 273.08 \n", + "5 0.3554 0.0051 0.0786 0.0016 73.395 0.088 273.08 \n", + "6 0.3758 0.0052 0.0978 0.0019 81.061 0.088 273.08 \n", + "7 0.3903 0.0053 0.1190 0.0023 90.706 0.088 273.08 \n", + "8 0.3914 0.0053 0.1477 0.0028 100.966 0.088 273.08 \n", + "9 0.3879 0.0053 0.1614 0.0030 104.806 0.088 273.08 \n", + "10 0.3724 0.0052 0.1875 0.0033 110.846 0.088 273.08\n", + "11 0.3550 0.0051 0.2068 0.0036 114.105 0.088 273.08\n", + "12 0.2727 0.0044 0.2531 0.0041 118.020 0.088 273.08\n", + "13 0.3343 0.0049 0.2268 0.0038 116.295 0.088 273.08\n", + "1 0.2048 0.0038 0.0106 0.0003 25.754 0.050 253.05 \n", + "2 0.3019 0.0049 0.0217 0.0005 30.479 0.050 253.05 \n", + "3 0.4638 0.0056 0.0436 0.0010 45.352 0.050 253.05 \n", + "4 0.5319 0.0056 0.0647 0.0014 58.188 0.050 253.05 \n", + "5 0.5854 0.0054 0.1077 0.0021 78.315 0.084 253.05 \n", + "6 0.5979 0.0054 0.1497 0.0028 98.276 0.084 253.05\n", + "7 0.5898 0.0054 0.1801 0.0032 109.241 0.084 253.05\n", + "8 0.5042 0.0057 0.0570 0.0012 51.343 0.084 253.05\n", + "9 0.5644 0.0055 0.0861 0.0017 67.594 0.084 253.05\n", + "10 0.5949 0.0054 0.1267 0.0024 86.883 0.084 253.05\n", + "11 0.5826 0.0054 0.2015 0.0035 116.614 0.084 253.05\n", + "12 0.5537 0.0055 0.2431 0.0040 129.873 0.084 253.05\n", + "13 0.4973 0.0055 0.2971 0.0046 139.161 0.084 253.05\n", + "14 0.4971 0.0055 0.2972 0.0046 139.261 0.084 253.05\n", + "1 0.7076 0.0050 0.0257 0.0006 27.983 0.056 223.10\n", + "2 0.7774 0.0041 0.0522 0.0011 44.918 0.056 223.10\n", + "3 0.8077 0.0036 0.0930 0.0019 64.906 0.081 223.10\n", + "4 0.8131 0.0035 0.1261 0.0024 84.799 0.081 223.10\n", + "5 0.8057 0.0035 0.1584 0.0029 104.410 0.081 223.10\n", + "6 0.7843 0.0038 0.1982 0.0035 125.782 0.081 223.10\n", + "7 0.7533 0.0041 0.2380 0.0040 144.287 0.081 223.10\n", + "8 0.7150 0.0045 0.2813 0.0044 159.015 0.081 223.10\n", + "9 0.6942 0.0047 0.3064 0.0047 165.347 0.081 223.10\n", + "\"\"\"), sep='\\s+', engine='python')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "49c01826", + "metadata": {}, + "outputs": [], + "source": [ + "# Model from Lasala, FPE, 2016: https://doi.org/10.1016/j.fluid.2016.05.015\n", + "j = {\n", + " \"kind\": \"advancedPRaEres\",\n", + " \"model\": {\n", + " \"Tcrit / K\": [304.21, 126.19],\n", + " \"pcrit / Pa\": [7.383e6, 3395800.0],\n", + " \"alphas\": [{\"type\": \"PR78\", \"acentric\": 0.22394}, {\"type\": \"PR78\", \"acentric\": 0.0372}],\n", + " \"aresmodel\": {\"type\": \"Wilson\", \"m\": [[0.0, -3.4768], [3.5332, 0.0]], \"n\": [[0.0, 825], [-585, 0.0]]},\n", + " \"options\": {\"s\": 2.0, \"brule\": \"Quadratic\", \"CEoS\": -0.52398}\n", + " }\n", + "}\n", + "\n", + "model = teqp.make_model(j)\n", + "for T in [223.15, 253.05, 273.08, 293.1]:\n", + " ipure = 0\n", + "\n", + " [rhoL0, rhoV0] = model.superanc_rhoLV(T, ipure)\n", + "\n", + " rhovecL0 = np.array([0.0, 0.0]); rhovecL0[ipure] = rhoL0\n", + " rhovecV0 = np.array([0.0, 0.0]); rhovecV0[ipure] = rhoV0\n", + "\n", + " J = model.trace_VLE_isotherm_binary(T, rhovecL0, rhovecV0, opt)\n", + " df = pandas.DataFrame(J)\n", + " plt.plot(df['xL_0 / mole frac.'], df['pL / Pa']/1e6,'k')\n", + " plt.plot(df['xV_0 / mole frac.'], df['pV / Pa']/1e6,'k')\n", + "\n", + "plt.plot(1-dat['x1'], dat['p/bar']/10, 'o')\n", + "plt.plot(1-dat['y1'], dat['p/bar']/10, '^')\n", + "\n", + "plt.gca().set(xlabel='$x_1$ / mole frac.', ylabel='$p$ / MPa', ylim=(0, 25))\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 +} diff --git a/doc/source/models/index.rst b/doc/source/models/index.rst index 56450765..6d8314d5 100644 --- a/doc/source/models/index.rst +++ b/doc/source/models/index.rst @@ -8,6 +8,7 @@ Models make_model cubics QuantumPR + AdvancedCubicMixing model_potentials CPA PCSAFT diff --git a/include/teqp/models/cubics.hpp b/include/teqp/models/cubics.hpp index 31940341..33555880 100644 --- a/include/teqp/models/cubics.hpp +++ b/include/teqp/models/cubics.hpp @@ -94,6 +94,44 @@ class MathiasCopemanAlphaFunction { using AlphaFunctionOptions = std::variant, TwuAlphaFunction, MathiasCopemanAlphaFunction>; +template +auto build_alpha_functions(const TC& Tc_K, const nlohmann::json& jalphas){ + std::vector alphas; + std::size_t i = 0; + if (jalphas.size() != Tc_K.size()){ + throw teqp::InvalidArgument("alpha must be the same length as components"); + } + for (auto alpha : jalphas){ + std::string type = alpha.at("type"); + if (type == "Twu"){ + std::valarray c_ = alpha.at("c"); + Eigen::Array3d c = Eigen::Map(&(c_[0]), c_.size()); + alphas.emplace_back(TwuAlphaFunction(Tc_K[i], c)); + } + else if (type == "Mathias-Copeman"){ + std::valarray c_ = alpha.at("c"); + Eigen::Array3d c = Eigen::Map(&(c_[0]), c_.size()); + alphas.emplace_back(MathiasCopemanAlphaFunction(Tc_K[i], c)); + } + else if (type == "PR78"){ + double acentric = alpha.at("acentric"); + double mi = 0.0; + if (acentric < 0.491) { + mi = 0.37464 + 1.54226*acentric - 0.26992*pow2(acentric); + } + else { + mi = 0.379642 + 1.48503*acentric -0.164423*pow2(acentric) + 0.016666*pow3(acentric); + } + alphas.emplace_back( BasicAlphaFunction(Tc_K[i], mi)); + } + else{ + throw teqp::InvalidArgument("alpha type is not understood: "+type); + } + i++; + } + return alphas; +}; + template class GenericCubic { protected: @@ -321,28 +359,6 @@ inline auto make_generalizedcubic(const nlohmann::json& spec){ double Delta1, Delta2, OmegaA, OmegaB; std::string kind = "custom"; - auto add_alphas = [&](const nlohmann::json& jalphas){ - std::size_t i = 0; - if (jalphas.size() != Tc_K.size()){ - throw teqp::InvalidArgument("alpha must be the same length as components"); - } - for (auto alpha : jalphas){ - std::string type = alpha.at("type"); - std::valarray c_ = alpha.at("c"); - Eigen::Array3d c = Eigen::Map(&(c_[0]), c_.size()); - if (type == "Twu"){ - alphas.emplace_back(TwuAlphaFunction(Tc_K[i], c)); - } - else if (type == "Mathias-Copeman"){ - alphas.emplace_back(MathiasCopemanAlphaFunction(Tc_K[i], c)); - } - else{ - throw teqp::InvalidArgument("alpha type is not understood: "+type); - } - i++; - } - }; - if (spec.at("type") == "PR" ){ Delta1 = 1+sqrt(2.0); Delta2 = 1-sqrt(2.0); @@ -368,7 +384,7 @@ inline auto make_generalizedcubic(const nlohmann::json& spec){ if (!spec["alpha"].is_array()){ throw teqp::InvalidArgument("alpha must be array of objects"); } - add_alphas(spec.at("alpha")); + alphas = build_alpha_functions(Tc_K, spec.at("alpha")); } } else if (spec.at("type") == "SRK"){ @@ -384,7 +400,7 @@ inline auto make_generalizedcubic(const nlohmann::json& spec){ if (!spec["alpha"].is_array()){ throw teqp::InvalidArgument("alpha must be array of objects"); } - add_alphas(spec.at("alpha")); + alphas = build_alpha_functions(Tc_K, spec.at("alpha")); } // See https://doi.org/10.1021/acs.iecr.1c00847 OmegaA = 1.0 / (9.0 * (cbrt(2) - 1)); @@ -416,12 +432,23 @@ inline auto make_generalizedcubic(const nlohmann::json& spec){ enum class AdvancedPRaEMixingRules {knotspecified, kLinear, kQuadratic}; +NLOHMANN_JSON_SERIALIZE_ENUM( AdvancedPRaEMixingRules, { + {AdvancedPRaEMixingRules::knotspecified, nullptr}, + {AdvancedPRaEMixingRules::kLinear, "Linear"}, + {AdvancedPRaEMixingRules::kQuadratic, "Quadratic"}, +}) + struct AdvancedPRaEOptions{ AdvancedPRaEMixingRules brule = AdvancedPRaEMixingRules::kQuadratic; double s = 2.0; double CEoS = -sqrt(2.0)/2.0*log(1.0 + sqrt(2.0)); }; +inline void from_json(const json& j, AdvancedPRaEOptions& o) { + j.at("brule").get_to(o.brule); + j.at("s").get_to(o.s); + j.at("CEoS").get_to(o.CEoS); +} /** A residual Helmholtz term that returns nothing (the empty term) @@ -522,7 +549,7 @@ class WilsonResidualHelmholtzOverRT { } }; -using ResidualHelmholtzOverRTOptions = std::variant, WilsonResidualHelmholtzOverRT>; +using ResidualHelmholtzOverRTVariant = std::variant, WilsonResidualHelmholtzOverRT>; /** Cubic EOS with advanced mixing rules, the EoS/aE method of Jaubert and co-workers @@ -538,17 +565,21 @@ class AdvancedPRaEres { const NumType OmegaA = 0.45723552892138218938; const NumType OmegaB = 0.077796073903888455972; const int superanc_code = CubicSuperAncillary::PR_CODE; - const double CEoS; + protected: + std::valarray Tc_K, pc_Pa; + std::valarray ai, bi; const AlphaFunctions alphas; - const ResidualHelmholtzOverRTOptions ares; + const ResidualHelmholtzOverRTVariant ares; Eigen::ArrayXXd lmat; - const double s; + const AdvancedPRaEMixingRules brule; + const double s; + const double CEoS; nlohmann::json meta; @@ -575,8 +606,8 @@ class AdvancedPRaEres { }; public: - AdvancedPRaEres(const std::valarray& Tc_K, const std::valarray& pc_Pa, const AlphaFunctions& alphas, const ResidualHelmholtzOverRTOptions& ares, const Eigen::ArrayXXd& lmat, const AdvancedPRaEOptions& options = {}) - : alphas(alphas), ares(ares), lmat(lmat), s(options.s), brule(options.brule), CEoS(options.CEoS) + AdvancedPRaEres(const std::valarray& Tc_K, const std::valarray& pc_Pa, const AlphaFunctions& alphas, const ResidualHelmholtzOverRTVariant& ares, const Eigen::ArrayXXd& lmat, const AdvancedPRaEOptions& options = {}) + : Tc_K(Tc_K), pc_Pa(pc_Pa), alphas(alphas), ares(ares), lmat(lmat), brule(options.brule), s(options.s), CEoS(options.CEoS) { ai.resize(Tc_K.size()); bi.resize(Tc_K.size()); @@ -590,6 +621,9 @@ class AdvancedPRaEres { void set_meta(const nlohmann::json& j) { meta = j; } auto get_meta() const { return meta; } auto get_lmat() const { return lmat; } + auto get_Tc_K() const { return Tc_K; } + auto get_pc_Pa() const { return pc_Pa; } + static double get_bi(double Tc_K, double pc_Pa){ const NumType OmegaB = 0.077796073903888455972; const NumType R = 8.31446261815324; @@ -691,6 +725,38 @@ class AdvancedPRaEres { } }; +inline auto make_AdvancedPRaEres(const nlohmann::json& j){ + + std::valarray Tc_K = j.at("Tcrit / K"); + std::valarray pc_Pa = j.at("pcrit / Pa"); + + std::vector alphas = build_alpha_functions(Tc_K, j.at("alphas")); + + auto get_ares_model = [&](const nlohmann::json& armodel) -> ResidualHelmholtzOverRTVariant { + + std::string type = armodel.at("type"); + if (type == "Wilson"){ + std::vector b; + for (auto i = 0; i < Tc_K.size(); ++i){ + b.push_back(teqp::AdvancedPRaEres::get_bi(Tc_K[i], pc_Pa[i])); + } + auto mWilson = build_square_matrix(armodel.at("m")); + auto nWilson = build_square_matrix(armodel.at("n")); + return WilsonResidualHelmholtzOverRT(b, mWilson, nWilson); + } + else{ + throw teqp::InvalidArgument("bad type of ares model: " + type); + } + }; + auto aresmodel = get_ares_model(j.at("aresmodel")); + + AdvancedPRaEOptions options = j.at("options"); + auto model = teqp::AdvancedPRaEres(Tc_K, pc_Pa, alphas, aresmodel, Eigen::ArrayXXd::Zero(2, 2), options); + return model; +} + +using advancedPRaEres_t = decltype(make_AdvancedPRaEres({})); + /** The quantum corrected Peng-Robinson model as developed in diff --git a/interface/CPP/teqp_impl_factory.cpp b/interface/CPP/teqp_impl_factory.cpp index d00e0b3a..860f8e89 100644 --- a/interface/CPP/teqp_impl_factory.cpp +++ b/interface/CPP/teqp_impl_factory.cpp @@ -23,6 +23,7 @@ namespace teqp { {"SRK", [](const nlohmann::json& spec){ return make_owned(make_canonicalSRK(spec));}}, {"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));}}, {"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 a0e36fcb..3ada3225 100644 --- a/interface/pybind11_wrapper.cpp +++ b/interface/pybind11_wrapper.cpp @@ -109,6 +109,7 @@ const std::type_index SW_EspindolaHeredia2009_i{std::type_index(typeid(SW_Espind const std::type_index EXP6_Kataoka1992_i{std::type_index(typeid(EXP6_Kataoka1992_t))}; 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))}; /** At runtime we can add additional model-specific methods that only apply for a particular model. We take in a Python-wrapped @@ -176,6 +177,15 @@ 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 == advancedPRaEres_i){ + setattr("get_a", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed(o).get_a(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj)); + setattr("get_b", MethodType(py::cpp_function([](py::object& o, double T, REArrayd& molefrac){ return get_typed(o).get_b(T, molefrac); }, "self"_a, "T"_a, "molefrac"_a), obj)); + setattr("superanc_rhoLV", MethodType(py::cpp_function([](py::object& o, double T, std::size_t i){ return get_typed(o).superanc_rhoLV(T, i); }, "self"_a, "T"_a, "ifluid"_a), 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_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); diff --git a/src/tests/catch_test_cubics.cxx b/src/tests/catch_test_cubics.cxx index 34c5fdf6..b8acb201 100644 --- a/src/tests/catch_test_cubics.cxx +++ b/src/tests/catch_test_cubics.cxx @@ -677,3 +677,17 @@ TEST_CASE("Advanced cubic EOS", "[AdvancedPR]"){ } } } + +TEST_CASE("Advanced cubic EOS w/ make_model", "[AdvancedPR]"){ + auto j = R"({ + "kind": "advancedPRaEres", + "model": { + "Tcrit / K": [304.21, 126.19], + "pcrit / Pa": [7.383e6, 3395800.0], + "alphas": [{"type": "PR78", "acentric": 0.22394}, {"type": "PR78", "acentric": 0.0372}], + "aresmodel": {"type": "Wilson", "m": [[0.0, -3.4768], [3.5332, 0.0]], "n": [[0.0, 825], [-585, 0.0]]}, + "options": {"s": 2.0, "brule": "Quadratic", "CEoS": -0.52398} + } + })"_json; + auto model = make_model(j); +}