Newer
Older
template<typename X> auto POW2(X x) { return x * x; };
template<typename X> auto POW3(X x) { return x * POW2(x); };
enum class association_classes {not_set, a1A, a2B, a3B, a4C, not_associating};
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);
}
}
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
enum class radial_dist { CS, KG, OT };
/// Function that calculates the association binding strength between site A of molecule i and site B on molecule j
template<typename BType, typename TType, typename RhoType, typename VecType>
auto get_DeltaAB_pure(radial_dist dist, double epsABi, double betaABi, BType b_cubic, TType T, RhoType rhomolar, const VecType& molefrac) {
using eta_type = std::common_type_t<decltype(rhomolar), decltype(b_cubic)>;
eta_type eta;
eta_type g_vm_ref;
// Calculate the contact value of the radial distribution function g(v)
switch (dist) {
case radial_dist::CS: {
// Carnahan - Starling EOS, given by Kontogeorgis et al., Ind.Eng.Chem.Res. 2006, 45, 4855 - 4868, Eq. 4a and 4b:
eta = (rhomolar / 4.0) * b_cubic;
g_vm_ref = (2.0 - eta) / (2.0 * POW3(1.0 - eta));
break;
}
case radial_dist::KG: {
// Function proposed by Kontogeorgis, G.M.; Yakoumis, I.V.; Meijer, H.; Hendriks, E.M.; Moorwood, T., Fluid Phase Equilib. 1999, 158 - 160, 201.
eta = (rhomolar / 4.0) * b_cubic;
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;
}
default: {
throw std::invalid_argument("Bad radial_dist");
}
}
double R_gas = 8.3144598;
// Calculate the association strength between site Ai and Bi for a pure compent
auto DeltaAiBj = forceeval(g_vm_ref*(exp(epsABi /(T*R_gas)) - 1.0)*b_cubic* betaABi);
return DeltaAiBj;
};
/// Routine that calculates the fractions of sites Ai not bound to other active sites for pure fluids
/// Some association schemes are explicitly solvable for self - associating compounds, see Huang and Radosz, Ind.Eng.Chem.Res., 29 (11), 1990
/// So far implemented association schemes : 1A, 2B, 3B, 4C (see Kontogeorgis et al., Ind. Eng. Chem. Res. 2006, 45, 4855 - 4868)
///
template<typename BType, typename TType, typename RhoType, typename VecType>
auto XA_calc_pure(int N_sites, association_classes scheme, double epsABi, double betaABi, const BType b_cubic, const TType T, const RhoType rhomolar, const VecType& molefrac) {
// Matrix XA(A, j) that contains all of the fractions of sites A not bonded to other active sites for each molecule i
// Start values for the iteration(set all sites to non - bonded, = 1)
Eigen::Array<RhoType, Eigen::Dynamic, Eigen::Dynamic> XA; // A maximum of 4 association sites(A, B, C, D)
XA.resize(N_sites, 1);
XA.setOnes();
// Get the association strength between the associating sites
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
auto DeltaAiBj = get_DeltaAB_pure(dist, epsABi, betaABi, b_cubic, T, rhomolar, molefrac);
if (scheme == association_classes::a1A) { // Acids
// Only one association site "A" (OH - group with C = O - group)
XA(0, 0) = forceeval((-1.0 + sqrt(1.0 + 4.0 * rhomolar * DeltaAiBj)) / (2.0 * rhomolar * DeltaAiBj));
}
else if (scheme == association_classes::a2B) { // Alcohols
// Two association sites "A" and "B"
XA(0, 0) = forceeval((-1.0 + sqrt(1.0 + 4.0 * rhomolar * DeltaAiBj)) / (2.0 * rhomolar * DeltaAiBj));
XA(1, 0) = XA(0, 0); // XB = XA;
}
else if (scheme == association_classes::a3B) { // Glycols
// Three association sites "A", "B", "C"
XA(0, 0) = forceeval((-(1.0 - rhomolar * DeltaAiBj) + sqrt(POW2(1.0 + rhomolar * DeltaAiBj) + 4.0 * rhomolar * DeltaAiBj)) / (4.0 * rhomolar * DeltaAiBj));
XA(1, 0) = XA(0, 0); // XB = XA
XA(2, 0) = 2 * XA(0, 0) - 1; // XC = 2XA - 1
}
else if (scheme == association_classes::a4C) { // Water
// Four association sites "A", "B", "C", "D"
XA(0, 0) = forceeval((-1.0 + sqrt(1.0 + 8.0 * rhomolar * DeltaAiBj)) / (4.0 * rhomolar * DeltaAiBj));
XA(1, 0) = XA(0, 0); // XB = XA
XA(2, 0) = XA(0, 0); // XC = XA
XA(3, 0) = XA(0, 0); // XD = XA
}
else if (scheme == association_classes::not_associating) { // non - associating compound
XA(0, 0) = 1;
XA(1, 0) = 1;
XA(2, 0) = 1;
XA(3, 0) = 1;
}
else {
throw std::invalid_argument("Bad scheme");
}
return XA;
};
enum class cubic_flag {not_set, PR, SRK};
auto get_cubic_flag(const std::string& s) {
if (s == "PR") { return cubic_flag::PR; }
else if (s == "SRK") { return cubic_flag::SRK; }
else {
throw std::invalid_argument("bad cubic flag:" + s);
}
}
class CPACubic {
private:
std::valarray<double> a0, bi, c1, Tc;
double delta_1, delta_2;
std::valarray<std::valarray<double>> k_ij;
public:
CPACubic(cubic_flag flag, const std::valarray<double> &a0, const std::valarray<double> &bi, const std::valarray<double> &c1, const std::valarray<double> &Tc) : a0(a0), bi(bi), c1(c1), Tc(Tc) {
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
switch (flag) {
case cubic_flag::PR:
{ delta_1 = 1 + sqrt(2); delta_2 = 1 - sqrt(2); break; }
case cubic_flag::SRK:
{ delta_1 = 0; delta_2 = 1; break; }
default:
throw std::invalid_argument("Bad cubic flag");
}
k_ij.resize(Tc.size()); for (auto i = 0; i < k_ij.size(); ++i) { k_ij[i].resize(Tc.size()); }
};
template<typename TType>
auto get_ai(TType T, int i) const {
return a0[i] * POW2(1 + c1[i]*(1 - sqrt(T / Tc[i])));
}
template<typename TType, typename VecType>
auto get_ab(const TType T, const VecType& molefrac) const {
using return_type = std::common_type_t<decltype(T), decltype(molefrac[0])>;
return_type asummer = 0.0, bsummer = 0.0;
for (auto i = 0; i < molefrac.size(); ++i) {
bsummer += molefrac[i] * bi[i];
auto ai = get_ai(T, i);
for (auto j = 0; j < molefrac.size(); ++j) {
auto aj = get_ai(T, j);
auto a_ij = (1.0 - k_ij[i][j]) * sqrt(ai * aj);
asummer += molefrac[i] * molefrac[j] * a_ij;
}
}
return std::make_tuple(asummer, bsummer);
}
template<typename TType, typename RhoType, typename VecType>
auto alphar(const TType T, const RhoType rhomolar, const VecType& molefrac) const {
auto [a_cubic, b_cubic] = get_ab(T, molefrac);
auto R_gas = 8.3144598;
return forceeval(-log(1 - b_cubic * rhomolar) - a_cubic / R_gas / T * log((delta_1 * b_cubic * rhomolar + 1) / (delta_2 * b_cubic * rhomolar + 1)) / b_cubic / (delta_1 - delta_2));
}
};
template<typename Cubic>
class CPAAssociation {
private:
const std::valarray<double> epsABi, betaABi;
auto get_N_sites(const std::vector<association_classes> &classes) {
auto get_N = [](auto cl) {
switch (cl) {
case association_classes::a1A: return 1;
case association_classes::a2B: return 2;
case association_classes::a3B: return 3;
case association_classes::a4C: return 4;
default: throw std::invalid_argument("Bad association class");
}
};
for (auto cl : classes) {
N_sites_out.push_back(get_N(cl));
}
return N_sites_out;
}
public:
CPAAssociation(const Cubic &&cubic, const std::vector<association_classes>& classes, const std::valarray<double> &epsABi, const std::valarray<double> &betaABi)
: cubic(cubic), classes(classes), epsABi(epsABi), betaABi(betaABi), N_sites(get_N_sites(classes)) {};
template<typename TType, typename RhoType, typename VecType>
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 the fraction of sites not bonded with other active sites
auto XA = XA_calc_pure(N_sites[0], classes[0], epsABi[0], betaABi[0], b_cubic, T, rhomolar, molefrac);
using return_type = std::common_type_t<decltype(T), decltype(rhomolar), decltype(molefrac[0])>;
return_type alpha_r_asso = 0.0;
auto i = 0;
for (auto xi : molefrac){ // loop over all components
auto XAi = XA.col(i);
alpha_r_asso += forceeval(xi * (log(XAi) - XAi / 2).sum());
alpha_r_asso += xi*static_cast<double>(N_sites[i])/2;
}
};
template <typename Cubic, typename Assoc>
class CPA {
public:
CPA(Cubic &&cubic, Assoc &&assoc) : cubic(cubic), assoc(assoc) {
}
/// Residual dimensionless Helmholtz energy from the SRK or PR core and contribution due to association
/// alphar = a/(R*T) where a and R are both molar quantities
template<typename TType, typename RhoType, typename VecType>
auto alphar(const TType& T, const RhoType& rhomolar, const VecType& molefrac) const {
// Calculate the contribution to alphar from the conventional cubic EOS
auto alpha_r_cubic = cubic.alphar(T, rhomolar, molefrac);
// Calculate the contribution to alphar from association
auto alpha_r_assoc = assoc.alphar(T, rhomolar, molefrac);
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
};
/// A factory function to return an instantiated CPA instance given
/// the JSON representation of the model
auto CPAfactory(const nlohmann::json &j){
auto build_cubic = [](const auto& j) {
auto N = j["pures"].size();
std::valarray<double> a0i(N), bi(N), c1(N), Tc(N);
std::size_t i = 0;
for (auto p : j["pures"]) {
a0i[i] = p["a0i / Pa m^6/mol^2"];
bi[i] = p["bi / m^3/mol"];
c1[i] = p["c1"];
Tc[i] = p["Tc / K"];
i++;
}
return CPACubic(get_cubic_flag(j["cubic"]), a0i, bi, c1, Tc);
};
auto build_assoc = [](const auto &&cubic, const auto& j) {
auto N = j["pures"].size();
std::vector<association_classes> classes;
std::valarray<double> 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++;
}
return CPAAssociation(std::move(cubic), classes, epsABi, betaABi);
};
return CPA(build_cubic(j), build_assoc(build_cubic(j), j));
}
}; /* namespace CPA */