CoolProp  6.6.0
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 
5 double 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 
31 double 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 
77 double 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 
120 AbstractCubic::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 
131 void 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 
147 double 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 }
156 double 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 }
171 double 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 
180 double 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 
185 double 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 }
192 double 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 }
199 double 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 }
202 double 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 
206 double AbstractCubic::cm_term() {
207  return cm;
208 }
209 
210 double AbstractCubic::aii_term(double tau, std::size_t i, std::size_t itau) {
211  return alpha[i]->term(tau, itau);
212 }
213 double 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 }
232 double 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 }
254 double 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 }
274 double 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 }
296 double 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 }
322 double 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 }
342 double 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 }
360 double 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 }
379 double 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 }
405 double 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 }
425 double 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 }
443 double 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 }
477 double 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 }
541 double 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 
569 double 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 }
576 double 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 }
583 double 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 }
591 double 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 }
600 double 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 }
603 double 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 }
610 double 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 }
619 double 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 
636 double 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 }
641 double 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 }
646 double 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 
653 double 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 }
657 double PengRobinson::b0_ii(std::size_t i) {
658  double b = 0.07780 * R_u * Tc[i] / pc[i];
659  return b;
660 }
661 double 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 Wed Nov 29 2023 14:26:23 for CoolProp by doxygen 1.9.1