CoolProp 6.8.1dev
An open-source fluid property and humid air property database
GeneralizedCubic.cpp
Go to the documentation of this file.
1#include "GeneralizedCubic.h"
2#include "CPnumerics.h"
3#include <cmath>
4
5double BasicMathiasCopemanAlphaFunction::term(double tau, std::size_t itau) {
6
7 // If we are not using the full Mathias-Copeman formulation for a_ii,
8 // we just use the simple results from the supplemental information because
9 // they are much more computationally efficient
10
11 // All derivatives have a common bracketed term, so we factor it out
12 // and calculate it here
13 double B = 1 + m * (1 - sqrt_Tr_Tci * sqrt(1 / tau));
14
15 switch (itau) {
16 case 0:
17 return a0 * B * B;
18 case 1:
19 return a0 * m * B / pow(tau, 3.0 / 2.0) * sqrt_Tr_Tci;
20 case 2:
21 return a0 * m / 2.0 * (m / pow(tau, 3) * Tr_over_Tci - 3 * B / pow(tau, 5.0 / 2.0) * sqrt_Tr_Tci);
22 case 3:
23 return (3.0 / 4.0) * a0 * m * (-3.0 * m / pow(tau, 4) * Tr_over_Tci + 5 * B / pow(tau, 7.0 / 2.0) * sqrt_Tr_Tci);
24 case 4:
25 return (3.0 / 8.0) * a0 * m * (29.0 * m / pow(tau, 5) * Tr_over_Tci - 35 * B / pow(tau, 9.0 / 2.0) * sqrt_Tr_Tci);
26 default:
27 throw -1;
28 }
29}
30
31double MathiasCopemanAlphaFunction::term(double tau, std::size_t itau) {
32 // Here we are using the full Mathias-Copeman formulation, introducing
33 // some additional computational effort, so we only evaluate the parameters that
34 // we actually need to evaluate, otherwise we just set their value to zero
35 // See info on the conditional (ternary) operator : http://www.cplusplus.com/articles/1AUq5Di1/
36 // Furthermore, this should help with branch prediction
37 double Di = 1 - sqrt_Tr_Tci / sqrt(tau);
38 double dDi_dtau = (itau >= 1) ? (1.0 / 2.0) * sqrt_Tr_Tci / (pow(tau, 1.5)) : 0;
39 double d2Di_dtau2 = (itau >= 2) ? -(3.0 / 4.0) * sqrt_Tr_Tci / (pow(tau, 2.5)) : 0;
40 double d3Di_dtau3 = (itau >= 3) ? (15.0 / 8.0) * sqrt_Tr_Tci / (pow(tau, 3.5)) : 0;
41 double d4Di_dtau4 = (itau >= 4) ? -(105.0 / 16.0) * sqrt_Tr_Tci / (pow(tau, 4.5)) : 0;
42
43 double Bi = 1, dBi_dtau = 0, d2Bi_dtau2 = 0, d3Bi_dtau3 = 0, d4Bi_dtau4 = 0;
44 for (int n = 1; n <= 3; ++n) {
45 Bi += c[n - 1] * pow(Di, n);
46 dBi_dtau += (itau < 1) ? 0 : (n * c[n - 1] * pow(Di, n - 1) * dDi_dtau);
47 d2Bi_dtau2 += (itau < 2) ? 0 : n * c[n - 1] * ((n - 1) * pow(dDi_dtau, 2) + Di * d2Di_dtau2) * pow(Di, n - 2);
48 d3Bi_dtau3 += (itau < 3)
49 ? 0
50 : n * c[n - 1] * (3 * (n - 1) * Di * dDi_dtau * d2Di_dtau2 + (n * n - 3 * n + 2) * pow(dDi_dtau, 3) + pow(Di, 2) * d3Di_dtau3)
51 * pow(Di, n - 3);
52 d4Bi_dtau4 +=
53 (itau < 4)
54 ? 0
55 : n * c[n - 1]
56 * (6 * (n * n - 3 * n + 2) * Di * pow(dDi_dtau, 2) * d2Di_dtau2 + (n * n * n - 6 * n * n + 11 * n - 6) * pow(dDi_dtau, 4)
57 + (4 * n * dDi_dtau * d3Di_dtau3 + 3 * n * pow(d2Di_dtau2, 2) - 4 * dDi_dtau * d3Di_dtau3 - 3 * pow(d2Di_dtau2, 2)) * pow(Di, 2)
58 + pow(Di, 3) * d4Di_dtau4)
59 * pow(Di, n - 4);
60 }
61 switch (itau) {
62 case 0:
63 return a0 * Bi * Bi;
64 case 1:
65 return 2 * a0 * Bi * dBi_dtau;
66 case 2:
67 return 2 * a0 * (Bi * d2Bi_dtau2 + dBi_dtau * dBi_dtau);
68 case 3:
69 return 2 * a0 * (Bi * d3Bi_dtau3 + 3 * dBi_dtau * d2Bi_dtau2);
70 case 4:
71 return 2 * a0 * (Bi * d4Bi_dtau4 + 4 * dBi_dtau * d3Bi_dtau3 + 3 * pow(d2Bi_dtau2, 2));
72 default:
73 throw -1;
74 }
75}
76
77double TwuAlphaFunction::term(double tau, std::size_t itau) {
78 // Here we are using the Twu formulation, introducing
79 // some additional computational effort, so we only evaluate the parameters that
80 // we actually need to evaluate, otherwise we just set their value to zero
81 // See info on the conditional (ternary) operator : http://www.cplusplus.com/articles/1AUq5Di1/
82 // Furthermore, this should help with branch prediction
83
84 const double L = c[0], M = c[1], N = c[2];
85 double A = pow(Tr_over_Tci / tau, M * N);
86 double B1 = (itau < 1) ? 0 : N / tau * (L * M * A - M + 1);
87 double dB1_dtau = (itau < 2) ? 0 : N / powInt(tau, 2) * (-L * M * M * N * A - L * M * A + M - 1);
88 double d2B1_dtau2 = (itau < 3) ? 0 : N / powInt(tau, 3) * (L * M * M * M * N * N * A + 3 * L * M * M * N * A + 2 * L * M * A - 2 * M + 2);
89 double d3B1_dtau3 =
90 (itau < 4) ? 0
91 : -N / powInt(tau, 4)
92 * (L * powInt(M, 4) * powInt(N, 3) * A + 6 * L * M * M * M * N * N * A + 11 * L * M * M * N * A + 6 * L * M * A - 6 * M + 6);
93
94 double dam_dtau, d2am_dtau2, d3am_dtau3, d4am_dtau4;
95 double am = a0 * pow(Tr_over_Tci / tau, N * (M - 1)) * exp(L * (1 - A));
96
97 if (itau == 0) {
98 return am;
99 } else {
100 // Calculate terms as needed
101 dam_dtau = am * B1;
102 d2am_dtau2 = (itau < 2) ? 0 : B1 * dam_dtau + am * dB1_dtau;
103 d3am_dtau3 = (itau < 3) ? 0 : B1 * d2am_dtau2 + am * d2B1_dtau2 + 2 * dB1_dtau * dam_dtau;
104 d4am_dtau4 = (itau < 4) ? 0 : B1 * d3am_dtau3 + am * d3B1_dtau3 + 3 * dB1_dtau * d2am_dtau2 + 3 * d2B1_dtau2 * dam_dtau;
105 }
106 switch (itau) {
107 case 1:
108 return dam_dtau;
109 case 2:
110 return d2am_dtau2;
111 case 3:
112 return d3am_dtau3;
113 case 4:
114 return d4am_dtau4;
115 default:
116 throw -1;
117 }
118}
119
120AbstractCubic::AbstractCubic(std::vector<double> Tc, std::vector<double> pc, std::vector<double> acentric, double R_u, double Delta_1, double Delta_2,
121 std::vector<double> C1, std::vector<double> C2, std::vector<double> C3)
122 : Tc(Tc), pc(pc), acentric(acentric), R_u(R_u), Delta_1(Delta_1), Delta_2(Delta_2) {
123 N = static_cast<int>(Tc.size());
124 k.resize(N, std::vector<double>(N, 0));
125 cm = 0.;
126 alpha.resize(N);
127 T_r = 1.0;
128 rho_r = 1.0;
129};
130
131void AbstractCubic::set_alpha(const std::vector<double>& C1, const std::vector<double>& C2, const std::vector<double>& C3) {
133 alpha.resize(Tc.size());
135 if (C1.empty() && C2.empty() && C3.empty()) {
136 for (std::size_t i = 0; i < Tc.size(); ++i) {
137 alpha[i].reset(new BasicMathiasCopemanAlphaFunction(a0_ii(i), m_ii(i), T_r / Tc[i]));
138 }
139 } else {
141 for (std::size_t i = 0; i < Tc.size(); ++i) {
142 alpha[i].reset(new MathiasCopemanAlphaFunction(a0_ii(i), C1[i], C2[i], C3[i], T_r / Tc[i]));
143 }
144 }
145}
146
147double AbstractCubic::am_term(double tau, const std::vector<double>& x, std::size_t itau) {
148 double summer = 0;
149 for (int i = N - 1; i >= 0; --i) {
150 for (int j = N - 1; j >= 0; --j) {
151 summer += x[i] * x[j] * aij_term(tau, i, j, itau);
152 }
153 }
154 return summer;
155}
156double AbstractCubic::d_am_term_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent) {
157 if (xN_independent) {
158 double summer = 0;
159 for (int j = N - 1; j >= 0; --j) {
160 summer += x[j] * aij_term(tau, i, j, itau);
161 }
162 return 2 * summer;
163 } else {
164 double summer = 0;
165 for (int k = N - 2; k >= 0; --k) {
166 summer += x[k] * (aij_term(tau, i, k, itau) - aij_term(tau, k, N - 1, itau));
167 }
168 return 2 * (summer + x[N - 1] * (aij_term(tau, N - 1, i, itau) - aij_term(tau, N - 1, N - 1, itau)));
169 }
170}
171double AbstractCubic::d2_am_term_dxidxj(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j,
172 bool xN_independent) {
173 if (xN_independent) {
174 return 2 * aij_term(tau, i, j, itau);
175 } else {
176 return 2 * (aij_term(tau, i, j, itau) - aij_term(tau, j, N - 1, itau) - aij_term(tau, N - 1, i, itau) + aij_term(tau, N - 1, N - 1, itau));
177 }
178}
179
180double AbstractCubic::d3_am_term_dxidxjdxk(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j, std::size_t k,
181 bool xN_independent) {
182 return 0;
183}
184
185double AbstractCubic::bm_term(const std::vector<double>& x) {
186 double summer = 0;
187 for (int i = N - 1; i >= 0; --i) {
188 summer += x[i] * b0_ii(i);
189 }
190 return summer;
191}
192double AbstractCubic::d_bm_term_dxi(const std::vector<double>& x, std::size_t i, bool xN_independent) {
193 if (xN_independent) {
194 return b0_ii(i);
195 } else {
196 return b0_ii(i) - b0_ii(N - 1);
197 }
198}
199double AbstractCubic::d2_bm_term_dxidxj(const std::vector<double>& x, std::size_t i, std::size_t j, bool xN_independent) {
200 return 0;
201}
202double AbstractCubic::d3_bm_term_dxidxjdxk(const std::vector<double>& x, std::size_t i, std::size_t j, std::size_t k, bool xN_independent) {
203 return 0;
204}
205
206double AbstractCubic::cm_term() {
207 return cm;
208}
209
210double AbstractCubic::aii_term(double tau, std::size_t i, std::size_t itau) {
211 return alpha[i]->term(tau, itau);
212}
213double AbstractCubic::u_term(double tau, std::size_t i, std::size_t j, std::size_t itau) {
214 double aii = aii_term(tau, i, 0), ajj = aii_term(tau, j, 0);
215 switch (itau) {
216 case 0:
217 return aii * ajj;
218 case 1:
219 return aii * aii_term(tau, j, 1) + ajj * aii_term(tau, i, 1);
220 case 2:
221 return (aii * aii_term(tau, j, 2) + 2 * aii_term(tau, i, 1) * aii_term(tau, j, 1) + ajj * aii_term(tau, i, 2));
222 case 3:
223 return (aii * aii_term(tau, j, 3) + 3 * aii_term(tau, i, 1) * aii_term(tau, j, 2) + 3 * aii_term(tau, i, 2) * aii_term(tau, j, 1)
224 + ajj * aii_term(tau, i, 3));
225 case 4:
226 return (aii * aii_term(tau, j, 4) + 4 * aii_term(tau, i, 1) * aii_term(tau, j, 3) + 6 * aii_term(tau, i, 2) * aii_term(tau, j, 2)
227 + 4 * aii_term(tau, i, 3) * aii_term(tau, j, 1) + ajj * aii_term(tau, i, 4));
228 default:
229 throw -1;
230 }
231}
232double AbstractCubic::aij_term(double tau, std::size_t i, std::size_t j, std::size_t itau) {
233 double u = u_term(tau, i, j, 0);
234
235 switch (itau) {
236 case 0:
237 return (1 - k[i][j]) * sqrt(u);
238 case 1:
239 return (1 - k[i][j]) / (2.0 * sqrt(u)) * u_term(tau, i, j, 1);
240 case 2:
241 return (1 - k[i][j]) / (4.0 * pow(u, 3.0 / 2.0)) * (2 * u * u_term(tau, i, j, 2) - pow(u_term(tau, i, j, 1), 2));
242 case 3:
243 return (1 - k[i][j]) / (8.0 * pow(u, 5.0 / 2.0))
244 * (4 * pow(u, 2) * u_term(tau, i, j, 3) - 6 * u * u_term(tau, i, j, 1) * u_term(tau, i, j, 2) + 3 * pow(u_term(tau, i, j, 1), 3));
245 case 4:
246 return (1 - k[i][j]) / (16.0 * pow(u, 7.0 / 2.0))
247 * (-4 * pow(u, 2) * (4 * u_term(tau, i, j, 1) * u_term(tau, i, j, 3) + 3 * pow(u_term(tau, i, j, 2), 2))
248 + 8 * pow(u, 3) * u_term(tau, i, j, 4) + 36 * u * pow(u_term(tau, i, j, 1), 2) * u_term(tau, i, j, 2)
249 - 15 * pow(u_term(tau, i, j, 1), 4));
250 default:
251 throw -1;
252 }
253}
254double AbstractCubic::psi_minus(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta) {
255 if (itau > 0) return 0.0;
256 double bmc = bm_term(x) - cm_term(); // appears only in the form (b-c) in the equations
257 double bracket = 1 - bmc * delta * rho_r;
258
259 switch (idelta) {
260 case 0:
261 return -log(bracket);
262 case 1:
263 return bmc * rho_r / bracket;
264 case 2:
265 return pow(bmc * rho_r / bracket, 2);
266 case 3:
267 return 2 * pow(bmc * rho_r / bracket, 3);
268 case 4:
269 return 6 * pow(bmc * rho_r / bracket, 4);
270 default:
271 throw -1;
272 }
273}
274double AbstractCubic::d_psi_minus_dxi(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
275 bool xN_independent) {
276 if (itau > 0) return 0.0;
277 double bmc = bm_term(x) - cm_term(); // appears only in the form (b-c) in the equations
278 double db_dxi = d_bm_term_dxi(x, i, xN_independent);
279 double bracket = 1 - bmc * delta * rho_r;
280
281 switch (idelta) {
282 case 0:
283 return delta * rho_r * db_dxi / bracket;
284 case 1:
285 return rho_r * db_dxi / pow(bracket, 2);
286 case 2:
287 return 2 * pow(rho_r, 2) * bmc * db_dxi / pow(bracket, 3);
288 case 3:
289 return 6 * pow(rho_r, 3) * pow(bmc, 2) * db_dxi / pow(bracket, 4);
290 case 4:
291 return 24 * pow(rho_r, 4) * pow(bmc, 3) * db_dxi / pow(bracket, 5);
292 default:
293 throw -1;
294 }
295}
296double AbstractCubic::d2_psi_minus_dxidxj(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
297 std::size_t j, bool xN_independent) {
298 if (itau > 0) return 0.0;
299 double bmc = bm_term(x) - cm_term(); // appears only in the form (b-c) in the equations
300 double db_dxi = d_bm_term_dxi(x, i, xN_independent), db_dxj = d_bm_term_dxi(x, j, xN_independent),
301 d2b_dxidxj = d2_bm_term_dxidxj(x, i, j, xN_independent);
302 double bracket = 1 - bmc * delta * rho_r;
303
304 switch (idelta) {
305 case 0:
306 return pow(delta * rho_r, 2) * db_dxi * db_dxj / pow(bracket, 2) + delta * rho_r * d2b_dxidxj / bracket;
307 case 1:
308 return 2 * delta * pow(rho_r, 2) * db_dxi * db_dxj / pow(bracket, 3) + rho_r * d2b_dxidxj / pow(bracket, 2);
309 case 2:
310 return 2 * pow(rho_r, 2) * db_dxi * db_dxj / pow(bracket, 4) * (2 * delta * rho_r * bmc + 1)
311 + 2 * pow(rho_r, 2) * bmc * d2b_dxidxj / pow(bracket, 3);
312 case 3:
313 return 12 * pow(rho_r, 3) * bmc * db_dxi * db_dxj / pow(bracket, 5) * (delta * rho_r * bmc + 1)
314 + 6 * pow(rho_r, 3) * pow(bmc, 2) * d2b_dxidxj / pow(bracket, 4);
315 case 4:
316 return 24 * pow(rho_r, 4) * pow(bmc, 2) * db_dxi * db_dxj / pow(bracket, 6) * (2 * delta * rho_r * bmc + 3)
317 + 24 * pow(rho_r, 4) * pow(bmc, 3) * d2b_dxidxj / pow(bracket, 5);
318 default:
319 throw -1;
320 }
321}
322double AbstractCubic::d3_psi_minus_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
323 std::size_t j, std::size_t k, bool xN_independent) {
324 if (itau > 0) return 0.0;
325 double bmc = bm_term(x) - cm_term(); // appears only in the form (b-c) in the equations
326 double db_dxi = d_bm_term_dxi(x, i, xN_independent), db_dxj = d_bm_term_dxi(x, j, xN_independent), db_dxk = d_bm_term_dxi(x, k, xN_independent),
327 d2b_dxidxj = d2_bm_term_dxidxj(x, i, j, xN_independent), d2b_dxidxk = d2_bm_term_dxidxj(x, i, k, xN_independent),
328 d2b_dxjdxk = d2_bm_term_dxidxj(x, j, k, xN_independent), d3b_dxidxjdxk = d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent);
329 double bracket = 1 - bmc * delta * rho_r;
330
331 switch (idelta) {
332 case 0:
333 return delta * rho_r * d3b_dxidxjdxk / bracket + 2 * pow(delta * rho_r, 3) * db_dxi * db_dxj * db_dxk / pow(bracket, 3)
334 + pow(delta * rho_r, 2) / pow(bracket, 2) * (db_dxi * d2b_dxjdxk + db_dxj * d2b_dxidxk + db_dxk * d2b_dxidxj);
335 case 1:
336 return rho_r * d3b_dxidxjdxk / pow(bracket, 2) + 6 * pow(delta, 2) * pow(rho_r, 3) * db_dxi * db_dxj * db_dxk / pow(bracket, 4)
337 + 2 * delta * pow(rho_r, 2) / pow(bracket, 3) * (db_dxi * d2b_dxjdxk + db_dxj * d2b_dxidxk + db_dxk * d2b_dxidxj);
338 default:
339 throw -1;
340 }
341}
342double AbstractCubic::PI_12(double delta, const std::vector<double>& x, std::size_t idelta) {
343 double bm = bm_term(x);
344 double cm = cm_term();
345 switch (idelta) {
346 case 0:
347 return (1 + (Delta_1 * bm + cm) * rho_r * delta) * (1 + (Delta_2 * bm + cm) * rho_r * delta);
348 case 1:
349 return rho_r * (2 * (bm * Delta_1 + cm) * (bm * Delta_2 + cm) * delta * rho_r + (Delta_1 + Delta_2) * bm + 2 * cm);
350 case 2:
351 return 2 * (Delta_1 * bm + cm) * (Delta_2 * bm + cm) * pow(rho_r, 2);
352 case 3:
353 return 0;
354 case 4:
355 return 0;
356 default:
357 throw -1;
358 }
359}
360double AbstractCubic::d_PI_12_dxi(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, bool xN_independent) {
361 double bm = bm_term(x);
362 double cm = cm_term();
363 double db_dxi = d_bm_term_dxi(x, i, xN_independent);
364 switch (idelta) {
365 case 0:
366 return delta * rho_r * db_dxi * (2 * Delta_1 * Delta_2 * bm * delta * rho_r + (Delta_1 + Delta_2) * (1 + cm * delta * rho_r));
367 case 1:
368 return rho_r * db_dxi * (4 * Delta_1 * Delta_2 * bm * delta * rho_r + (Delta_1 + Delta_2) * (1 + 2 * cm * delta * rho_r));
369 case 2:
370 return 2 * pow(rho_r, 2) * (2 * Delta_1 * Delta_2 * bm + (Delta_1 + Delta_2) * cm) * db_dxi;
371 case 3:
372 return 0;
373 case 4:
374 return 0;
375 default:
376 throw -1;
377 }
378}
379double AbstractCubic::d2_PI_12_dxidxj(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j,
380 bool xN_independent) {
381 double bm = bm_term(x);
382 double cm = cm_term();
383 double db_dxi = d_bm_term_dxi(x, i, xN_independent), db_dxj = d_bm_term_dxi(x, j, xN_independent),
384 d2b_dxidxj = d2_bm_term_dxidxj(x, i, j, xN_independent);
385 switch (idelta) {
386 case 0:
387 return delta * rho_r
388 * (2 * Delta_1 * Delta_2 * delta * rho_r * db_dxi * db_dxj
389 + (2 * Delta_1 * Delta_2 * bm * delta * rho_r + (Delta_1 + Delta_2) * (1 + cm * delta * rho_r)) * d2b_dxidxj);
390 case 1:
391 return rho_r
392 * (4 * Delta_1 * Delta_2 * delta * rho_r * db_dxi * db_dxj
393 + (4 * Delta_1 * Delta_2 * bm * delta * rho_r + (Delta_1 + Delta_2) * (1 + 2 * cm * delta * rho_r)) * d2b_dxidxj);
394 case 2:
395 return 2 * pow(rho_r, 2)
396 * (2 * Delta_1 * Delta_2 * db_dxi * db_dxj + (2 * Delta_1 * Delta_2 * bm + (Delta_1 + Delta_2) * cm) * d2b_dxidxj);
397 case 3:
398 return 0;
399 case 4:
400 return 0;
401 default:
402 throw -1;
403 }
404}
405double AbstractCubic::d3_PI_12_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j, std::size_t k,
406 bool xN_independent) {
407 double bm = bm_term(x);
408 double cm = cm_term();
409 double db_dxi = d_bm_term_dxi(x, i, xN_independent), db_dxj = d_bm_term_dxi(x, j, xN_independent), db_dxk = d_bm_term_dxi(x, k, xN_independent),
410 d2b_dxidxj = d2_bm_term_dxidxj(x, i, j, xN_independent), d2b_dxidxk = d2_bm_term_dxidxj(x, i, k, xN_independent),
411 d2b_dxjdxk = d2_bm_term_dxidxj(x, j, k, xN_independent), d3b_dxidxjdxk = d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent);
412 switch (idelta) {
413 case 0:
414 return delta * rho_r
415 * ((2 * Delta_1 * Delta_2 * bm * delta * rho_r + (Delta_1 + Delta_2) * (1 + cm * delta * rho_r)) * d3b_dxidxjdxk
416 + 2 * Delta_1 * Delta_2 * delta * rho_r * (db_dxi * d2b_dxjdxk + db_dxj * d2b_dxidxk + db_dxk * d2b_dxidxj));
417 case 1:
418 return rho_r
419 * ((4. * Delta_1 * Delta_2 * bm * delta * rho_r + (Delta_1 + Delta_2) * (1 + 2 * cm * delta * rho_r)) * d3b_dxidxjdxk
420 + 4 * Delta_1 * Delta_2 * delta * rho_r * (db_dxi * d2b_dxjdxk + db_dxj * d2b_dxidxk + db_dxk * d2b_dxidxj));
421 default:
422 throw -1;
423 }
424}
425double AbstractCubic::psi_plus(double delta, const std::vector<double>& x, std::size_t idelta) {
426 switch (idelta) {
427 case 0:
428 return A_term(delta, x) * c_term(x) / (Delta_1 - Delta_2);
429 case 1:
430 return rho_r / PI_12(delta, x, 0);
431 case 2:
432 return -rho_r / pow(PI_12(delta, x, 0), 2) * PI_12(delta, x, 1);
433 case 3:
434 return rho_r * (-PI_12(delta, x, 0) * PI_12(delta, x, 2) + 2 * pow(PI_12(delta, x, 1), 2)) / pow(PI_12(delta, x, 0), 3);
435 case 4:
436 // Term -PI_12(delta,x,0)*PI_12(delta,x,3) in the numerator is zero (and removed) since PI_12(delta,x,3) = 0
437 return rho_r * (6 * PI_12(delta, x, 0) * PI_12(delta, x, 1) * PI_12(delta, x, 2) - 6 * pow(PI_12(delta, x, 1), 3))
438 / pow(PI_12(delta, x, 0), 4);
439 default:
440 throw -1;
441 }
442}
443double AbstractCubic::d_psi_plus_dxi(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, bool xN_independent) {
444 double bracket = 0;
445 if (idelta == 0) {
446 return (A_term(delta, x) * d_c_term_dxi(x, i, xN_independent) + c_term(x) * d_A_term_dxi(delta, x, i, xN_independent)) / (Delta_1 - Delta_2);
447 }
448 // All the terms with at least one delta derivative are multiplied by a common term of -rhor/PI12^2
449 // So we just evaluate the bracketed term and then multiply by the common factor in the front
450 switch (idelta) {
451 case 1:
452 bracket = d_PI_12_dxi(delta, x, 0, i, xN_independent);
453 break;
454 case 2:
455 bracket = (d_PI_12_dxi(delta, x, 1, i, xN_independent)
456 + 2 / rho_r * PI_12(delta, x, 0) * PI_12(delta, x, 1) * d_psi_plus_dxi(delta, x, 1, i, xN_independent));
457 break;
458 case 3: {
459 bracket =
460 (d_PI_12_dxi(delta, x, 2, i, xN_independent)
461 + 2 / rho_r * (pow(PI_12(delta, x, 1), 2) + PI_12(delta, x, 0) * PI_12(delta, x, 2)) * d_psi_plus_dxi(delta, x, 1, i, xN_independent)
462 + 4 / rho_r * PI_12(delta, x, 0) * PI_12(delta, x, 1) * d_psi_plus_dxi(delta, x, 2, i, xN_independent));
463 break;
464 }
465 case 4:
466 // d_PI_12_dxi(delta, x, 3, i, xN_independent) = 0, and PI_12(delta,x,0)*PI_12(delta,x,3) = 0, so removed from sum
467 bracket =
468 (6 / rho_r * PI_12(delta, x, 1) * PI_12(delta, x, 2) * d_psi_plus_dxi(delta, x, 1, i, xN_independent)
469 + 6 / rho_r * (pow(PI_12(delta, x, 1), 2) + PI_12(delta, x, 0) * PI_12(delta, x, 2)) * d_psi_plus_dxi(delta, x, 2, i, xN_independent)
470 + 6 / rho_r * PI_12(delta, x, 0) * PI_12(delta, x, 1) * d_psi_plus_dxi(delta, x, 3, i, xN_independent));
471 break;
472 default:
473 throw -1;
474 }
475 return -rho_r / pow(PI_12(delta, x, 0), 2) * bracket;
476}
477double AbstractCubic::d2_psi_plus_dxidxj(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j,
478 bool xN_independent) {
479 double bracket = 0;
480 double PI12 = PI_12(delta, x, 0);
481 if (idelta == 0) {
482 return (A_term(delta, x) * d2_c_term_dxidxj(x, i, j, xN_independent) + c_term(x) * d2_A_term_dxidxj(delta, x, i, j, xN_independent)
483 + d_A_term_dxi(delta, x, i, xN_independent) * d_c_term_dxi(x, j, xN_independent)
484 + d_A_term_dxi(delta, x, j, xN_independent) * d_c_term_dxi(x, i, xN_independent))
485 / (Delta_1 - Delta_2);
486 }
487 // All the terms with at least one delta derivative have a common factor of -1/PI_12^2 out front
488 // so we just calculate the bracketed term and then multiply later on
489 switch (idelta) {
490 case 1:
491 bracket = (rho_r * d2_PI_12_dxidxj(delta, x, 0, i, j, xN_independent)
492 + 2 * PI12 * d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_psi_plus_dxi(delta, x, 1, i, xN_independent));
493 break;
494 case 2:
495 bracket = (rho_r * d2_PI_12_dxidxj(delta, x, 1, i, j, xN_independent)
496 + 2 * (PI12 * d_PI_12_dxi(delta, x, 1, j, xN_independent) + PI_12(delta, x, 1) * d_PI_12_dxi(delta, x, 0, j, xN_independent))
497 * d_psi_plus_dxi(delta, x, 1, i, xN_independent)
498 + 2 * PI12 * PI_12(delta, x, 1) * d2_psi_plus_dxidxj(delta, x, 1, i, j, xN_independent)
499 + 2 * PI12 * d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_psi_plus_dxi(delta, x, 2, i, xN_independent));
500 break;
501 case 3: {
502 bracket =
503 (rho_r * d2_PI_12_dxidxj(delta, x, 2, i, j, xN_independent)
504 + 2 * (PI12 * PI_12(delta, x, 2) + pow(PI_12(delta, x, 1), 2)) * d2_psi_plus_dxidxj(delta, x, 1, i, j, xN_independent)
505 + 4 * (PI12 * d_PI_12_dxi(delta, x, 1, j, xN_independent) + PI_12(delta, x, 1) * d_PI_12_dxi(delta, x, 0, j, xN_independent))
506 * d_psi_plus_dxi(delta, x, 2, i, xN_independent)
507 + 2
508 * (PI12 * d_PI_12_dxi(delta, x, 2, j, xN_independent) + 2 * PI_12(delta, x, 1) * d_PI_12_dxi(delta, x, 1, j, xN_independent)
509 + d_PI_12_dxi(delta, x, 0, j, xN_independent) * PI_12(delta, x, 2))
510 * d_psi_plus_dxi(delta, x, 1, i, xN_independent)
511 + 4 * PI12 * PI_12(delta, x, 1) * d2_psi_plus_dxidxj(delta, x, 2, i, j, xN_independent)
512 + 2 * PI12 * d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_psi_plus_dxi(delta, x, 3, i, xN_independent));
513 break;
514 }
515 case 4:
516 // rho_r*d2_PI_12_dxidxj(delta, x, 3, i, j, xN_independent) = 0
517 // PI_12(delta, x, 3) = 0
518 // PI12*d_PI_12_dxi(delta, x, 3, j, xN_independent) = 0
519 // d_PI_12_dxi(delta, x, 0, j, xN_independent)*PI_12(delta, x, 3) = 0
520 bracket =
521 (+6 * (PI12 * PI_12(delta, x, 2) + pow(PI_12(delta, x, 1), 2)) * d2_psi_plus_dxidxj(delta, x, 2, i, j, xN_independent)
522 + 6 * PI_12(delta, x, 1) * PI_12(delta, x, 2) * d2_psi_plus_dxidxj(delta, x, 1, i, j, xN_independent)
523 + 6 * (PI12 * d_PI_12_dxi(delta, x, 1, j, xN_independent) + PI_12(delta, x, 1) * d_PI_12_dxi(delta, x, 0, j, xN_independent))
524 * d_psi_plus_dxi(delta, x, 3, i, xN_independent)
525 + 6
526 * (PI12 * d_PI_12_dxi(delta, x, 2, j, xN_independent) + 2 * PI_12(delta, x, 1) * d_PI_12_dxi(delta, x, 1, j, xN_independent)
527 + d_PI_12_dxi(delta, x, 0, j, xN_independent) * PI_12(delta, x, 2))
528 * d_psi_plus_dxi(delta, x, 2, i, xN_independent)
529 + 6
530 * (PI_12(delta, x, 1) * d_PI_12_dxi(delta, x, 2, j, xN_independent)
531 + PI_12(delta, x, 2) * d_PI_12_dxi(delta, x, 1, j, xN_independent))
532 * d_psi_plus_dxi(delta, x, 1, i, xN_independent)
533 + 6 * PI12 * PI_12(delta, x, 1) * d2_psi_plus_dxidxj(delta, x, 3, i, j, xN_independent)
534 + 2 * PI12 * d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_psi_plus_dxi(delta, x, 4, i, xN_independent));
535 break;
536 default:
537 throw -1;
538 }
539 return -1 / pow(PI12, 2) * bracket;
540}
541double AbstractCubic::d3_psi_plus_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j,
542 std::size_t k, bool xN_independent) {
543 double PI12 = PI_12(delta, x, 0);
544 switch (idelta) {
545 case 0:
546 return (A_term(delta, x) * d3_c_term_dxidxjdxk(x, i, j, k, xN_independent)
547 + c_term(x) * d3_A_term_dxidxjdxk(delta, x, i, j, k, xN_independent)
548 + d_A_term_dxi(delta, x, i, xN_independent) * d2_c_term_dxidxj(x, j, k, xN_independent)
549 + d_A_term_dxi(delta, x, j, xN_independent) * d2_c_term_dxidxj(x, i, k, xN_independent)
550 + d_A_term_dxi(delta, x, k, xN_independent) * d2_c_term_dxidxj(x, i, j, xN_independent)
551 + d_c_term_dxi(x, i, xN_independent) * d2_A_term_dxidxj(delta, x, j, k, xN_independent)
552 + d_c_term_dxi(x, j, xN_independent) * d2_A_term_dxidxj(delta, x, i, k, xN_independent)
553 + d_c_term_dxi(x, k, xN_independent) * d2_A_term_dxidxj(delta, x, i, j, xN_independent))
554 / (Delta_1 - Delta_2);
555 case 1:
556 return -1 / pow(PI12, 2)
557 * (rho_r * d3_PI_12_dxidxjdxk(delta, x, 0, i, j, k, xN_independent)
558 + 2
559 * (PI12 * d2_PI_12_dxidxj(delta, x, 0, j, k, xN_independent)
560 + d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_PI_12_dxi(delta, x, 0, k, xN_independent))
561 * d_psi_plus_dxi(delta, x, 1, i, xN_independent)
562 + 2 * PI12 * d_PI_12_dxi(delta, x, 0, j, xN_independent) * d2_psi_plus_dxidxj(delta, x, 1, i, k, xN_independent)
563 + 2 * PI12 * d_PI_12_dxi(delta, x, 0, k, xN_independent) * d2_psi_plus_dxidxj(delta, x, 1, i, j, xN_independent));
564 default:
565 throw -1;
566 }
567}
568
569double AbstractCubic::tau_times_a(double tau, const std::vector<double>& x, std::size_t itau) {
570 if (itau == 0) {
571 return tau * am_term(tau, x, 0);
572 } else {
573 return tau * am_term(tau, x, itau) + itau * am_term(tau, x, itau - 1);
574 }
575}
576double AbstractCubic::d_tau_times_a_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent) {
577 if (itau == 0) {
578 return tau * d_am_term_dxi(tau, x, 0, i, xN_independent);
579 } else {
580 return tau * d_am_term_dxi(tau, x, itau, i, xN_independent) + itau * d_am_term_dxi(tau, x, itau - 1, i, xN_independent);
581 }
582}
583double AbstractCubic::d2_tau_times_a_dxidxj(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j,
584 bool xN_independent) {
585 if (itau == 0) {
586 return tau * d2_am_term_dxidxj(tau, x, 0, i, j, xN_independent);
587 } else {
588 return tau * d2_am_term_dxidxj(tau, x, itau, i, j, xN_independent) + itau * d2_am_term_dxidxj(tau, x, itau - 1, i, j, xN_independent);
589 }
590}
591double AbstractCubic::d3_tau_times_a_dxidxjdxk(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j,
592 std::size_t k, bool xN_independent) {
593 if (itau == 0) {
594 return tau * d3_am_term_dxidxjdxk(tau, x, 0, i, j, k, xN_independent);
595 } else {
596 return tau * d3_am_term_dxidxjdxk(tau, x, itau, i, j, k, xN_independent)
597 + itau * d3_am_term_dxidxjdxk(tau, x, itau - 1, i, j, k, xN_independent);
598 }
599}
600double AbstractCubic::alphar(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta) {
601 return psi_minus(delta, x, itau, idelta) - tau_times_a(tau, x, itau) / (R_u * T_r) * psi_plus(delta, x, idelta);
602}
603double AbstractCubic::d_alphar_dxi(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
604 bool xN_independent) {
605 return (d_psi_minus_dxi(delta, x, itau, idelta, i, xN_independent)
606 - 1 / (R_u * T_r)
607 * (d_tau_times_a_dxi(tau, x, itau, i, xN_independent) * psi_plus(delta, x, idelta)
608 + tau_times_a(tau, x, itau) * d_psi_plus_dxi(delta, x, idelta, i, xN_independent)));
609}
610double AbstractCubic::d2_alphar_dxidxj(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
611 std::size_t j, bool xN_independent) {
612 return (d2_psi_minus_dxidxj(delta, x, itau, idelta, i, j, xN_independent)
613 - 1 / (R_u * T_r)
614 * (d2_tau_times_a_dxidxj(tau, x, itau, i, j, xN_independent) * psi_plus(delta, x, idelta)
615 + d_tau_times_a_dxi(tau, x, itau, i, xN_independent) * d_psi_plus_dxi(delta, x, idelta, j, xN_independent)
616 + d_tau_times_a_dxi(tau, x, itau, j, xN_independent) * d_psi_plus_dxi(delta, x, idelta, i, xN_independent)
617 + tau_times_a(tau, x, itau) * d2_psi_plus_dxidxj(delta, x, idelta, i, j, xN_independent)));
618}
619double AbstractCubic::d3_alphar_dxidxjdxk(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
620 std::size_t j, std::size_t k, bool xN_independent) {
621 return (d3_psi_minus_dxidxjdxk(delta, x, itau, idelta, i, j, k, xN_independent)
622 - 1 / (R_u * T_r)
623 * (d2_tau_times_a_dxidxj(tau, x, itau, i, j, xN_independent) * d_psi_plus_dxi(delta, x, idelta, k, xN_independent)
624 + d3_tau_times_a_dxidxjdxk(tau, x, itau, i, j, k, xN_independent) * psi_plus(delta, x, idelta)
625
626 + d_tau_times_a_dxi(tau, x, itau, i, xN_independent) * d2_psi_plus_dxidxj(delta, x, idelta, j, k, xN_independent)
627 + d2_tau_times_a_dxidxj(tau, x, itau, i, k, xN_independent) * d_psi_plus_dxi(delta, x, idelta, j, xN_independent)
628
629 + d_tau_times_a_dxi(tau, x, itau, j, xN_independent) * d2_psi_plus_dxidxj(delta, x, idelta, i, k, xN_independent)
630 + d2_tau_times_a_dxidxj(tau, x, itau, j, k, xN_independent) * d_psi_plus_dxi(delta, x, idelta, i, xN_independent)
631
632 + tau_times_a(tau, x, itau) * d3_psi_plus_dxidxjdxk(delta, x, idelta, i, j, k, xN_independent)
633 + d_tau_times_a_dxi(tau, x, itau, k, xN_independent) * d2_psi_plus_dxidxj(delta, x, idelta, i, j, xN_independent)));
634}
635
636double SRK::a0_ii(std::size_t i) {
637 // Values from Soave, 1972 (Equilibrium constants from a ..)
638 double a = 0.42747 * R_u * R_u * Tc[i] * Tc[i] / pc[i];
639 return a;
640}
641double SRK::b0_ii(std::size_t i) {
642 // Values from Soave, 1972 (Equilibrium constants from a ..)
643 double b = 0.08664 * R_u * Tc[i] / pc[i];
644 return b;
645}
646double SRK::m_ii(std::size_t i) {
647 // Values from Soave, 1972 (Equilibrium constants from a ..)
648 double omega = acentric[i];
649 double m = 0.480 + 1.574 * omega - 0.176 * omega * omega;
650 return m;
651}
652
653double PengRobinson::a0_ii(std::size_t i) {
654 double a = 0.45724 * R_u * R_u * Tc[i] * Tc[i] / pc[i];
655 return a;
656}
657double PengRobinson::b0_ii(std::size_t i) {
658 double b = 0.07780 * R_u * Tc[i] / pc[i];
659 return b;
660}
661double PengRobinson::m_ii(std::size_t i) {
662 double omega = acentric[i];
663 double m = 0.37464 + 1.54226 * omega - 0.26992 * omega * omega;
664 return m;
665}
Generated on Thu Jun 5 2025 04:10:24 for CoolProp by doxygen 1.9.4