diff --git a/doc/source/models/CPA.ipynb b/doc/source/models/CPA.ipynb index 47e2d936..62cd599c 100644 --- a/doc/source/models/CPA.ipynb +++ b/doc/source/models/CPA.ipynb @@ -7,8 +7,64 @@ "source": [ "# Cubic Plus Association (CPA)\n", "\n", - "The combination of a cubic EOS with association" + "The combination of a cubic EOS with association with the association term. The sum of the terms goes like:\n", + "$$\n", + "\\alpha^{\\rm r} = \\alpha^{\\rm r}_{\\rm cub} + \\alpha^{\\rm r}_{\\rm assoc}\n", + "$$\n", + "\n", + "## Cubic part\n", + "\n", + "The residual contribution to $\\alpha$ is expressed as the sum :\n", + "$$\n", + "\\alpha^{\\rm r}_{\\rm cub,rep} +\\alpha^{\\rm r}_{\\rm cub,att}\n", + "$$ where the cubic parts come from\n", + "\n", + "The repulsive part of the cubic EOS contribution:\n", + "$$\n", + "\\alpha^{\\rm r}_{\\rm cub,rep} = -\\ln(1 - b_{\\rm mix}\\rho) \n", + "$$\n", + "The attractive part of the cubic EOS contribution:\n", + "$$\n", + "\\alpha^{\\rm r}_{\\rm cub,att} = -\\frac{a_{\\rm mix}}{RT}\\dfrac{\\ln\\left(\\frac{\\Delta_1 b_{\\rm mix}\\rho + 1}{\\Delta_2b_{\\rm mix}\\rho + 1}\\right)}{b_{\\rm mix}\\cdot(\\Delta_1 - \\Delta_2)}\n", + "$$\n", + "with the coefficients depending on the cubic type:\n", + "\n", + "SRK: $\\Delta_1=1$, $\\Delta_2=0$\n", + "\n", + "PR: $\\Delta_1=1+\\sqrt{2}$, $\\Delta_2=1-\\sqrt{2}$\n", + "\n", + "The mixture models used for the $a_{\\rm mix}$ and $b_{\\rm mix}$ are the classical ones:\n", + "\n", + "$$\n", + "a_{\\rm mix} = \\sum_i\\sum_jx_ix_j(1-k_{ij})a_{ij}(T)\n", + "$$\n", + "with x the mole fraction, $k_{ij}$ a weighting parameter\n", + "$$\n", + "a_{ij}(T) = \\sqrt{a_ia_j}\n", + "$$\n", + "and\n", + "$$\n", + "a_{i}(T) = a_{0i}\\left[1+c_{1i}(1-\\sqrt{T/T_{{\\rm crit},i}})\\right]^2\n", + "$$\n", + "and for $b$:\n", + "$$\n", + "b_{\\rm mix} = \\sum_ix_ib_i\n", + "$$\n", + "so there are three cubic parameters per fluid that need to be obtained though fitting: $b_{i}$, $a_{0i}$, $c_{1i}$. The value of $a_{\\rm ij}$ depends on temperature while $b_{\\rm mix}$ does not.\n", + "\n", + "## Association part\n", + "\n", + "For the association, one must have a solid understanding of the association approach that is being applied. To this end, a short discussion of the general approach is required. \n", + "\n" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2180ef67", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -27,7 +83,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.4" + "version": "3.10.13" } }, "nbformat": 4, diff --git a/doc/source/models/Molecule.drawio.svg b/doc/source/models/Molecule.drawio.svg new file mode 100644 index 00000000..c8d3a371 --- /dev/null +++ b/doc/source/models/Molecule.drawio.svg @@ -0,0 +1,4 @@ + + + +
B
B
B
B
A
A
C
C
A
A
A
AText is not SVG - cannot display \ No newline at end of file diff --git a/include/teqp/models/CPA.hpp b/include/teqp/models/CPA.hpp index 1ad88242..7bafb265 100644 --- a/include/teqp/models/CPA.hpp +++ b/include/teqp/models/CPA.hpp @@ -3,6 +3,9 @@ #include "nlohmann/json.hpp" #include #include "teqp/types.hpp" +#include "teqp/exceptions.hpp" +#include "teqp/models/association/association.hpp" +#include "teqp/models/association/association_types.hpp" namespace teqp { @@ -11,29 +14,10 @@ namespace CPA { template auto POW2(X x) { return x * x; }; template auto POW3(X x) { return x * POW2(x); }; -enum class association_classes {not_set, a1A, a2B, a3B, a4C, not_associating}; - -inline auto get_association_classes(const std::string& s) { - if (s == "1A") { return association_classes::a1A; } - else if (s == "2B") { return association_classes::a2B; } - else if (s == "2B") { return association_classes::a2B; } - else if (s == "3B") { return association_classes::a3B; } - else if (s == "4C") { return association_classes::a4C; } - else { - throw std::invalid_argument("bad association flag:" + s); - } -} - -enum class radial_dist { CS, KG, OT }; - -inline auto get_radial_dist(const std::string& s) { - if (s == "CS") { return radial_dist::CS; } - else if (s == "KG") { return radial_dist::KG; } - else if (s == "OT") { return radial_dist::OT; } - else { - throw std::invalid_argument("bad association flag:" + s); - } -} +using association::radial_dist; +using association::association_classes; +using association::get_radial_dist; +using association::get_association_classes; /// Function that calculates the association binding strength between site A of molecule i and site B on molecule j template @@ -57,10 +41,10 @@ inline auto get_DeltaAB_pure(radial_dist dist, double epsABi, double betaABi, BT g_vm_ref = 1.0 / (1.0 - 1.9 * eta); break; } - case radial_dist::OT: { - g_vm_ref = 1.0 / (1.0 - 0.475 * rhomolar * b_cubic); - break; - } +// case radial_dist::OT: { // This is identical to KG +// g_vm_ref = 1.0 / (1.0 - 0.475 * rhomolar * b_cubic); +// break; +// } default: { throw std::invalid_argument("Bad radial_dist"); } @@ -136,13 +120,12 @@ inline auto get_cubic_flag(const std::string& s) { class CPACubic { private: - std::valarray a0, bi, c1, Tc; + const std::valarray a0, bi, c1, Tc; double delta_1, delta_2; - std::valarray> k_ij; - double R_gas; - + const double R_gas; + const std::optional>> kmat; public: - CPACubic(cubic_flag flag, const std::valarray &a0, const std::valarray &bi, const std::valarray &c1, const std::valarray &Tc, double R_gas) : a0(a0), bi(bi), c1(c1), Tc(Tc), R_gas(R_gas) { + CPACubic(cubic_flag flag, const std::valarray &a0, const std::valarray &bi, const std::valarray &c1, const std::valarray &Tc, const double R_gas, const std::optional>> & kmat = std::nullopt) : a0(a0), bi(bi), c1(c1), Tc(Tc), R_gas(R_gas), kmat(kmat) { switch (flag) { case cubic_flag::PR: { delta_1 = 1 + sqrt(2); delta_2 = 1 - sqrt(2); break; } @@ -151,8 +134,9 @@ class CPACubic { default: throw std::invalid_argument("Bad cubic flag"); } - k_ij.resize(Tc.size()); for (auto i = 0U; i < k_ij.size(); ++i) { k_ij[i].resize(Tc.size()); } }; + + std::size_t size() const {return a0.size(); } template auto R(const VecType& /*molefrac*/) const { return R_gas; } @@ -171,7 +155,8 @@ class CPACubic { auto ai = get_ai(T, i); for (auto j = 0U; j < molefrac.size(); ++j) { auto aj = get_ai(T, j); - auto a_ij = (1.0 - k_ij[i][j]) * sqrt(ai * aj); + double kij = (kmat) ? kmat.value()[i][j] : 0.0; + auto a_ij = (1.0 - kij) * sqrt(ai * aj); asummer += molefrac[i] * molefrac[j] * a_ij; } } @@ -188,14 +173,14 @@ class CPACubic { } }; -template +/** Implement the association approach of Huang & Radosz for pure fluids + */ class CPAAssociation { private: - const Cubic cubic; const std::vector classes; const radial_dist dist; - const std::valarray epsABi, betaABi; - const std::vector N_sites; + const std::valarray epsABi, betaABi, bi; + const std::vector N_sites; const double R_gas; auto get_N_sites(const std::vector &the_classes) { @@ -216,13 +201,13 @@ class CPAAssociation { } public: - CPAAssociation(const Cubic &&cubic, const std::vector& classes, const radial_dist dist, const std::valarray &epsABi, const std::valarray &betaABi, double R_gas) - : cubic(cubic), classes(classes), dist(dist), epsABi(epsABi), betaABi(betaABi), N_sites(get_N_sites(classes)), R_gas(R_gas) {}; + CPAAssociation(const std::vector& classes, const radial_dist dist, const std::valarray &epsABi, const std::valarray &betaABi, const std::valarray &bi, double R_gas) + : classes(classes), dist(dist), epsABi(epsABi), betaABi(betaABi), bi(bi), N_sites(get_N_sites(classes)), R_gas(R_gas) {}; template auto alphar(const TType& T, const RhoType& rhomolar, const VecType& molefrac) const { - // Calculate a and b of the mixture - auto [a_cubic, b_cubic] = cubic.get_ab(T, molefrac); + // Calculate b of the mixture + auto b_cubic = (Eigen::Map(&bi[0], bi.size())*molefrac).sum(); // Calculate the fraction of sites not bonded with other active sites auto RT = forceeval(R_gas * T); // R times T @@ -259,6 +244,9 @@ class CPAEOS { /// alphar = a/(R*T) where a and R are both molar quantities template auto alphar(const TType& T, const RhoType& rhomolar, const VecType& molefrac) const { + if (static_cast(molefrac.size()) != cubic.size()){ + throw teqp::InvalidArgument("Mole fraction size is not correct; should be " + std::to_string(cubic.size())); + } // Calculate the contribution to alphar from the conventional cubic EOS auto alpha_r_cubic = cubic.alphar(T, rhomolar, molefrac); @@ -270,12 +258,43 @@ class CPAEOS { } }; +struct AssociationVariantWrapper{ + using vartype = std::variant; + const vartype holder; + + AssociationVariantWrapper(const vartype& holder) : holder(holder) {}; + + template + auto alphar(const TType& T, const RhoType& rhomolar, const MoleFracsType& molefracs) const{ + return std::visit([&](auto& h){ return h.alphar(T, rhomolar, molefracs); }, holder); + } +}; + /// A factory function to return an instantiated CPA instance given /// the JSON representation of the model inline auto CPAfactory(const nlohmann::json &j){ auto build_cubic = [](const auto& j) { auto N = j["pures"].size(); std::valarray a0i(N), bi(N), c1(N), Tc(N); + std::vector> kmat; + if (j.contains("kmat")){ + kmat = j.at("kmat"); + std::string kmaterr = "The kmat is the wrong size. It should be square with dimension " + std::to_string(N); + if (kmat.size() != N){ + throw teqp::InvalidArgument(kmaterr); + } + else{ + for (auto& krow: kmat){ + if(krow.size() != N){ + throw teqp::InvalidArgument(kmaterr); + } + } + } + } + else{ + kmat.resize(N); for (auto i = 0U; i < N; ++i){ kmat[i].resize(N); for (auto j = 0U; j < N; ++j){kmat[i][j] = 0.0;} } + } + std::size_t i = 0; for (auto p : j["pures"]) { a0i[i] = p["a0i / Pa m^6/mol^2"]; @@ -284,23 +303,79 @@ inline auto CPAfactory(const nlohmann::json &j){ Tc[i] = p["Tc / K"]; i++; } - return CPACubic(get_cubic_flag(j["cubic"]), a0i, bi, c1, Tc, j["R_gas / J/mol/K"]); + return CPACubic(get_cubic_flag(j["cubic"]), a0i, bi, c1, Tc, j["R_gas / J/mol/K"], kmat); }; - auto build_assoc = [](const auto &&cubic, const auto& j) { + + auto build_assoc_pure = [](const auto& j) -> AssociationVariantWrapper{ auto N = j["pures"].size(); - std::vector classes; - radial_dist dist = get_radial_dist(j.at("radial_dist")); - std::valarray epsABi(N), betaABi(N); - std::size_t i = 0; - for (auto p : j["pures"]) { - epsABi[i] = p["epsABi / J/mol"]; - betaABi[i] = p["betaABi"]; - classes.push_back(get_association_classes(p["class"])); - i++; + if (N == 1 && j.at("pures").contains("class") ){ + // This is the backwards compatible approach + // with the old style of defining the association class {1,2B...} + std::vector classes; + radial_dist dist = get_radial_dist(j.at("radial_dist")); + std::valarray epsABi(N), betaABi(N), bi(N); + std::size_t i = 0; + for (auto p : j.at("pures")) { + epsABi[i] = p.at("epsABi / J/mol"); + betaABi[i] = p.at("betaABi"); + bi[i] = p.at("bi / m^3/mol"); + classes.push_back(get_association_classes(p.at("class"))); + i++; + } + return AssociationVariantWrapper{CPAAssociation(classes, dist, epsABi, betaABi, bi, j["R_gas / J/mol/K"])}; + } + else{ + // This branch uses the new code + Eigen::ArrayXd b_m3mol(N), beta(N), epsilon_Jmol(N); + association::AssociationOptions opt; + opt.radial_dist = get_radial_dist(j.at("radial_dist")); + if (j.contains("options")){ + opt = j.at("options"); // Pulls in the options that are POD types + } + + std::vector> molecule_sites; + std::size_t i = 0; + std::set unique_site_types; + for (auto p : j["pures"]) { + epsilon_Jmol[i] = p.at("epsABi / J/mol"); + beta[i] = p.at("betaABi"); + b_m3mol[i] = p.at("bi / m^3/mol"); + molecule_sites.push_back(p.at("sites")); + for (auto & s : molecule_sites.back()){ + unique_site_types.insert(s); + } + i++; + } + if (j.contains("options") && j.at("options").contains("interaction_partners")){ + opt.interaction_partners = j.at("options").at("interaction_partners"); + for (auto [k,partners] : opt.interaction_partners){ + if (unique_site_types.count(k) == 0){ + throw teqp::InvalidArgument("Site is invalid in interaction_partners: " + k); + } + for (auto& partner : partners){ + if (unique_site_types.count(partner) == 0){ + throw teqp::InvalidArgument("Partner " + partner + " is invalid for site " + k); + } + } + } + } + else{ + // Every site type is assumed to interact with every other site type, except for itself + for (auto& site1 : unique_site_types){ + std::vector partners; + for (auto& site2: unique_site_types){ + if (site1 != site2){ + partners.push_back(site2); + } + } + opt.interaction_partners[site1] = partners; + } + } + + return AssociationVariantWrapper{association::Association(b_m3mol, beta, epsilon_Jmol, molecule_sites, opt)}; } - return CPAAssociation(std::move(cubic), classes, dist, epsABi, betaABi, j["R_gas / J/mol/K"]); }; - return CPAEOS(build_cubic(j), build_assoc(build_cubic(j), j)); + return CPAEOS(build_cubic(j), build_assoc_pure( j)); } }; /* namespace CPA */ diff --git a/include/teqp/models/association/association.hpp b/include/teqp/models/association/association.hpp new file mode 100644 index 00000000..7872a59a --- /dev/null +++ b/include/teqp/models/association/association.hpp @@ -0,0 +1,269 @@ +/** +General routines for the calculation of association + +The implementation follows the approach of Langenbach for the index compression, + + Many helpful hints from Andres Riedemann + +*/ + +#pragma once + +#include +#include + +#include "teqp/constants.hpp" +#include "teqp/types.hpp" + +#include + +#include "teqp/models/association/association_types.hpp" + +namespace teqp{ +namespace association{ + +struct AssociationOptions{ + std::map> interaction_partners; + std::vector site_order; + association::radial_dist radial_dist; + double alpha = 0.5; + double rtol = 1e-12, atol = 1e-12; + int max_iters = 10; +}; +inline void from_json(const nlohmann::json& j, AssociationOptions& o) { + if (j.contains("alpha")){ j.at("alpha").get_to(o.alpha); } + if (j.contains("rtol")){ j.at("rtol").get_to(o.rtol); } + if (j.contains("atol")){ j.at("atol").get_to(o.atol); } + if (j.contains("max_iters")){ j.at("max_iters").get_to(o.max_iters); } +} + +/*** + A general class for calculating association fractions and association energy for mixtures + + A mixture is formed of multiple components, each component has an index i. + For each component, it may multiple unique sites, and each unique site has a multiplicity associated with it. + + For instance, for 4C water, there are two "e" sites and two "H" sites in the nomenclature of Clapeyron. + Thus the "e" site has multiplicity of 2 and the "H" site has multiplicity of 2 + */ +class Association{ +public: + using CompSite = std::tuple; + using CompCIndex = std::tuple; + struct IndexMapper{ + std::map to_CompSite; + std::map to_siteid; + std::map CompCIndex_to_siteid; + Eigen::ArrayXi counts; ///< An array of counts of each siteid for the entire mixture + Eigen::ArrayXi N_sites; ///< How many total sites are on each molecule + Eigen::ArrayXi N_unique_sites; ///< How many unique types of sites are on each molecule + Eigen::ArrayXi comp_index; ///< The pure component indices associated with each siteid + std::vector> molecule_sites; + }; +private: + IndexMapper make_mapper(const std::vector>& molecule_sites, const AssociationOptions& options) const { + IndexMapper ind; + ind.counts.resize(1000); + ind.comp_index.resize(1000); + ind.N_sites.resize(molecule_sites.size()); + ind.N_unique_sites.resize(molecule_sites.size()); + ind.molecule_sites = molecule_sites; + + // Build the maps from siteid->(component index, site name) and (component index, site name)->siteid + std::size_t icomp = 0; + std::size_t siteid = 0; + for (auto& molecule: molecule_sites){ + /// Count up how many times each site is present in the list of sites + std::map site_counts; + for (auto& site: molecule){ + ++site_counts[site]; + } + auto unique_sites_on_molecule = std::set(molecule.begin(), molecule.end()); + if (!options.site_order.empty()){ + // TODO: enforce sites to appear in the order matching the specification + // TODO: this would be required to for instance check the D matrix of Langenbach and Enders + } + std::size_t Cindex = 0; + for (auto& site: unique_sites_on_molecule){ + CompSite cs{icomp, site}; + ind.to_CompSite[siteid] = cs; + ind.to_siteid[cs] = siteid; + ind.CompCIndex_to_siteid[CompCIndex{icomp, Cindex}] = siteid; + ind.counts[siteid] = site_counts[site]; + ind.comp_index[siteid] = static_cast(icomp); + Cindex++; + siteid++; + } + ind.N_sites[icomp] = static_cast(molecule.size()); + ind.N_unique_sites[icomp] = static_cast(unique_sites_on_molecule.size()); + icomp++; + } + + ind.counts.conservativeResize(siteid); + ind.comp_index.conservativeResize(siteid); + return ind; + } + + /*** + Construct the counting matrix \f$ D_{IJ} \f$ as given by Langenbach and Enders + */ + auto make_D(const IndexMapper& ind, const AssociationOptions& options ) const{ + + auto get_DIJ = [&ind, &options](std::size_t I, std::size_t J) -> int { + /** Return the value of an entry in the D_{IJ} matrix + + For a given unique site, look at all other sites on all other molecules + */ + auto [ph1, typei] = ind.to_CompSite.at(I); + auto [ph2, typej] = ind.to_CompSite.at(J); + auto contains = [](auto& container, const auto& val){ return std::find(container.begin(), container.end(), val) != container.end(); }; + /// If interaction parameters are not provided, assume conservatively that all sites can interact with all other sites + if (options.interaction_partners.empty() || (contains(options.interaction_partners.at(typei), typej))){ + return ind.counts[J]; + } + return 0; + }; + + int Ngroups = static_cast(ind.to_siteid.size()); + Eigen::ArrayXXi D(Ngroups, Ngroups); + // I and J are the numerical indices of the sites + for (int I = 0; I < Ngroups; ++I){ + for (int J = 0; J < Ngroups; ++J){ + D(I, J) = get_DIJ(I, J); + } + } + return D; + } + +public: + const Eigen::ArrayXd b_m3mol, ///< The covolume b, in m^3/mol + beta, ///< The volume factor, dimensionless + epsilon_Jmol; ///< The association energy of each molecule, in J/mol + const AssociationOptions options; + const IndexMapper mapper; + const Eigen::ArrayXXi D; + const radial_dist m_radial_dist; + + Association(const Eigen::ArrayXd& b_m3mol, const Eigen::ArrayXd& beta, const Eigen::ArrayXd& epsilon_Jmol, const std::vector>& molecule_sites, const AssociationOptions& options) : b_m3mol(b_m3mol), beta(beta), epsilon_Jmol(epsilon_Jmol), options(options), mapper(make_mapper(molecule_sites, options)), D(make_D(mapper, options)), m_radial_dist(options.radial_dist){ + } + + /** + Build the Delta matrix, where entries are given by + \f[ + \Delta_{IJ} = ... + \f] + */ + template + auto get_Delta(const TType& T, const RhoType& rhomolar, const MoleFracsType& molefracs) const { + + using resulttype = std::common_type_t; // Type promotion, without the const-ness + using Mat = Eigen::Array; + auto Ngroups = mapper.to_CompSite.size(); + auto bmix = (molefracs*b_m3mol).sum(); + auto eta = bmix*rhomolar/4.0; + + decltype(forceeval(eta)) g; + switch(m_radial_dist){ + case radial_dist::CS: + g = (2.0-eta)/(2.0*(1.0-eta)*(1.0-eta)*(1.0-eta)); break; + case radial_dist::KG: + g = 1.0 / (1.0 - 1.9*eta); break; + default: + throw std::invalid_argument("Bad radial distribution"); + } + + Mat Delta = Mat::Zero(Ngroups, Ngroups); + for (auto I = 0U; I < Ngroups; ++I){ + auto i = std::get<0>(mapper.to_CompSite.at(I)); + for (auto J = 0U; J < Ngroups; ++J){ + auto j = std::get<0>(mapper.to_CompSite.at(J)); + + using namespace teqp::constants; + // The CR1 rule is used to calculate the cross contributions + if (true){ // Is CR1 // TODO: also allow the other combining rule + auto b_ij = (b_m3mol[i] + b_m3mol[j])/2.0; + auto epsilon_ij_Jmol = (epsilon_Jmol[i] + epsilon_Jmol[j])/2.0; + auto beta_ij = sqrt(beta[i]*beta[j]); + Delta(I, J) = g*b_ij*beta_ij*(exp(epsilon_ij_Jmol/(R_CODATA2017*T))-1.0)/N_A; + } + } + } + return Delta; + } + + template + auto successive_substitution(const TType& T, const RhoType& rhomolar, const MoleFracsType& molefracs, const XType& X_init) const { + + using resulttype = std::common_type_t; // Type promotion, without the const-ness + using Mat = Eigen::Array; + + const Mat Delta = get_Delta(T, rhomolar, molefracs); +// const Mat DD = Delta.array()*D.cast().array(); // coefficient-wise product, with D upcast from int to type of Delta matrix + + auto Ngroups = mapper.to_CompSite.size(); + Eigen::RowVectorX xj(Ngroups); // Mole fractions of the component containing each siteid + for (auto I = 0U; I< Ngroups; ++I){ + xj(I) = molefracs(std::get<0>(mapper.to_CompSite.at(I))); + } + + using rDDXtype = std::decay_t>; // Type promotion, without the const-ness + Eigen::MatrixX rDDX = rhomolar*N_A*(Delta.array()*D.cast().array()).matrix(); + for (auto j = 0; j < rDDX.rows(); ++j){ + rDDX.row(j).array() = rDDX.row(j).array()*xj.array(); + } +// rDDX.rowwise() *= xj; + + Eigen::ArrayX> X = X_init, Xnew; + + for (auto counter = 0; counter < options.max_iters; ++counter){ + // calculate the new array of non-bonded site fractions X + Xnew = options.alpha*X + (1.0-options.alpha)/(1.0+(rDDX*X.matrix()).array()); + // These unaryExpr extract the numerical value from an Eigen array of generic type, allowing for comparison. + // Otherwise for instance it is imposible to compare two complex numbers (if you are using complex step derivatives) + auto diff = (Xnew-X).eval().cwiseAbs().unaryExpr([](const auto&x){return getbaseval(x); }).eval(); + auto tol = (options.rtol*Xnew + options.atol).unaryExpr([](const auto&x){return getbaseval(x); }).eval(); + if ((diff < tol).all()){ + break; + } + X = Xnew; + } + return X; + } + + /** + \brief Calculate the contribution \f$\alpha = a/(RT)\f$, where the Helmholtz energy \f$a\f$ is on a molar basis, making \f$\alpha\f$ dimensionless. + */ + template + auto alphar(const TType& T, const RhoType& rhomolar, const MoleFracsType& molefracs) const { + + // Do the sucessive substitution to obtain the non-bonded fractions for each unique site + Eigen::ArrayXd X_init = Eigen::ArrayXd::Ones(mapper.to_siteid.size()); + auto X_A = successive_substitution(T, rhomolar, molefracs, X_init); + + // Calculate the contribution alpha based on the "naive" summation like in Clapeyron + using resulttype = std::common_type_t; // Type promotion, without the const-ness + resulttype alpha_r_asso = 0.0; + for (auto icomp = 0; icomp < molefracs.size(); ++icomp){ // loop over all components + resulttype summer = 0.0; + for (auto jsite = 0; jsite < mapper.N_unique_sites[icomp]; ++jsite){ + // Get the siteid given the (component index, site index on the component) + const auto I = mapper.CompCIndex_to_siteid.at({icomp, jsite}); + + // each contribution (per unique site on a molecule) is multiplied by its multiplicity within the molecule + summer += (log(X_A(I))- X_A(I)/2.0 + 0.5)*static_cast(mapper.counts(I)); + } + alpha_r_asso += molefracs[icomp]*summer; + } + return alpha_r_asso; + +// // And do the summation to calculate the contribution to alpha=a/(RT) +// using V = std::common_type_t; +// auto xj = molefracs(mapper.comp_index); // mole fractions of the components attached at each site +// auto alpha_r_asso = (xj.template cast()*(log(X_A) - X_A/2.0 + 0.5).template cast()*mapper.counts.cast()).sum(); +// return forceeval(alpha_r_asso); + } + +}; + +} +} diff --git a/include/teqp/models/association/association_types.hpp b/include/teqp/models/association/association_types.hpp new file mode 100644 index 00000000..ab49a302 --- /dev/null +++ b/include/teqp/models/association/association_types.hpp @@ -0,0 +1,30 @@ +#pragma once + +namespace teqp { +namespace association{ + +enum class association_classes {not_set, a1A, a2B, a3B, a4C, not_associating}; + +inline auto get_association_classes(const std::string& s) { + if (s == "1A") { return association_classes::a1A; } + else if (s == "2B") { return association_classes::a2B; } + else if (s == "2B") { return association_classes::a2B; } + else if (s == "3B") { return association_classes::a3B; } + else if (s == "4C") { return association_classes::a4C; } + else { + throw std::invalid_argument("bad association flag:" + s); + } +} + +enum class radial_dist { CS, KG }; + +inline auto get_radial_dist(const std::string& s) { + if (s == "CS") { return radial_dist::CS; } + else if (s == "KG") { return radial_dist::KG; } + else { + throw std::invalid_argument("bad association flag:" + s); + } +} + +} +} diff --git a/src/tests/catch_test_association.cxx b/src/tests/catch_test_association.cxx new file mode 100644 index 00000000..12f5ab9c --- /dev/null +++ b/src/tests/catch_test_association.cxx @@ -0,0 +1,142 @@ +#include +#include +#include +using Catch::Approx; + +#include + +#include "teqp/cpp/teqpcpp.hpp" +#include "teqp/models/association/association.hpp" + +using namespace teqp; + +TEST_CASE("Test making the indices and D like in Langenbach", "[association]"){ + auto b_m3mol = (Eigen::ArrayXd(2) << 0.0491/1e3, 0.0145/1e3).finished(); + auto beta = (Eigen::ArrayXd(2) << 8e-3, 69.2e-3).finished(); + auto epsilon_Jmol = (Eigen::ArrayXd(2) << 215.00*100, 166.55*100).finished(); + std::vector> molecule_sites = {{"B"},{"N","N","P"},{"P"}}; + association::AssociationOptions opt; + opt.interaction_partners = {{"B", {"N", "P", "B"}}, {"N", {"P", "B"}}, {"P", {"N", "B"}}}; + opt.site_order = {"B","P","N"}; + + association::Association a(b_m3mol, beta, epsilon_Jmol, molecule_sites, opt); +} + +TEST_CASE("Test ethanol + water association", "[association]"){ + auto b_m3mol = (Eigen::ArrayXd(2) << 0.0491/1e3, 0.0145/1e3).finished(); + auto beta = (Eigen::ArrayXd(2) << 8e-3, 69.2e-3).finished(); + auto epsilon_Jmol = (Eigen::ArrayXd(2) << 215.00*100, 166.55*100).finished(); + + std::vector> molecule_sites = {{"e", "H"}, {"e", "e", "H", "H"}}; + association::AssociationOptions opt; + opt.interaction_partners = {{"e", {"H",}}, {"H", {"e",}}}; + + association::Association a(b_m3mol, beta, epsilon_Jmol, molecule_sites, opt); + + auto molefracs = (Eigen::ArrayXd(2) << 0.3, 0.7).finished(); + auto Ngroups = a.mapper.to_siteid.size(); + Eigen::ArrayXd X_init = Eigen::ArrayXd::Ones(Ngroups); + + double v = 3.0680691201961814e-5; + Eigen::ArrayXd X_A = a.successive_substitution(303.15, 1/v, molefracs, X_init); + CHECK(X_A[0] == Approx(0.06258400385436955)); + CHECK(X_A[3] == Approx(0.10938445109190545)); + + BENCHMARK("SS"){ + return a.successive_substitution(303.15, 1/v, molefracs, X_init); + }; + BENCHMARK("alphar"){ + return a.alphar(303.15, 1/v, molefracs); + }; +} + +TEST_CASE("Ethanol + water with CPA", "[association]"){ + nlohmann::json water = { + {"a0i / Pa m^6/mol^2", 0.12277}, {"bi / m^3/mol", 0.0000145}, {"c1", 0.6736}, {"Tc / K", 647.13}, + {"epsABi / J/mol", 16655.0}, {"betaABi", 0.0692}, {"sites", {"e","e","H","H"}} + }; + nlohmann::json ethanol = { + {"a0i / Pa m^6/mol^2", 0.85164}, {"bi / m^3/mol", 0.0491e-3}, {"c1", 0.7502}, {"Tc / K", 513.92}, + {"epsABi / J/mol", 21500.0}, {"betaABi", 0.008}, {"sites", {"e","H"}} + }; + nlohmann::json jCPA = { + {"cubic", "SRK"}, + {"radial_dist", "CS"}, +// {"combining", "CR1"}, + {"options", {}}, + {"pures", {ethanol, water}}, + {"R_gas / J/mol/K", 8.31446261815324} + }; + nlohmann::json j = { + {"kind", "CPA"}, + {"model", jCPA} + }; + auto model = teqp::cppinterface::make_model(j, false); + + auto molefracs = (Eigen::ArrayXd(2) << 0.3, 0.7).finished(); + double T = 303.15, rhomolar=1/3.0680691201961814e-5; + double R = model->get_R(molefracs); + + auto alphar = model->get_Ar00(303.15, rhomolar, molefracs); + CHECK(alphar == Approx(-8.333844120879878)); + + double p = T*R*rhomolar*(1+model->get_Ar01(T, rhomolar, molefracs)); + CHECK(p == Approx(1e5)); + BENCHMARK("p(T,rho)"){ + return T*R*rhomolar*(1+model->get_Ar01(T, rhomolar, molefracs)); + }; +} + +TEST_CASE("Trace ethanol + water isobaric VLE with CPA", "[associationVLE]"){ + nlohmann::json water = { + {"a0i / Pa m^6/mol^2", 0.12277}, {"bi / m^3/mol", 0.0000145}, {"c1", 0.6736}, {"Tc / K", 647.13}, + {"epsABi / J/mol", 16655.0}, {"betaABi", 0.0692}, {"sites", {"e","e","H","H"}} + }; + nlohmann::json ethanol = { + {"a0i / Pa m^6/mol^2", 0.85164}, {"bi / m^3/mol", 0.0491e-3}, {"c1", 0.7502}, {"Tc / K", 513.92}, + {"epsABi / J/mol", 21500.0}, {"betaABi", 0.008}, {"sites", {"e","H"}} + }; + nlohmann::json options = { + {"alpha", 0.4999} + }; + nlohmann::json jCPA = { + {"cubic", "SRK"}, + {"radial_dist", "CS"}, + {"options", options}, +// {"combining", "CR1"}, + {"pures", {ethanol, water}}, + {"R_gas / J/mol/K", 8.31446261815324} + }; + + nlohmann::json j = { + {"kind", "CPA"}, + {"validate", false}, + {"model", jCPA} + }; + auto model = teqp::cppinterface::make_model(j, false); + + nlohmann::json jeth = { + {"kind", "CPA"}, + {"validate", false}, + {"model", { + {"cubic", "SRK"}, + {"radial_dist", "CS"}, + {"pures", {ethanol}}, + {"R_gas / J/mol/K", 8.31446261815324} + }} + }; + auto modeleth = teqp::cppinterface::make_model(jeth, false); + + double T = 351; + auto rhoLrhoV = modeleth->pure_VLE_T(T, 15985.159939940693, 35.82761304572266, 30); + double rhoL = rhoLrhoV[0]; + auto z = (Eigen::ArrayXd(1) << 1.0).finished(); + auto p = rhoL*modeleth->R(z)*T*(1+modeleth->get_Ar01(T, rhoL, z)); + + auto rhovecL = (Eigen::ArrayXd(2) << rhoLrhoV[0], 0).finished(); + auto rhovecV = (Eigen::ArrayXd(2) << rhoLrhoV[1], 0).finished(); + + PVLEOptions opt; opt.verbosity = 100; + auto o = model->trace_VLE_isobar_binary(p, T, rhovecL, rhovecV, opt); + CHECK(o.size() > 10); +} diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx index d6151356..4298dcb7 100644 --- a/src/tests/catch_tests.cxx +++ b/src/tests/catch_tests.cxx @@ -450,7 +450,7 @@ TEST_CASE("Test water Clapeyron.jl", "[CPA]") { std::vector schemes = { CPA::association_classes::a4C }; std::valarray epsAB = { 16655 }, betaAB = { 0.0692 }; CPA::radial_dist dist = CPA::radial_dist::CS; - CPA::CPAAssociation cpaa(std::move(cub), schemes, dist, epsAB, betaAB, R); + CPA::CPAAssociation cpaa(schemes, dist, epsAB, betaAB, bi, R); CPA::CPAEOS cpa(std::move(cub), std::move(cpaa)); using tdc = TDXDerivatives; @@ -487,7 +487,7 @@ TEST_CASE("Test water", "[CPA]") { std::vector schemes = { CPA::association_classes::a4C }; std::valarray epsAB = { 16655 }, betaAB = { 0.0692 }; - CPA::CPAAssociation cpaa(std::move(cub), schemes, CPA::radial_dist::KG, epsAB, betaAB, R); + CPA::CPAAssociation cpaa(schemes, CPA::radial_dist::KG, epsAB, betaAB, bi, R); CPA::CPAEOS cpa(std::move(cub), std::move(cpaa)); using tdc = TDXDerivatives; @@ -512,7 +512,7 @@ TEST_CASE("Test [C2mim][NTF2]", "[CPA]") { std::vector schemes = { CPA::association_classes::a2B }; std::valarray epsAB = { 177388/10.0 }, betaAB = { 0.00266 }; - CPA::CPAAssociation cpaa(std::move(cub), schemes, CPA::radial_dist::KG, epsAB, betaAB, R); + CPA::CPAAssociation cpaa(schemes, CPA::radial_dist::KG, epsAB, betaAB, bi, R); CPA::CPAEOS cpa(std::move(cub), std::move(cpaa)); using tdc = TDXDerivatives;