From a0fcb6f26c7630973af2b2acca1254ff28970ed8 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Fri, 28 Apr 2023 16:50:04 -0400
Subject: [PATCH] Implementation of Gubbins and Twu + Luckas
Appears to be working, additional testing required. Generalizing to handle G&V and G&T correlation integrals remains
---
include/teqp/models/pcsaft.hpp | 258 +---------
.../models/saft/correlation_integrals.hpp | 18 +-
include/teqp/models/saft/polar_terms.hpp | 447 ++++++++++++++++++
include/teqp/models/saftvrmie.hpp | 15 +-
src/tests/catch_test_SAFTpolar.cxx | 43 +-
5 files changed, 511 insertions(+), 270 deletions(-)
create mode 100644 include/teqp/models/saft/polar_terms.hpp
diff --git a/include/teqp/models/pcsaft.hpp b/include/teqp/models/pcsaft.hpp
index cec7537..324c134 100644
--- a/include/teqp/models/pcsaft.hpp
+++ b/include/teqp/models/pcsaft.hpp
@@ -10,6 +10,7 @@
#include "teqp/types.hpp"
#include "teqp/exceptions.hpp"
#include "teqp/constants.hpp"
+#include "teqp/models/saft/polar_terms.hpp"
#include <optional>
namespace teqp {
@@ -146,74 +147,6 @@ auto get_I2(const Eta& eta, MbarType mbar) {
return std::make_tuple(forceeval(summer_I2), forceeval(summer_etadI2deta));
}
-/// Eq. 10 from Gross and Vrabec
-template <typename Eta, typename MType, typename TType>
-auto get_JDD_2ij(const Eta& eta, const MType& mij, const TType& Tstarij) {
- static Eigen::ArrayXd a_0 = (Eigen::ArrayXd(5) << 0.3043504, -0.1358588, 1.4493329, 0.3556977, -2.0653308).finished();
- static Eigen::ArrayXd a_1 = (Eigen::ArrayXd(5) << 0.9534641, -1.8396383, 2.0131180, -7.3724958, 8.2374135).finished();
- static Eigen::ArrayXd a_2 = (Eigen::ArrayXd(5) << -1.1610080, 4.5258607, 0.9751222, -12.281038, 5.9397575).finished();
-
- static Eigen::ArrayXd b_0 = (Eigen::ArrayXd(5) << 0.2187939, -1.1896431, 1.1626889, 0, 0).finished();
- static Eigen::ArrayXd b_1 = (Eigen::ArrayXd(5) << -0.5873164, 1.2489132, -0.5085280, 0, 0).finished();
- static Eigen::ArrayXd b_2 = (Eigen::ArrayXd(5) << 3.4869576, -14.915974, 15.372022, 0, 0).finished();
-
- std::common_type_t<Eta, MType> summer = 0.0;
- for (auto n = 0; n < 5; ++n){
- auto anij = a_0[n] + (mij-1)/mij*a_1[n] + (mij-1)/mij*(mij-2)/mij*a_2[n]; // Eq. 12
- auto bnij = b_0[n] + (mij-1)/mij*b_1[n] + (mij-1)/mij*(mij-2)/mij*b_2[n]; // Eq. 13
- summer += (anij + bnij/Tstarij)*pow(eta, n);
- }
- return forceeval(summer);
-}
-
-/// Eq. 11 from Gross and Vrabec
-template <typename Eta, typename MType>
-auto get_JDD_3ijk(const Eta& eta, const MType& mijk) {
- static Eigen::ArrayXd c_0 = (Eigen::ArrayXd(5) << -0.0646774, 0.1975882, -0.8087562, 0.6902849, 0.0).finished();
- static Eigen::ArrayXd c_1 = (Eigen::ArrayXd(5) << -0.9520876, 2.9924258, -2.3802636, -0.2701261, 0.0).finished();
- static Eigen::ArrayXd c_2 = (Eigen::ArrayXd(5) << -0.6260979, 1.2924686, 1.6542783, -3.4396744, 0.0).finished();
- std::common_type_t<Eta, MType> summer = 0.0;
- for (auto n = 0; n < 5; ++n){
- auto cnijk = c_0[n] + (mijk-1)/mijk*c_1[n] + (mijk-1)/mijk*(mijk-2)/mijk*c_2[n]; // Eq. 14
- summer += cnijk*pow(eta, n);
- }
- return forceeval(summer);
-}
-
-/// Eq. 12 from Gross and Vrabec, AICHEJ
-template <typename Eta, typename MType, typename TType>
-auto get_JQQ_2ij(const Eta& eta, const MType& mij, const TType& Tstarij) {
- static Eigen::ArrayXd a_0 = (Eigen::ArrayXd(5) << 1.2378308, 2.4355031, 1.6330905, -1.6118152, 6.9771185).finished();
- static Eigen::ArrayXd a_1 = (Eigen::ArrayXd(5) << 1.2854109, -11.465615, 22.086893, 7.4691383, -17.197772).finished();
- static Eigen::ArrayXd a_2 = (Eigen::ArrayXd(5) << 1.7942954, 0.7695103, 7.2647923, 94.486699, -77.148458).finished();
-
- static Eigen::ArrayXd b_0 = (Eigen::ArrayXd(5) << 0.4542718, -4.5016264, 3.5858868, 0.0, 0.0).finished();
- static Eigen::ArrayXd b_1 = (Eigen::ArrayXd(5) << -0.8137340, 10.064030, -10.876631, 0.0, 0.0).finished();
- static Eigen::ArrayXd b_2 = (Eigen::ArrayXd(5) << 6.8682675, -5.1732238, -17.240207, 0.0, 0.0).finished();
-
- std::common_type_t<Eta, MType> summer = 0.0;
- for (auto n = 0; n < 5; ++n){
- auto anij = a_0[n] + (mij-1)/mij*a_1[n] + (mij-1)/mij*(mij-2)/mij*a_2[n]; // Eq. 12
- auto bnij = b_0[n] + (mij-1)/mij*b_1[n] + (mij-1)/mij*(mij-2)/mij*b_2[n]; // Eq. 13
- summer += (anij + bnij/Tstarij)*pow(eta, n);
- }
- return forceeval(summer);
-}
-
-/// Eq. 13 from Gross and Vrabec, AICHEJ
-template <typename Eta, typename MType>
-auto get_JQQ_3ijk(const Eta& eta, const MType& mijk) {
- static Eigen::ArrayXd c_0 = (Eigen::ArrayXd(5) << 0.5000437, 6.5318692, -16.014780, 14.425970, 0.0).finished();
- static Eigen::ArrayXd c_1 = (Eigen::ArrayXd(5) << 2.0002094, -6.7838658, 20.383246, -10.895984, 0.0).finished();
- static Eigen::ArrayXd c_2 = (Eigen::ArrayXd(5) << 3.1358271, 7.2475888, 3.0759478, 0.0, 0.0).finished();
- std::common_type_t<Eta, MType> summer = 0.0;
- for (auto n = 0; n < 5; ++n){
- auto cnijk = c_0[n] + (mijk-1)/mijk*c_1[n] + (mijk-1)/mijk*(mijk-2)/mijk*c_2[n]; // Eq. 14
- summer += cnijk*pow(eta, n);
- }
- return forceeval(summer);
-}
-
/**
Sum up three array-like objects that can each have different container types and value types
*/
@@ -337,197 +270,16 @@ public:
}
};
-/***
- * \brief The dipolar contribution given in Gross and Vrabec
- */
-class PCSAFTDipolarContribution {
-private:
- const Eigen::ArrayXd m, sigma_Angstrom, epsilon_over_k, mustar2, nmu;
- template<typename A> auto POW2(const A& x) const { return forceeval(x*x); }
- template<typename A> auto POW3(const A& x) const { return forceeval(POW2(x)*x); }
-public:
- PCSAFTDipolarContribution(const Eigen::ArrayX<double> &m, const Eigen::ArrayX<double> &sigma_Angstrom, const Eigen::ArrayX<double> &epsilon_over_k, const Eigen::ArrayX<double> &mustar2, const Eigen::ArrayX<double> &nmu) : m(m), sigma_Angstrom(sigma_Angstrom), epsilon_over_k(epsilon_over_k), mustar2(mustar2), nmu(nmu) {
- // Check lengths match
- if (m.size() != mustar2.size()){
- throw teqp::InvalidArgument("bad size of mustar2");
- }
- if (m.size() != nmu.size()){
- throw teqp::InvalidArgument("bad size of n");
- }
- }
-
- /// Eq. 8 from Gross and Vrabec
- template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
- auto get_alpha2DD(const TTYPE& T, const RhoType& rhoN_A3, const EtaType& eta, const VecType& mole_fractions) const{
- const auto& x = mole_fractions; // concision
- const auto& sigma = sigma_Angstrom; // concision
- const auto N = mole_fractions.size();
- std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summer = 0.0;
- for (auto i = 0; i < N; ++i){
- for (auto j = 0; j < N; ++j){
- auto ninj = nmu[i]*nmu[j];
- if (ninj > 0){
- // Lorentz-Berthelot mixing rules
- auto epskij = sqrt(epsilon_over_k[i]*epsilon_over_k[j]);
- auto sigmaij = (sigma[i]+sigma[j])/2;
-
- auto Tstarij = T/epskij;
- auto mij = std::max(sqrt(m[i]*m[j]), 2);
- summer += x[i]*x[j]*epsilon_over_k[i]/T*epsilon_over_k[j]/T*POW3(sigma[i]*sigma[j]/sigmaij)*ninj*mustar2[i]*mustar2[j]*get_JDD_2ij(eta, mij, Tstarij);
- }
- }
- }
- return forceeval(-static_cast<double>(EIGEN_PI)*rhoN_A3*summer);
- }
-
- /// Eq. 9 from Gross and Vrabec
- template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
- auto get_alpha3DD(const TTYPE& T, const RhoType& rhoN_A3, const EtaType& eta, const VecType& mole_fractions) const{
- const auto& x = mole_fractions; // concision
- const auto& sigma = sigma_Angstrom; // concision
- const auto N = mole_fractions.size();
- std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summer = 0.0;
- for (auto i = 0; i < N; ++i){
- for (auto j = 0; j < N; ++j){
- for (auto k = 0; k < N; ++k){
- auto ninjnk = nmu[i]*nmu[j]*nmu[k];
- if (ninjnk > 0){
- // Lorentz-Berthelot mixing rules for sigma
- auto sigmaij = (sigma[i]+sigma[j])/2;
- auto sigmaik = (sigma[i]+sigma[k])/2;
- auto sigmajk = (sigma[j]+sigma[k])/2;
-
- auto mijk = std::max(pow(m[i]*m[j]*m[k], 1.0/3.0), 2.0);
- summer += x[i]*x[j]*x[k]*epsilon_over_k[i]/T*epsilon_over_k[j]/T*epsilon_over_k[k]/T*POW3(sigma[i]*sigma[j]*sigma[k])/(sigmaij*sigmaik*sigmajk)*ninjnk*mustar2[i]*mustar2[j]*mustar2[k]*get_JDD_3ijk(eta, mijk);
- }
- }
- }
- }
- return forceeval(-4.0*POW2(static_cast<double>(EIGEN_PI))/3.0*POW2(rhoN_A3)*summer);
- }
-
- /***
- * \brief Get the dipolar contribution to \f$ \alpha = A/(NkT) \f$
- */
- template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
- auto eval(const TTYPE& T, const RhoType& rho_A3, const EtaType& eta, const VecType& mole_fractions) const {
- auto alpha2 = get_alpha2DD(T, rho_A3, eta, mole_fractions);
- auto alpha3 = get_alpha3DD(T, rho_A3, eta, mole_fractions);
- auto alpha = forceeval(alpha2/(1.0-alpha3/alpha2));
-
- using alpha2_t = decltype(alpha2);
- using alpha3_t = decltype(alpha3);
- using alpha_t = decltype(alpha);
- struct PCSAFTDipolarContributionTerms{
- alpha2_t alpha2;
- alpha3_t alpha3;
- alpha_t alpha;
- };
- return PCSAFTDipolarContributionTerms{alpha2, alpha3, alpha};
- }
-};
-
-/***
- * \brief The quadrupolar contribution from Gross and Vrabec
- *
- */
-class PCSAFTQuadrupolarContribution {
-private:
- const Eigen::ArrayXd m, sigma_Angstrom, epsilon_over_k, Qstar2, nQ;
- template<typename A> auto POW2(const A& x) const { return forceeval(x*x); }
- template<typename A> auto POW3(const A& x) const { return forceeval(POW2(x)*x); }
- template<typename A> auto POW5(const A& x) const { return forceeval(POW2(x)*POW3(x)); }
- template<typename A> auto POW7(const A& x) const { return forceeval(POW2(x)*POW5(x)); }
-public:
- PCSAFTQuadrupolarContribution(const Eigen::ArrayX<double> &m, const Eigen::ArrayX<double> &sigma_Angstrom, const Eigen::ArrayX<double> &epsilon_over_k, const Eigen::ArrayX<double> &Qstar2, const Eigen::ArrayX<double> &nQ) : m(m), sigma_Angstrom(sigma_Angstrom), epsilon_over_k(epsilon_over_k), Qstar2(Qstar2), nQ(nQ) {
- // Check lengths match
- if (m.size() != Qstar2.size()){
- throw teqp::InvalidArgument("bad size of mustar2");
- }
- if (m.size() != nQ.size()){
- throw teqp::InvalidArgument("bad size of n");
- }
- }
- PCSAFTQuadrupolarContribution& operator=( const PCSAFTQuadrupolarContribution& ) = delete; // non copyable
-
- /// Eq. 9 from Gross and Vrabec
- template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
- auto get_alpha2QQ(const TTYPE& T, const RhoType& rhoN_A3, const EtaType& eta, const VecType& mole_fractions) const{
- const auto& x = mole_fractions; // concision
- const auto& sigma = sigma_Angstrom; // concision
- const auto N = mole_fractions.size();
- std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summer = 0.0;
- for (auto i = 0; i < N; ++i){
- for (auto j = 0; j < N; ++j){
- auto ninj = nQ[i]*nQ[j];
- if (ninj > 0){
- // Lorentz-Berthelot mixing rules
- auto epskij = sqrt(epsilon_over_k[i]*epsilon_over_k[j]);
- auto sigmaij = (sigma[i]+sigma[j])/2;
-
- auto Tstarij = T/epskij;
- auto mij = sqrt(m[i]*m[j]);
- summer += x[i]*x[j]*epsilon_over_k[i]/T*epsilon_over_k[j]/T*POW5(sigma[i]*sigma[j])/POW7(sigmaij)*ninj*Qstar2[i]*Qstar2[j]*get_JQQ_2ij(eta, mij, Tstarij);
- }
- }
- }
- return forceeval(-static_cast<double>(EIGEN_PI)*POW2(3.0/4.0)*rhoN_A3*summer);
- }
-
- /// Eq. 1 from Gross and Vrabec
- template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
- auto get_alpha3QQ(const TTYPE& T, const RhoType& rhoN_A3, const EtaType& eta, const VecType& mole_fractions) const{
- const auto& x = mole_fractions; // concision
- const auto& sigma = sigma_Angstrom; // concision
- const auto N = mole_fractions.size();
- std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summer = 0.0;
- for (auto i = 0; i < N; ++i){
- for (auto j = 0; j < N; ++j){
- for (auto k = 0; k < N; ++k){
- auto ninjnk = nQ[i]*nQ[j]*nQ[k];
- if (ninjnk > 0){
- // Lorentz-Berthelot mixing rules for sigma
- auto sigmaij = (sigma[i]+sigma[j])/2;
- auto sigmaik = (sigma[i]+sigma[k])/2;
- auto sigmajk = (sigma[j]+sigma[k])/2;
-
- auto mijk = pow(m[i]*m[j]*m[k], 1.0/3.0);
- summer += x[i]*x[j]*x[k]*epsilon_over_k[i]/T*epsilon_over_k[j]/T*epsilon_over_k[k]/T*POW5(sigma[i]*sigma[j]*sigma[k])/POW3(sigmaij*sigmaik*sigmajk)*ninjnk*Qstar2[i]*Qstar2[j]*Qstar2[k]*get_JDD_3ijk(eta, mijk);
- }
- }
- }
- }
- return forceeval(-4.0*POW2(static_cast<double>(EIGEN_PI))/3.0*POW3(3.0/4.0)*POW2(rhoN_A3)*summer);
- }
-
- /***
- * \brief Get the quadrupolar contribution to \f$ \alpha = A/(NkT) \f$
- */
- template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
- auto eval(const TTYPE& T, const RhoType& rho_A3, const EtaType& eta, const VecType& mole_fractions) const {
- auto alpha2 = get_alpha2QQ(T, rho_A3, eta, mole_fractions);
- auto alpha3 = get_alpha3QQ(T, rho_A3, eta, mole_fractions);
- auto alpha = forceeval(alpha2/(1.0-alpha3/alpha2));
-
- using alpha2_t = decltype(alpha2);
- using alpha3_t = decltype(alpha3);
- using alpha_t = decltype(alpha);
- struct PCSAFTQuadrupolarContributionTerms{
- alpha2_t alpha2;
- alpha3_t alpha3;
- alpha_t alpha;
- };
- return PCSAFTQuadrupolarContributionTerms{alpha2, alpha3, alpha};
- }
-};
-
/** A class used to evaluate mixtures using PC-SAFT model
This is the classical Gross and Sadowski model from 2001: https://doi.org/10.1021/ie0003887
-with the error fixed as noted in a comment: https://doi.org/10.1021/acs.iecr.9b01515
+with the errors fixed as noted in a comment: https://doi.org/10.1021/acs.iecr.9b01515
*/
class PCSAFTMixture {
+public:
+ using PCSAFTDipolarContribution = SAFTpolar::DipolarContributionGrossVrabec;
+ using PCSAFTQuadrupolarContribution = SAFTpolar::QuadrupolarContributionGrossVrabec;
protected:
Eigen::ArrayX<double> m, ///< number of segments
mminus1, ///< m-1
diff --git a/include/teqp/models/saft/correlation_integrals.hpp b/include/teqp/models/saft/correlation_integrals.hpp
index d2fa92c..6529a1c 100644
--- a/include/teqp/models/saft/correlation_integrals.hpp
+++ b/include/teqp/models/saft/correlation_integrals.hpp
@@ -21,14 +21,14 @@ static const std::map<int, std::array<double, 12>> Luckas_J_coeffs = {
{15, {-2.14335965e-01, 2.24240261e-02, 2.68094773e-01, 8.11899188e00, 2.15506735e-01, -8.11465705e00, -1.88310645e01, -4.18309476e-01, 1.88679367e01, 1.02033085e01, 2.37674032e-01, -1.00120648e01}},
};
-template<int n>
+
class LuckasJIntegral{
public:
-
+ const int n;
const std::array<double, 12> a;
double a00,a01,a02,a10,a11,a12,a20,a21,a22,a30,a31,a32;
- LuckasJIntegral() : a(Luckas_J_coeffs.at(n)){
+ LuckasJIntegral(int n) : n(n), a(Luckas_J_coeffs.at(n)){
a00 = a[0]; a01 = a[1]; a02 = a[2];
a10 = a[3]; a11 = a[4]; a12 = a[5];
a20 = a[6]; a21 = a[7]; a22 = a[8];
@@ -36,13 +36,15 @@ public:
}
template<typename TType, typename RhoType>
- auto get_J(const TType& Tstar, const RhoType& rhostar){
+ auto get_J(const TType& Tstar, const RhoType& rhostar) const{
auto Z_1 = 0.3 + 0.05*n;
auto Z_2 = 1.0/n;
auto A_0 = a00 + a10*rhostar + a20*rhostar*rhostar + a30*rhostar*rhostar*rhostar;
auto A_1 = a01 + a11*rhostar + a21*rhostar*rhostar + a31*rhostar*rhostar*rhostar;
auto A_2 = a02 + a12*rhostar + a22*rhostar*rhostar + a32*rhostar*rhostar*rhostar;
- std::common_type_t<TType, RhoType> out = (A_0 + A_1*pow(Tstar, Z_1) + A_2*pow(Tstar, Z_2))*exp(1.0/(Tstar + 4.0/pow(std::abs(log(rhostar/sqrt(2.0))), 3.0)));
+ // |x| = sqrt(x^2), the latter is well-suited to differentiation
+ auto differentiable_abs = [](const auto& x){ return sqrt(x*x); };
+ std::common_type_t<TType, RhoType> out = (A_0 + A_1*pow(Tstar, Z_1) + A_2*pow(Tstar, Z_2))*exp(1.0/(Tstar + 4.0/pow(differentiable_abs(log(rhostar/sqrt(2.0))), 3.0)));
return out;
}
};
@@ -54,14 +56,14 @@ static const std::map<std::tuple<int, int>, std::array<double, 16>> Luckas_K_coe
{{444,555},{ 1.07862851e-03, 2.88647260e-04, 2.07369614e-03, 1.79600129e-04, 7.57084957e-02, -2.28761197e-03, -7.33670248e-02, 3.12736172e-02, 1.73518964e-01, 4.54174801e-03, -5.14488122e-01, 2.85483011e-01, 1.71670939e-01, -1.07503182e-03, -5.51880129e-01, 4.16887229e-01}}
};
-template<int n1, int n2>
class LuckasKIntegral{
public:
+ const int n1, n2;
const std::array<double, 16> a;
double a00, a01, a02, a03, a10, a11, a12, a13, a20, a21, a22, a23, a30, a31, a32, a33;
- LuckasKIntegral() : a(Luckas_K_coeffs.at({n1, n2})){
+ LuckasKIntegral(const int n1, const int n2) : n1(n1), n2(n2), a(Luckas_K_coeffs.at({n1, n2})){
a00 = a[0]; a01 = a[1]; a02 = a[2]; a03 = a[3];
a10 = a[4]; a11 = a[5]; a12 = a[6]; a13 = a[7];
a20 = a[8]; a21 = a[9]; a22 = a[10]; a23 = a[11];
@@ -69,7 +71,7 @@ public:
}
template<typename TType, typename RhoType>
- auto get_K(const TType& Tstar, const RhoType& rhostar){
+ auto get_K(const TType& Tstar, const RhoType& rhostar) const{
auto Z_1 = 2.0;
auto Z_2 = 3.0;
auto Z_3 = 4.0;
diff --git a/include/teqp/models/saft/polar_terms.hpp b/include/teqp/models/saft/polar_terms.hpp
new file mode 100644
index 0000000..173a2ce
--- /dev/null
+++ b/include/teqp/models/saft/polar_terms.hpp
@@ -0,0 +1,447 @@
+#pragma once
+
+/**
+ This header contains methods that pertain to polar contributions to SAFT models
+
+ Initially the contribution of Gross and Vrabec were implemented for PC-SAFT, but they can be used with other
+ non-polar base models as well, so this header collects all the things in one place
+ */
+
+#include "teqp/types.hpp"
+#include "teqp/exceptions.hpp"
+#include "correlation_integrals.hpp"
+#include <Eigen/Dense>
+
+namespace teqp{
+
+namespace SAFTpolar{
+
+/// Eq. 10 from Gross and Vrabec
+template <typename Eta, typename MType, typename TType>
+auto get_JDD_2ij(const Eta& eta, const MType& mij, const TType& Tstarij) {
+ static Eigen::ArrayXd a_0 = (Eigen::ArrayXd(5) << 0.3043504, -0.1358588, 1.4493329, 0.3556977, -2.0653308).finished();
+ static Eigen::ArrayXd a_1 = (Eigen::ArrayXd(5) << 0.9534641, -1.8396383, 2.0131180, -7.3724958, 8.2374135).finished();
+ static Eigen::ArrayXd a_2 = (Eigen::ArrayXd(5) << -1.1610080, 4.5258607, 0.9751222, -12.281038, 5.9397575).finished();
+
+ static Eigen::ArrayXd b_0 = (Eigen::ArrayXd(5) << 0.2187939, -1.1896431, 1.1626889, 0, 0).finished();
+ static Eigen::ArrayXd b_1 = (Eigen::ArrayXd(5) << -0.5873164, 1.2489132, -0.5085280, 0, 0).finished();
+ static Eigen::ArrayXd b_2 = (Eigen::ArrayXd(5) << 3.4869576, -14.915974, 15.372022, 0, 0).finished();
+
+ std::common_type_t<Eta, MType> summer = 0.0;
+ for (auto n = 0; n < 5; ++n){
+ auto anij = a_0[n] + (mij-1)/mij*a_1[n] + (mij-1)/mij*(mij-2)/mij*a_2[n]; // Eq. 12
+ auto bnij = b_0[n] + (mij-1)/mij*b_1[n] + (mij-1)/mij*(mij-2)/mij*b_2[n]; // Eq. 13
+ summer += (anij + bnij/Tstarij)*pow(eta, n);
+ }
+ return forceeval(summer);
+}
+
+/// Eq. 11 from Gross and Vrabec
+template <typename Eta, typename MType>
+auto get_JDD_3ijk(const Eta& eta, const MType& mijk) {
+ static Eigen::ArrayXd c_0 = (Eigen::ArrayXd(5) << -0.0646774, 0.1975882, -0.8087562, 0.6902849, 0.0).finished();
+ static Eigen::ArrayXd c_1 = (Eigen::ArrayXd(5) << -0.9520876, 2.9924258, -2.3802636, -0.2701261, 0.0).finished();
+ static Eigen::ArrayXd c_2 = (Eigen::ArrayXd(5) << -0.6260979, 1.2924686, 1.6542783, -3.4396744, 0.0).finished();
+ std::common_type_t<Eta, MType> summer = 0.0;
+ for (auto n = 0; n < 5; ++n){
+ auto cnijk = c_0[n] + (mijk-1)/mijk*c_1[n] + (mijk-1)/mijk*(mijk-2)/mijk*c_2[n]; // Eq. 14
+ summer += cnijk*pow(eta, n);
+ }
+ return forceeval(summer);
+}
+
+/// Eq. 12 from Gross and Vrabec, AICHEJ
+template <typename Eta, typename MType, typename TType>
+auto get_JQQ_2ij(const Eta& eta, const MType& mij, const TType& Tstarij) {
+ static Eigen::ArrayXd a_0 = (Eigen::ArrayXd(5) << 1.2378308, 2.4355031, 1.6330905, -1.6118152, 6.9771185).finished();
+ static Eigen::ArrayXd a_1 = (Eigen::ArrayXd(5) << 1.2854109, -11.465615, 22.086893, 7.4691383, -17.197772).finished();
+ static Eigen::ArrayXd a_2 = (Eigen::ArrayXd(5) << 1.7942954, 0.7695103, 7.2647923, 94.486699, -77.148458).finished();
+
+ static Eigen::ArrayXd b_0 = (Eigen::ArrayXd(5) << 0.4542718, -4.5016264, 3.5858868, 0.0, 0.0).finished();
+ static Eigen::ArrayXd b_1 = (Eigen::ArrayXd(5) << -0.8137340, 10.064030, -10.876631, 0.0, 0.0).finished();
+ static Eigen::ArrayXd b_2 = (Eigen::ArrayXd(5) << 6.8682675, -5.1732238, -17.240207, 0.0, 0.0).finished();
+
+ std::common_type_t<Eta, MType> summer = 0.0;
+ for (auto n = 0; n < 5; ++n){
+ auto anij = a_0[n] + (mij-1)/mij*a_1[n] + (mij-1)/mij*(mij-2)/mij*a_2[n]; // Eq. 12
+ auto bnij = b_0[n] + (mij-1)/mij*b_1[n] + (mij-1)/mij*(mij-2)/mij*b_2[n]; // Eq. 13
+ summer += (anij + bnij/Tstarij)*pow(eta, n);
+ }
+ return forceeval(summer);
+}
+
+/// Eq. 13 from Gross and Vrabec, AICHEJ
+template <typename Eta, typename MType>
+auto get_JQQ_3ijk(const Eta& eta, const MType& mijk) {
+ static Eigen::ArrayXd c_0 = (Eigen::ArrayXd(5) << 0.5000437, 6.5318692, -16.014780, 14.425970, 0.0).finished();
+ static Eigen::ArrayXd c_1 = (Eigen::ArrayXd(5) << 2.0002094, -6.7838658, 20.383246, -10.895984, 0.0).finished();
+ static Eigen::ArrayXd c_2 = (Eigen::ArrayXd(5) << 3.1358271, 7.2475888, 3.0759478, 0.0, 0.0).finished();
+ std::common_type_t<Eta, MType> summer = 0.0;
+ for (auto n = 0; n < 5; ++n){
+ auto cnijk = c_0[n] + (mijk-1)/mijk*c_1[n] + (mijk-1)/mijk*(mijk-2)/mijk*c_2[n]; // Eq. 14
+ summer += cnijk*pow(eta, n);
+ }
+ return forceeval(summer);
+}
+
+/***
+ * \brief The dipolar contribution given in Gross and Vrabec
+ */
+class DipolarContributionGrossVrabec {
+private:
+ const Eigen::ArrayXd m, sigma_Angstrom, epsilon_over_k, mustar2, nmu;
+ template<typename A> auto POW2(const A& x) const { return forceeval(x*x); }
+ template<typename A> auto POW3(const A& x) const { return forceeval(POW2(x)*x); }
+public:
+ DipolarContributionGrossVrabec(const Eigen::ArrayX<double> &m, const Eigen::ArrayX<double> &sigma_Angstrom, const Eigen::ArrayX<double> &epsilon_over_k, const Eigen::ArrayX<double> &mustar2, const Eigen::ArrayX<double> &nmu) : m(m), sigma_Angstrom(sigma_Angstrom), epsilon_over_k(epsilon_over_k), mustar2(mustar2), nmu(nmu) {
+ // Check lengths match
+ if (m.size() != mustar2.size()){
+ throw teqp::InvalidArgument("bad size of mustar2");
+ }
+ if (m.size() != nmu.size()){
+ throw teqp::InvalidArgument("bad size of n");
+ }
+ }
+
+ /// Eq. 8 from Gross and Vrabec
+ template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
+ auto get_alpha2DD(const TTYPE& T, const RhoType& rhoN_A3, const EtaType& eta, const VecType& mole_fractions) const{
+ const auto& x = mole_fractions; // concision
+ const auto& sigma = sigma_Angstrom; // concision
+ const auto N = mole_fractions.size();
+ std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summer = 0.0;
+ for (auto i = 0; i < N; ++i){
+ for (auto j = 0; j < N; ++j){
+ auto ninj = nmu[i]*nmu[j];
+ if (ninj > 0){
+ // Lorentz-Berthelot mixing rules
+ auto epskij = sqrt(epsilon_over_k[i]*epsilon_over_k[j]);
+ auto sigmaij = (sigma[i]+sigma[j])/2;
+
+ auto Tstarij = T/epskij;
+ auto mij = std::max(sqrt(m[i]*m[j]), 2.0);
+ summer += x[i]*x[j]*epsilon_over_k[i]/T*epsilon_over_k[j]/T*POW3(sigma[i]*sigma[j]/sigmaij)*ninj*mustar2[i]*mustar2[j]*get_JDD_2ij(eta, mij, Tstarij);
+ }
+ }
+ }
+ return forceeval(-static_cast<double>(EIGEN_PI)*rhoN_A3*summer);
+ }
+
+ /// Eq. 9 from Gross and Vrabec
+ template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
+ auto get_alpha3DD(const TTYPE& T, const RhoType& rhoN_A3, const EtaType& eta, const VecType& mole_fractions) const{
+ const auto& x = mole_fractions; // concision
+ const auto& sigma = sigma_Angstrom; // concision
+ const auto N = mole_fractions.size();
+ std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summer = 0.0;
+ for (auto i = 0; i < N; ++i){
+ for (auto j = 0; j < N; ++j){
+ for (auto k = 0; k < N; ++k){
+ auto ninjnk = nmu[i]*nmu[j]*nmu[k];
+ if (ninjnk > 0){
+ // Lorentz-Berthelot mixing rules for sigma
+ auto sigmaij = (sigma[i]+sigma[j])/2;
+ auto sigmaik = (sigma[i]+sigma[k])/2;
+ auto sigmajk = (sigma[j]+sigma[k])/2;
+
+ auto mijk = std::max(pow(m[i]*m[j]*m[k], 1.0/3.0), 2.0);
+ summer += x[i]*x[j]*x[k]*epsilon_over_k[i]/T*epsilon_over_k[j]/T*epsilon_over_k[k]/T*POW3(sigma[i]*sigma[j]*sigma[k])/(sigmaij*sigmaik*sigmajk)*ninjnk*mustar2[i]*mustar2[j]*mustar2[k]*get_JDD_3ijk(eta, mijk);
+ }
+ }
+ }
+ }
+ return forceeval(-4.0*POW2(static_cast<double>(EIGEN_PI))/3.0*POW2(rhoN_A3)*summer);
+ }
+
+ /***
+ * \brief Get the dipolar contribution to \f$ \alpha = A/(NkT) \f$
+ */
+ template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
+ auto eval(const TTYPE& T, const RhoType& rho_A3, const EtaType& eta, const VecType& mole_fractions) const {
+ auto alpha2 = get_alpha2DD(T, rho_A3, eta, mole_fractions);
+ auto alpha3 = get_alpha3DD(T, rho_A3, eta, mole_fractions);
+ auto alpha = forceeval(alpha2/(1.0-alpha3/alpha2));
+
+ using alpha2_t = decltype(alpha2);
+ using alpha3_t = decltype(alpha3);
+ using alpha_t = decltype(alpha);
+ struct DipolarContributionTerms{
+ alpha2_t alpha2;
+ alpha3_t alpha3;
+ alpha_t alpha;
+ };
+ return DipolarContributionTerms{alpha2, alpha3, alpha};
+ }
+};
+
+/***
+ * \brief The quadrupolar contribution from Gross and Vrabec
+ *
+ */
+class QuadrupolarContributionGrossVrabec {
+private:
+ const Eigen::ArrayXd m, sigma_Angstrom, epsilon_over_k, Qstar2, nQ;
+ template<typename A> auto POW2(const A& x) const { return forceeval(x*x); }
+ template<typename A> auto POW3(const A& x) const { return forceeval(POW2(x)*x); }
+ template<typename A> auto POW5(const A& x) const { return forceeval(POW2(x)*POW3(x)); }
+ template<typename A> auto POW7(const A& x) const { return forceeval(POW2(x)*POW5(x)); }
+public:
+ QuadrupolarContributionGrossVrabec(const Eigen::ArrayX<double> &m, const Eigen::ArrayX<double> &sigma_Angstrom, const Eigen::ArrayX<double> &epsilon_over_k, const Eigen::ArrayX<double> &Qstar2, const Eigen::ArrayX<double> &nQ) : m(m), sigma_Angstrom(sigma_Angstrom), epsilon_over_k(epsilon_over_k), Qstar2(Qstar2), nQ(nQ) {
+ // Check lengths match
+ if (m.size() != Qstar2.size()){
+ throw teqp::InvalidArgument("bad size of mustar2");
+ }
+ if (m.size() != nQ.size()){
+ throw teqp::InvalidArgument("bad size of n");
+ }
+ }
+ QuadrupolarContributionGrossVrabec& operator=( const QuadrupolarContributionGrossVrabec& ) = delete; // non copyable
+
+ /// Eq. 9 from Gross and Vrabec
+ template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
+ auto get_alpha2QQ(const TTYPE& T, const RhoType& rhoN_A3, const EtaType& eta, const VecType& mole_fractions) const{
+ const auto& x = mole_fractions; // concision
+ const auto& sigma = sigma_Angstrom; // concision
+ const auto N = mole_fractions.size();
+ std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summer = 0.0;
+ for (auto i = 0; i < N; ++i){
+ for (auto j = 0; j < N; ++j){
+ auto ninj = nQ[i]*nQ[j];
+ if (ninj > 0){
+ // Lorentz-Berthelot mixing rules
+ auto epskij = sqrt(epsilon_over_k[i]*epsilon_over_k[j]);
+ auto sigmaij = (sigma[i]+sigma[j])/2;
+
+ auto Tstarij = T/epskij;
+ auto mij = sqrt(m[i]*m[j]);
+ summer += x[i]*x[j]*epsilon_over_k[i]/T*epsilon_over_k[j]/T*POW5(sigma[i]*sigma[j])/POW7(sigmaij)*ninj*Qstar2[i]*Qstar2[j]*get_JQQ_2ij(eta, mij, Tstarij);
+ }
+ }
+ }
+ return forceeval(-static_cast<double>(EIGEN_PI)*POW2(3.0/4.0)*rhoN_A3*summer);
+ }
+
+ /// Eq. 1 from Gross and Vrabec
+ template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
+ auto get_alpha3QQ(const TTYPE& T, const RhoType& rhoN_A3, const EtaType& eta, const VecType& mole_fractions) const{
+ const auto& x = mole_fractions; // concision
+ const auto& sigma = sigma_Angstrom; // concision
+ const auto N = mole_fractions.size();
+ std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summer = 0.0;
+ for (auto i = 0; i < N; ++i){
+ for (auto j = 0; j < N; ++j){
+ for (auto k = 0; k < N; ++k){
+ auto ninjnk = nQ[i]*nQ[j]*nQ[k];
+ if (ninjnk > 0){
+ // Lorentz-Berthelot mixing rules for sigma
+ auto sigmaij = (sigma[i]+sigma[j])/2;
+ auto sigmaik = (sigma[i]+sigma[k])/2;
+ auto sigmajk = (sigma[j]+sigma[k])/2;
+
+ auto mijk = pow(m[i]*m[j]*m[k], 1.0/3.0);
+ summer += x[i]*x[j]*x[k]*epsilon_over_k[i]/T*epsilon_over_k[j]/T*epsilon_over_k[k]/T*POW5(sigma[i]*sigma[j]*sigma[k])/POW3(sigmaij*sigmaik*sigmajk)*ninjnk*Qstar2[i]*Qstar2[j]*Qstar2[k]*get_JDD_3ijk(eta, mijk);
+ }
+ }
+ }
+ }
+ return forceeval(-4.0*POW2(static_cast<double>(EIGEN_PI))/3.0*POW3(3.0/4.0)*POW2(rhoN_A3)*summer);
+ }
+
+ /***
+ * \brief Get the quadrupolar contribution to \f$ \alpha = A/(NkT) \f$
+ */
+ template<typename TTYPE, typename RhoType, typename EtaType, typename VecType>
+ auto eval(const TTYPE& T, const RhoType& rho_A3, const EtaType& eta, const VecType& mole_fractions) const {
+ auto alpha2 = get_alpha2QQ(T, rho_A3, eta, mole_fractions);
+ auto alpha3 = get_alpha3QQ(T, rho_A3, eta, mole_fractions);
+ auto alpha = forceeval(alpha2/(1.0-alpha3/alpha2));
+
+ using alpha2_t = decltype(alpha2);
+ using alpha3_t = decltype(alpha3);
+ using alpha_t = decltype(alpha);
+ struct QuadrupolarContributionTerms{
+ alpha2_t alpha2;
+ alpha3_t alpha3;
+ alpha_t alpha;
+ };
+ return QuadrupolarContributionTerms{alpha2, alpha3, alpha};
+ }
+};
+
+/**
+ \tparam JIntegral A type that can be indexed with a single integer n to give the J^{(n)} integral
+ \tparam KIntegral A type that can be indexed with a two integers a and b to give the K(a,b) integral
+
+ The flexibility was added to include J and K integrals from either Luckas et al. or Gubbins and Twu (or any others following the interface)
+ */
+template<class JIntegral, class KIntegral>
+class MultipolarContributionGubbinsTwu {
+private:
+ const Eigen::ArrayXd sigma_m, epsilon_over_k, mubar2, Qbar2;
+ template<typename A> auto POW2(const A& x) const { return forceeval(x*x); }
+ template<typename A> auto POW3(const A& x) const { return forceeval(POW2(x)*x); }
+ template<typename A> auto POW4(const A& x) const { return forceeval(POW2(x)*POW2(x)); }
+ template<typename A> auto POW5(const A& x) const { return forceeval(POW2(x)*POW3(x)); }
+ template<typename A> auto POW7(const A& x) const { return forceeval(POW2(x)*POW5(x)); }
+ template<typename A> auto POW8(const A& x) const { return forceeval(POW4(x)*POW4(x)); }
+ template<typename A> auto POW10(const A& x) const { return forceeval(POW2(x)*POW8(x)); }
+ template<typename A> auto POW12(const A& x) const { return forceeval(POW4(x)*POW8(x)); }
+
+ const JIntegral J6{6};
+ const JIntegral J8{8};
+ const JIntegral J10{10};
+ const JIntegral J11{11};
+ const JIntegral J13{13};
+ const JIntegral J15{15};
+ const KIntegral K222_333{222, 333};
+ const KIntegral K233_344{233, 344};
+ const KIntegral K334_445{334, 445};
+ const KIntegral K444_555{444, 555};
+ const double epsilon_0 = 8.8541878128e-12; // https://en.wikipedia.org/wiki/Vacuum_permittivity, in F/m, or C^2⋅N^−1⋅m^−2
+ const double PI_ = static_cast<double>(EIGEN_PI);
+
+public:
+ MultipolarContributionGubbinsTwu(const Eigen::ArrayX<double> &sigma_m, const Eigen::ArrayX<double> &epsilon_over_k, const Eigen::ArrayX<double> &mubar2, const Eigen::ArrayX<double> &Qbar2) : sigma_m(sigma_m), epsilon_over_k(epsilon_over_k), mubar2(mubar2), Qbar2(Qbar2) {
+ // Check lengths match
+ if (sigma_m.size() != mubar2.size()){
+ throw teqp::InvalidArgument("bad size of mubar2");
+ }
+ if (sigma_m.size() != Qbar2.size()){
+ throw teqp::InvalidArgument("bad size of Qbar2");
+ }
+ }
+ MultipolarContributionGubbinsTwu& operator=( const MultipolarContributionGubbinsTwu& ) = delete; // non copyable
+
+ template<typename TTYPE, typename RhoType, typename VecType>
+ auto get_alpha2(const TTYPE& T, const RhoType& rhoN, const RhoType& rhostar, const VecType& mole_fractions) const{
+ const auto& x = mole_fractions; // concision
+ const auto& sigma = sigma_m; // concision
+
+ const auto N = mole_fractions.size();
+ std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summer112 = 0.0, summer123 = 0.0, summer224 = 0.0;
+
+ for (auto i = 0; i < N; ++i){
+ for (auto j = 0; j < N; ++j){
+ // Lorentz-Berthelot mixing rules
+ auto epskij = sqrt(epsilon_over_k[i]*epsilon_over_k[j]);
+ auto sigmaij = (sigma[i]+sigma[j])/2;
+
+ auto Tstari = T/epsilon_over_k[i], Tstarj = T/epsilon_over_k[j];
+ auto leading = x[i]*x[j]/(Tstari*Tstarj); // common for all alpha_2 terms
+ auto Tstarij = forceeval(T/epskij);
+
+ summer112 += leading*POW3(sigma[i]*sigma[j])/POW3(sigmaij)*mubar2[i]*mubar2[j]*J6.get_J(Tstarij, rhostar);
+ summer123 += leading*POW3(sigma[i])*POW5(sigma[j])/POW5(sigmaij)*mubar2[i]*Qbar2[j]*J8.get_J(Tstarij, rhostar);
+ summer224 += leading*POW5(sigma[i]*sigma[j])/POW7(sigmaij)*Qbar2[i]*Qbar2[j]*J10.get_J(Tstarij, rhostar);
+ }
+ }
+
+ auto alpha2_112 = -2.0*PI_*rhoN*POW2(4*PI_*epsilon_0)/3.0*summer112;
+ auto alpha2_123 = -PI_*rhoN*POW2(4*PI_*epsilon_0)/3.0*summer123;
+ auto alpha2_224 = -14.0*PI_*rhoN*POW2(4*PI_*epsilon_0)/5.0*summer224;
+
+ return forceeval(alpha2_112 + 2.0*alpha2_123 + alpha2_224);
+ }
+
+ template<typename TTYPE, typename RhoType, typename VecType>
+ auto get_alpha3(const TTYPE& T, const RhoType& rhoN, const RhoType& rhostar, const VecType& mole_fractions) const{
+ const auto& x = mole_fractions; // concision
+ const auto& sigma = sigma_m; // concision
+ const auto N = mole_fractions.size();
+ std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summerA_112_112_224 = 0.0, summerA_112_123_213 = 0.0, summerA_123_123_224 = 0.0, summerA_224_224_224 = 0.0;
+ std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> summerB_112_112_112 = 0.0, summerB_112_123_123 = 0.0, summerB_123_123_224 = 0.0, summerB_224_224_224 = 0.0;
+ for (auto i = 0; i < N; ++i){
+ for (auto j = 0; j < N; ++j){
+
+ // Lorentz-Berthelot mixing rules
+ auto epskij = sqrt(epsilon_over_k[i]*epsilon_over_k[j]);
+ auto sigmaij = (sigma[i]+sigma[j])/2;
+
+ auto Tstari = T/epsilon_over_k[i], Tstarj = T/epsilon_over_k[j];
+ auto Tstarij = forceeval(T/epskij);
+
+ auto leading = x[i]*x[j]/pow(Tstari*Tstarj, 3.0/2.0); // common for all alpha_3A terms
+
+ summerA_112_112_224 += leading*pow(sigma[i]*sigma[j], 11.0/2.0)/POW8(sigmaij)*mubar2[i]*mubar2[j]*sqrt(Qbar2[i]*Qbar2[j])*J11.get_J(Tstarij, rhostar);
+ summerA_112_123_213 += leading*pow(sigma[i]*sigma[j], 11.0/2.0)/POW8(sigmaij)*mubar2[i]*mubar2[j]*sqrt(Qbar2[i]*Qbar2[j])*J11.get_J(Tstarij, rhostar);
+ summerA_123_123_224 += leading*pow(sigma[i], 11.0/2.0)*pow(sigma[j], 15.0/2.0)/POW10(sigmaij)*mubar2[i]*sqrt(Qbar2[i])*pow(Qbar2[j], 3.0/2.0)*J13.get_J(Tstarij, rhostar);
+ summerA_224_224_224 += leading*pow(sigma[i]*sigma[j], 15.0/2.0)/POW12(sigmaij)*Qbar2[i]*Qbar2[j]*J15.get_J(Tstarij, rhostar);
+
+ for (auto k = 0; k < N; ++k){
+ auto Tstark = T/epsilon_over_k[k];
+ auto epskik = sqrt(epsilon_over_k[i]*epsilon_over_k[k]);
+ auto epskjk = sqrt(epsilon_over_k[j]*epsilon_over_k[k]);
+ auto Tstarik = T/epskik;
+ auto Tstarjk = T/epskjk;
+
+ // Lorentz-Berthelot mixing rules for sigma
+ auto sigmaij = (sigma[i]+sigma[j])/2;
+ auto sigmaik = (sigma[i]+sigma[k])/2;
+ auto sigmajk = (sigma[j]+sigma[k])/2;
+ auto leadingijk = x[i]*x[j]*x[k]/(Tstari*Tstarj*Tstark);
+
+ auto K222333 = pow(K222_333.get_K(Tstarij, rhostar)*K222_333.get_K(Tstarik, rhostar)*K222_333.get_K(Tstarjk, rhostar), 1.0/3.0);
+ summerB_112_112_112 += leadingijk*POW3(sigma[i]*sigma[j]*sigma[k])/(sigmaij*sigmaik*sigmajk)*mubar2[i]*mubar2[j]*mubar2[k]*K222333;
+ auto K233344 = pow(K233_344.get_K(Tstarij, rhostar)*K233_344.get_K(Tstarik, rhostar)*K233_344.get_K(Tstarjk, rhostar), 1.0/3.0);
+ summerB_112_123_123 += leadingijk*POW3(sigma[i]*sigma[j])*POW5(sigma[k])/(sigmaij*POW2(sigmaik*sigmajk))*mubar2[i]*mubar2[j]*Qbar2[k]*K233344;
+ auto K334445 = pow(K334_445.get_K(Tstarij, rhostar)*K334_445.get_K(Tstarik, rhostar)*K334_445.get_K(Tstarjk, rhostar), 1.0/3.0);
+ summerB_123_123_224 += leadingijk*POW3(sigma[i])*POW5(sigma[j]*sigma[k])/(POW2(sigmaij*sigmaik)*POW3(sigmajk))*mubar2[i]*Qbar2[j]*Qbar2[k]*K334445;
+ auto K444555 = pow(K444_555.get_K(Tstarij, rhostar)*K444_555.get_K(Tstarik, rhostar)*K444_555.get_K(Tstarjk, rhostar), 1.0/3.0);
+ summerB_224_224_224 += leadingijk*POW5(sigma[i]*sigma[j]*sigma[k])/(POW3(sigmaij*sigmaik*sigmajk))*Qbar2[i]*Qbar2[j]*Qbar2[k]*K444555;
+ }
+ }
+ }
+ auto alpha3A_112_112_224 = 8.0*PI_*rhoN*POW3(4*PI_*epsilon_0)/25.0*summerA_112_112_224;
+ auto alpha3A_112_123_213 = 8.0*PI_*rhoN*POW3(4*PI_*epsilon_0)/75.0*summerA_112_123_213;
+ auto alpha3A_123_123_224 = 8.0*PI_*rhoN*POW3(4*PI_*epsilon_0)/35.0*summerA_123_123_224;
+ auto alpha3A_224_224_224 = 144.0*PI_*rhoN*POW3(4*PI_*epsilon_0)/245.0*summerA_224_224_224;
+
+ auto alpha3A = 3.0*alpha3A_112_112_224 + 6.0*alpha3A_112_123_213 + 6.0*alpha3A_123_123_224 + alpha3A_224_224_224;
+
+ auto alpha3B_112_112_112 = 32.0*POW3(PI_)*POW2(rhoN)*POW3(4*PI_*epsilon_0)/135.0*sqrt(14*PI_/5.0)*summerB_112_112_112;
+ auto alpha3B_112_123_123 = 64.0*POW3(PI_)*POW2(rhoN)*POW3(4*PI_*epsilon_0)/315.0*sqrt(3.0*PI_)*summerB_112_123_123;
+ auto alpha3B_123_123_224 = -32.0*POW3(PI_)*POW2(rhoN)*POW3(4*PI_*epsilon_0)/45.0*sqrt(22.0*PI_/63.0)*summerB_123_123_224;
+ auto alpha3B_224_224_224 = 32.0*POW3(PI_)*POW2(rhoN)*POW3(4*PI_*epsilon_0)/2025.0*sqrt(2002.0*PI_)*summerB_224_224_224;
+
+ auto alpha3B = alpha3B_112_112_112 + 3.0*alpha3B_112_123_123 + 3.0*alpha3B_123_123_224 + alpha3B_224_224_224;
+
+ return forceeval(alpha3A + alpha3B);
+ }
+
+ /***
+ * \brief Get the contribution to \f$ \alpha = A/(NkT) \f$
+ */
+ template<typename TTYPE, typename RhoType, typename VecType>
+ auto eval(const TTYPE& T, const RhoType& rhoN, const VecType& mole_fractions) const {
+
+ // Calculate the effective reduced diameter (cubed) to be used for evaluation
+ // Eq. 24 from Gubbins
+ std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])> sigma_x3 = 0.0;
+ auto N = mole_fractions.size();
+ for (auto i = 0; i < N; ++i){
+ for (auto j = 0; j < N; ++j){
+ auto sigmaij = (sigma_m[i] + sigma_m[j])/2;
+ sigma_x3 += mole_fractions[i]*mole_fractions[j]*POW3(sigmaij);
+ }
+ }
+ auto rhostar = rhoN*sigma_x3;
+
+ auto alpha2 = get_alpha2(T, rhoN, rhostar, mole_fractions);
+ auto alpha3 = get_alpha3(T, rhoN, rhostar, mole_fractions);
+ auto alpha = forceeval(alpha2/(1.0-alpha3/alpha2));
+
+ using alpha2_t = decltype(alpha2);
+ using alpha3_t = decltype(alpha3);
+ using alpha_t = decltype(alpha);
+ struct Terms{
+ alpha2_t alpha2;
+ alpha3_t alpha3;
+ alpha_t alpha;
+ };
+ return Terms{alpha2, alpha3, alpha};
+ }
+};
+using MCGTL = MultipolarContributionGubbinsTwu<LuckasJIntegral, LuckasKIntegral>;
+
+}
+}
diff --git a/include/teqp/models/saftvrmie.hpp b/include/teqp/models/saftvrmie.hpp
index b1f4b26..80b844b 100644
--- a/include/teqp/models/saftvrmie.hpp
+++ b/include/teqp/models/saftvrmie.hpp
@@ -704,12 +704,15 @@ struct SAFTVRMieChainContributionTerms{
\brief A class used to evaluate mixtures using the SAFT-VR-Mie model
*/
class SAFTVRMieMixture {
+public:
+ using PCSAFTDipolarContribution = SAFTpolar::DipolarContributionGrossVrabec;
+ using PCSAFTQuadrupolarContribution = SAFTpolar::QuadrupolarContributionGrossVrabec;
private:
std::vector<std::string> names;
const SAFTVRMieChainContributionTerms terms;
- std::optional<PCSAFT::PCSAFTDipolarContribution> dipolar; // Can be present or not
- std::optional<PCSAFT::PCSAFTQuadrupolarContribution> quadrupolar; // Can be present or not
+ std::optional<PCSAFTDipolarContribution> dipolar; // Can be present or not
+ std::optional<PCSAFTQuadrupolarContribution> quadrupolar; // Can be present or not
void check_kmat(const Eigen::ArrayXXd& kmat, std::size_t N) {
if (kmat.size() == 0){
@@ -750,7 +753,7 @@ private:
}
}
- auto build_dipolar(const std::vector<SAFTVRMieCoeffs> &coeffs) -> std::optional<PCSAFT::PCSAFTDipolarContribution>{
+ auto build_dipolar(const std::vector<SAFTVRMieCoeffs> &coeffs) -> std::optional<PCSAFTDipolarContribution>{
Eigen::ArrayXd mustar2(coeffs.size()), nmu(coeffs.size());
auto i = 0;
for (const auto &coeff : coeffs) {
@@ -762,9 +765,9 @@ private:
return std::nullopt; // No dipolar contribution is present
}
// The dispersive and hard chain initialization has already happened at this point
- return PCSAFT::PCSAFTDipolarContribution(terms.m, terms.sigma_A, terms.epsilon_over_k, mustar2, nmu);
+ return PCSAFTDipolarContribution(terms.m, terms.sigma_A, terms.epsilon_over_k, mustar2, nmu);
}
- auto build_quadrupolar(const std::vector<SAFTVRMieCoeffs> &coeffs) -> std::optional<PCSAFT::PCSAFTQuadrupolarContribution>{
+ auto build_quadrupolar(const std::vector<SAFTVRMieCoeffs> &coeffs) -> std::optional<PCSAFTQuadrupolarContribution>{
// The dispersive and hard chain initialization has already happened at this point
Eigen::ArrayXd Qstar2(coeffs.size()), nQ(coeffs.size());
auto i = 0;
@@ -776,7 +779,7 @@ private:
if ((Qstar2*nQ).cwiseAbs().sum() == 0){
return std::nullopt; // No quadrupolar contribution is present
}
- return PCSAFT::PCSAFTQuadrupolarContribution(terms.m, terms.sigma_A, terms.epsilon_over_k, Qstar2, nQ);
+ return PCSAFTQuadrupolarContribution(terms.m, terms.sigma_A, terms.epsilon_over_k, Qstar2, nQ);
}
public:
SAFTVRMieMixture(const std::vector<std::string> &names, const std::optional<Eigen::ArrayXXd>& kmat = std::nullopt) : SAFTVRMieMixture(get_coeffs_from_names(names), kmat){};
diff --git a/src/tests/catch_test_SAFTpolar.cxx b/src/tests/catch_test_SAFTpolar.cxx
index 966012a..9c91010 100644
--- a/src/tests/catch_test_SAFTpolar.cxx
+++ b/src/tests/catch_test_SAFTpolar.cxx
@@ -4,20 +4,57 @@
using Catch::Approx;
#include "teqp/models/saft/correlation_integrals.hpp"
+#include "teqp/models/saft/polar_terms.hpp"
+#include "teqp/types.hpp"
+using namespace teqp;
using namespace teqp::SAFTpolar;
TEST_CASE("Evaluation of J^{(n)}", "[LuckasJn]")
{
- LuckasJIntegral<12> J12;
+ LuckasJIntegral J12{12};
auto Jval = J12.get_J(3.0, 1.0);
}
TEST_CASE("Evaluation of K(xxx, yyy)", "[LuckasKnn]")
{
- auto Kval23 = LuckasKIntegral<222, 333>().get_K(1.0, 0.9);
+ auto Kval23 = LuckasKIntegral(222, 333).get_K(1.0, 0.9);
CHECK(Kval23 == Approx(0.03332).margin(0.02));
- auto Kval45 = LuckasKIntegral<444, 555>().get_K(1.0, 0.9);
+ auto Kval45 = LuckasKIntegral(444, 555).get_K(1.0, 0.9);
CHECK(Kval45 == Approx(0.01541).margin(0.02));
}
+
+TEST_CASE("Evaluation of Gubbins and Twu", "[GTLPolar]")
+{
+ auto sigma_m = (Eigen::ArrayXd(2) << 3e-10, 3.1e-10).finished();
+ auto epsilon_over_k= (Eigen::ArrayXd(2) << 200, 300).finished();
+ auto mubar2 = (Eigen::ArrayXd(2) << 0.0, 0.5).finished();
+ auto Qbar2 = (Eigen::ArrayXd(2) << 0.5, 0).finished();
+ MCGTL GTL{sigma_m, epsilon_over_k, mubar2, Qbar2};
+ auto z = (Eigen::ArrayXd(2) << 0.1, 0.9).finished();
+ auto rhoN = std::complex<double>(300, 1e-100);
+ GTL.eval(300.0, rhoN, z);
+}
+
+// This test is used to make sure that replacing std::abs with a more flexible function
+// that can handle differentation types like std::complex<double> is still ok
+
+TEST_CASE("Check derivative of |x|", "[diffabs]")
+{
+ double h = 1e-100;
+ auto safe_abs1 = [](const auto&x) { return sqrt(x*x); };
+ SECTION("|x|"){
+ CHECK(safe_abs1(3.1) == 3.1);
+ CHECK(safe_abs1(-3.1) == 3.1);
+ }
+ SECTION("sqrt(x^2)"){
+ CHECK(safe_abs1(std::complex<double>(3.1, h)).imag()/h == 1);
+ CHECK(safe_abs1(std::complex<double>(-3.1, h)).imag()/h == -1);
+ }
+ auto safe_abs2 = [](const auto&x) { return (getbaseval(x) < 0) ? -x : x; };
+ SECTION("check base"){
+ CHECK(safe_abs2(std::complex<double>(3.1, h)).imag()/h == 1);
+ CHECK(safe_abs2(std::complex<double>(-3.1, h)).imag()/h == -1);
+ }
+}
--
GitLab