Newer
Older
#pragma once
/// Coefficients for one fluid
struct SAFTCoeffs {
std::string name; ///< Name of fluid
double m, ///< number of segments
sigma_Angstrom, ///< [A] segment diameter
epsilon_over_k; ///< [K] depth of pair potential divided by Boltzman constant
std::string BibTeXKey; ///< The BibTeXKey for the reference for these coefficients
};
/// Manager class for PCSAFT coefficients
class PCSAFTLibrary {
std::map<std::string, SAFTCoeffs> coeffs;
public:
PCSAFTLibrary() {
insert_normal_fluid("Methane", 1.0000, 3.7039, 150.03, "Gross-IECR-2001");
insert_normal_fluid("Ethane", 1.6069, 3.5206, 191.42, "Gross-IECR-2001");
}
void insert_normal_fluid(const std::string& name, double m, const double sigma_Angstrom, const double epsilon_over_k, const std::string& BibTeXKey) {
SAFTCoeffs coeff;
coeff.name = name;
coeff.m = m;
coeff.sigma_Angstrom = sigma_Angstrom;
coeff.epsilon_over_k = epsilon_over_k;
coeff.BibTeXKey = BibTeXKey;
coeffs.insert(std::pair<std::string, SAFTCoeffs>(name, coeff));
}
const auto& get_normal_fluid(const std::string& name) {
auto it = coeffs.find(name);
if (it != coeffs.end()) {
return it->second;
}
else {
throw std::invalid_argument("Bad name:" + name);
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
}
}
};
/// Eqn. A.11
/// Erratum: should actually be 1/RHS of equation A.11 according to sample
/// FORTRAN code
template <typename Eta, typename Mbar>
auto C1(const Eta& eta, Mbar mbar) {
return 1.0 / (1.0
+ mbar * (8.0 * eta - 2.0 * eta * eta) / pow(1.0 - eta, 4)
+ (1.0 - mbar) * (20.0 * eta - 27.0 * eta * eta + 12.0 * pow(eta, 3) - 2.0 * pow(eta, 4)) / pow((1.0 - eta) * (2.0 - eta), 2));
}
/// Eqn. A.31
template <typename Eta, typename Mbar>
auto C2(const Eta& eta, Mbar mbar) {
return -pow(C1(eta, mbar), 2) * (
mbar * (-4.0 * eta * eta + 20.0 * eta + 8.0) / pow(1.0 - eta, 5)
+ (1.0 - mbar) * (2.0 * eta * eta * eta + 12.0 * eta * eta - 48.0 * eta + 40.0) / pow((1.0 - eta) * (2.0 - eta), 3)
);
}
/// Eqn. A.18
template<typename TYPE>
auto get_a(TYPE mbar) {
static Eigen::ArrayXd a_0 = (Eigen::ArrayXd(7) << 0.9105631445, 0.6361281449, 2.6861347891, -26.547362491, 97.759208784, -159.59154087, 91.297774084).finished();
static Eigen::ArrayXd a_1 = (Eigen::ArrayXd(7) << -0.3084016918, 0.1860531159, -2.5030047259, 21.419793629, -65.255885330, 83.318680481, -33.746922930).finished();
static Eigen::ArrayXd a_2 = (Eigen::ArrayXd(7) << -0.0906148351, 0.4527842806, 0.5962700728, -1.7241829131, -4.1302112531, 13.776631870, -8.6728470368).finished();
return forceeval(a_0 + ((mbar - 1.0) / mbar) * a_1 + ((mbar - 1.0) / mbar * (mbar - 2.0) / mbar) * a_2).eval();
}
/// Eqn. A.19
template<typename TYPE>
auto get_b(TYPE mbar) {
// See https://stackoverflow.com/a/35170514/1360263
static Eigen::ArrayXd b_0 = (Eigen::ArrayXd(7) << 0.724094694, 2.2382791861, -4.0025849485, -21.003576815, 26.855641363, 206.55133841, -355.60235612).finished();
static Eigen::ArrayXd b_1 = (Eigen::ArrayXd(7) << -0.5755498075, 0.6995095521, 3.8925673390, -17.215471648, 192.67226447, -161.82646165, -165.20769346).finished();
static Eigen::ArrayXd b_2 = (Eigen::ArrayXd(7) << 0.0976883116, -0.2557574982, -9.1558561530, 20.642075974, -38.804430052, 93.626774077, -29.666905585).finished();
return forceeval(b_0 + (mbar - 1.0) / mbar * b_1 + (mbar - 1.0) / mbar * (mbar - 2.0) / mbar * b_2).eval();
}
/// Residual contribution to alphar from hard-sphere (Eqn. A.6)
template<typename VecType>
auto get_alphar_hs(const VecType& zeta) {
auto Upsilon = 1.0 - zeta[3];
return forceeval(1.0 / zeta[0] * (3.0 * zeta[1] * zeta[2] / Upsilon
+ zeta[2] * zeta[2] * zeta[2] / zeta[3] / Upsilon / Upsilon
+ (zeta[2] * zeta[2] * zeta[2] / (zeta[3] * zeta[3]) - zeta[0]) * log(1.0 - zeta[3])
));
}
/// Residual contribution from hard-sphere (Eqn. A.26)
template<typename VecType>
auto Z_hs(const VecType& zeta) {
auto Upsilon = 1.0 - zeta[3];
return forceeval(zeta[3] / Upsilon
+ 3.0 * zeta[1] * zeta[2] / (zeta[0] * pow(Upsilon, 2))
+ (3.0 * pow(zeta[2], 3) - zeta[3] * pow(zeta[2], 3)) / (zeta[0] * pow(Upsilon, 3)));
}
/// Derivative term from Eqn. A.27
template<typename zVecType, typename dVecType>
auto rho_A3_dgij_HS_drhoA3(const zVecType& zeta, const dVecType& d,
std::size_t i, std::size_t j) {
auto Upsilon = 1.0 - zeta[3];
return forceeval(zeta[3] / pow(Upsilon, 2)
+ d[i] * d[j] / (d[i] + d[j]) * (3.0 * zeta[2] / pow(Upsilon, 2) + 6.0 * zeta[2] * zeta[3] / pow(Upsilon, 3))
+ pow(d[i] * d[j] / (d[i] + d[j]), 2) * (4.0 * pow(zeta[2], 2) / pow(Upsilon, 3) + 6.0 * pow(zeta[2], 2) * zeta[3] / pow(Upsilon, 4)));
}
/// Term from Eqn. A.7
template<typename zVecType, typename dVecType>
auto gij_HS(const zVecType& zeta, const dVecType& d,
std::size_t i, std::size_t j) {
auto Upsilon = 1.0 - zeta[3];
return forceeval(1.0 / (Upsilon)+d[i] * d[j] / (d[i] + d[j]) * 3.0 * zeta[2] / pow(Upsilon, 2)
+ pow(d[i] * d[j] / (d[i] + d[j]), 2) * 2.0 * pow(zeta[2], 2) / pow(Upsilon, 3));
}
/// Eqn. A.16, Eqn. A.29
template <typename Eta, typename MbarType>
auto get_I1(const Eta& eta, MbarType mbar) {
auto avec = get_a(mbar);
Eta summer_I1 = 0.0, summer_etadI1deta = 0.0;
for (std::size_t i = 0; i < 7; ++i) {
auto increment = avec(i) * pow(eta, static_cast<int>(i));
summer_I1 = summer_I1 + increment;
summer_etadI1deta = summer_etadI1deta + increment * (i + 1.0);
}
return std::make_tuple(forceeval(summer_I1), forceeval(summer_etadI1deta));
}
/// Eqn. A.17, Eqn. A.30
template <typename Eta, typename MbarType>
auto get_I2(const Eta& eta, MbarType mbar) {
auto bvec = get_b(mbar);
Eta summer_I2 = 0.0 * eta, summer_etadI2deta = 0.0 * eta;
for (std::size_t i = 0; i < 7; ++i) {
auto increment = bvec(i) * pow(eta, static_cast<int>(i));
summer_I2 = summer_I2 + increment;
summer_etadI2deta = summer_etadI2deta + increment * (i + 1.0);
}
return std::make_tuple(forceeval(summer_I2), forceeval(summer_etadI2deta));
}
PCSAFTLibrary library;
/**
Sum up three array-like objects that can each have different container types and value types
template<typename VecType1, typename NType>
auto powvec(const VecType1& v1, NType n) {
auto o = v1;
for (auto i = 0; i < v1.size(); ++i) {
}
/**
Sum up the coefficient-wise product of three array-like objects that can each have different container types and value types
*/
template<typename VecType1, typename VecType2, typename VecType3>
auto sumproduct(const VecType1& v1, const VecType2& v2, const VecType3& v3) {
using ResultType = typename std::common_type_t<decltype(v1[0]), decltype(v2[0]), decltype(v3[0])>;
return forceeval((v1.template cast<ResultType>() * v2.template cast<ResultType>() * v3.template cast<ResultType>()).sum());
}
/// Parameters for model evaluation
class SAFTCalc {
public:
// Just temperature dependent things
Eigen::ArrayX<NumType> d;
// These things also have composition dependence
ProductType m2_epsilon_sigma3_bar, ///< Eq. A. 12
m2_epsilon2_sigma3_bar; ///< Eq. A. 13
};
/// A class used to evaluate mixtures using PC-SAFT model
class PCSAFTMixture {
private:
Eigen::ArrayX<double> m, ///< number of segments
mminus1, ///< m-1
sigma_Angstrom, ///<
epsilon_over_k; ///< depth of pair potential divided by Boltzman constant
std::vector<std::string> names;
double k_ij; ///< binary interaction parameter
public:
PCSAFTMixture(const std::vector<std::string> &names) : names(names)
{
m.resize(names.size());
mminus1.resize(names.size());
sigma_Angstrom.resize(names.size());
epsilon_over_k.resize(names.size());
auto i = 0;
for (auto name : names) {
const SAFTCoeffs& coeff = library.get_normal_fluid(name);
m[i] = coeff.m;
mminus1[i] = m[i] - 1;
sigma_Angstrom[i] = coeff.sigma_Angstrom;
epsilon_over_k[i] = coeff.epsilon_over_k;
i++;
}
k_ij = 0;
};
PCSAFTMixture(const std::vector<SAFTCoeffs> &coeffs)
{
m.resize(coeffs.size());
mminus1.resize(coeffs.size());
sigma_Angstrom.resize(coeffs.size());
epsilon_over_k.resize(coeffs.size());
names.resize(coeffs.size());
auto i = 0;
for (const auto &coeff : coeffs) {
m[i] = coeff.m;
mminus1[i] = m[i] - 1;
sigma_Angstrom[i] = coeff.sigma_Angstrom;
epsilon_over_k[i] = coeff.epsilon_over_k;
names[i] = coeff.name;
i++;
}
k_ij = 0;
};
auto get_m(){ return m; }
auto get_sigma_Angstrom() { return sigma_Angstrom; }
auto get_epsilon_over_k_K() { return epsilon_over_k; }
void print_info() {
std::cout << "i m sigma / A e/kB / K \n ++++++++++++++" << std::endl;
for (auto i = 0; i < m.size(); ++i) {
std::cout << i << " " << m[i] << " " << sigma_Angstrom[i] << " " << epsilon_over_k[i] << std::endl;
}
}
template<typename VecType>
auto max_rhoN(double T, const VecType& mole_fractions) {
auto N = mole_fractions.size();
for (auto i = 0; i < N; ++i) {
d[i] = sigma_Angstrom[i] * (1.0 - 0.12 * exp(-3.0 * epsilon_over_k[i] / T));
}
return 6 * 0.74 / EIGEN_PI / (mole_fractions*m*powvec(d, 3)).sum()*1e30; // particles/m^3
}
const double R = get_R_gas<double>();
template<typename TTYPE, typename RhoType, typename VecType>
auto alphar(const TTYPE& T, const RhoType& rhomolar, const VecType& mole_fractions) const {
std::size_t N = m.size();
if (mole_fractions.size() != N) {
throw std::invalid_argument("Length of mole_fractions (" + std::to_string(mole_fractions.size()) + ") is not the length of components (" + std::to_string(N) + ")");
}
using TRHOType = std::common_type_t<decltype(T), decltype(mole_fractions[0]), decltype(m[0])>;
c.m2_epsilon_sigma3_bar = 0;
c.m2_epsilon2_sigma3_bar = 0;
c.d = sigma_Angstrom*(1.0 - 0.12*exp(-3.0*epsilon_over_k/T)); // [A]
for (std::size_t i = 0; i < N; ++i) {
for (std::size_t j = 0; j < N; ++j) {
// Eq. A.5
auto sigma_ij = 0.5 * sigma_Angstrom[i] + 0.5 * sigma_Angstrom[j];
auto eij_over_k = sqrt(epsilon_over_k[i] * epsilon_over_k[j]) * (1.0 - k_ij);
c.m2_epsilon_sigma3_bar = c.m2_epsilon_sigma3_bar + mole_fractions[i] * mole_fractions[j] * m[i] * m[j] * eij_over_k / T * pow(sigma_ij, 3);
c.m2_epsilon2_sigma3_bar = c.m2_epsilon2_sigma3_bar + mole_fractions[i] * mole_fractions[j] * m[i] * m[j] * pow(eij_over_k / T, 2) * pow(sigma_ij, 3);
}
}
/// Convert from molar density to number density in molecules/Angstrom^3
RhoType rho_A3 = rhomolar * N_A * 1e-30; //[molecules (not moles)/A^3]
constexpr double MY_PI = EIGEN_PI;
/// Evaluate the components of zeta
using ta = std::common_type_t<decltype(pi6), decltype(m[0]), decltype(c.d[0]), decltype(rho_A3)>;
std::vector<ta> zeta(4);
for (std::size_t n = 0; n < 4; ++n) {
// Eqn A.8
TRHOType xmdn = forceeval((mole_fractions*m*dn).sum());
}
/// Packing fraction is the 4-th value in zeta, at index 3
auto [I1, etadI1deta] = get_I1(eta, mbar);
auto [I2, etadI2deta] = get_I2(eta, mbar);
using tt = std::common_type_t<decltype(zeta[0]), decltype(c.d[0])>;
Eigen::ArrayX<tt> lngii_hs(mole_fractions.size());
for (std::size_t i = 0; i < lngii_hs.size(); ++i) {
auto alphar_hc = mbar * get_alphar_hs(zeta) - sumproduct(mole_fractions, mminus1, lngii_hs); // Eq. A.4
auto alphar_disp = -2 * MY_PI * rho_A3 * I1 * c.m2_epsilon_sigma3_bar - MY_PI * rho_A3 * mbar * C1(eta, mbar) * I2 * c.m2_epsilon2_sigma3_bar;
return forceeval(alphar_hc + alphar_disp);