Loading [MathJax]/extensions/TeX/AMSsymbols.js
CoolProp 6.8.1dev
An open-source fluid property and humid air property database
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Helmholtz.h
Go to the documentation of this file.
1
2#ifndef HELMHOLTZ_H
3#define HELMHOLTZ_H
4
5#include <vector>
6#include "rapidjson_include.h"
7//#include "Eigen/Core"
8#include "time.h"
9#include "CachedElement.h"
10#include "Backends/Cubics/GeneralizedCubic.h"
11#include "crossplatform_shared_ptr.h"
12
13#if ENABLE_CATCH
14# include "MultiComplex/MultiComplex.hpp"
15#endif
16
17namespace CoolProp {
18
19// #############################################################################
20// #############################################################################
21// #############################################################################
22// RESIDUAL TERMS
23// #############################################################################
24// #############################################################################
25// #############################################################################
26
27#define LIST_OF_DERIVATIVE_VARIABLES \
28 X(alphar) \
29 X(dalphar_ddelta) \
30 X(dalphar_dtau) \
31 X(d2alphar_ddelta2) \
32 X(d2alphar_dtau2) \
33 X(d2alphar_ddelta_dtau) \
34 X(d3alphar_ddelta3) \
35 X(d3alphar_ddelta_dtau2) \
36 X(d3alphar_ddelta2_dtau) \
37 X(d3alphar_dtau3) \
38 X(d4alphar_ddelta4) \
39 X(d4alphar_ddelta3_dtau) \
40 X(d4alphar_ddelta2_dtau2) \
41 X(d4alphar_ddelta_dtau3) \
42 X(d4alphar_dtau4) \
43 X(delta_x_dalphar_ddelta) \
44 X(tau_x_dalphar_dtau) \
45 X(delta2_x_d2alphar_ddelta2) \
46 X(deltatau_x_d2alphar_ddelta_dtau) \
47 X(tau2_x_d2alphar_dtau2)
48
49struct HelmholtzDerivatives
50{
51#define X(name) CoolPropDbl name;
52 LIST_OF_DERIVATIVE_VARIABLES
53#undef X
54 CoolPropDbl tau, delta, T_red, rhomolar_red;
55
56 void reset(CoolPropDbl v) {
57#define X(name) name = v;
58 LIST_OF_DERIVATIVE_VARIABLES
59#undef X
60 }
61 HelmholtzDerivatives operator+(const HelmholtzDerivatives& other) const {
62 HelmholtzDerivatives _new;
63#define X(name) _new.name = name + other.name;
64 LIST_OF_DERIVATIVE_VARIABLES
65#undef X
66 return _new;
67 }
68 HelmholtzDerivatives operator*(const CoolPropDbl& other) const {
69 HelmholtzDerivatives _new;
70#define X(name) _new.name = name * other;
71 LIST_OF_DERIVATIVE_VARIABLES
72#undef X
73 return _new;
74 }
75 HelmholtzDerivatives() {
76 reset(0.0);
77 T_red = _HUGE;
78 rhomolar_red = _HUGE;
79 };
81 double get(std::size_t itau, std::size_t idelta) {
82 if (itau == 0) {
83 if (idelta == 0) {
84 return alphar;
85 } else if (idelta == 1) {
86 return dalphar_ddelta;
87 } else if (idelta == 2) {
88 return d2alphar_ddelta2;
89 } else if (idelta == 3) {
90 return d3alphar_ddelta3;
91 } else if (idelta == 4) {
92 return d4alphar_ddelta4;
93 } else {
94 throw ValueError();
95 }
96 } else if (itau == 1) {
97 if (idelta == 0) {
98 return dalphar_dtau;
99 } else if (idelta == 1) {
100 return d2alphar_ddelta_dtau;
101 } else if (idelta == 2) {
102 return d3alphar_ddelta2_dtau;
103 } else if (idelta == 3) {
104 return d4alphar_ddelta3_dtau;
105 } else {
106 throw ValueError();
107 }
108 } else if (itau == 2) {
109 if (idelta == 0) {
110 return d2alphar_dtau2;
111 } else if (idelta == 1) {
112 return d3alphar_ddelta_dtau2;
113 } else if (idelta == 2) {
114 return d4alphar_ddelta2_dtau2;
115 } else {
116 throw ValueError();
117 }
118 } else if (itau == 3) {
119 if (idelta == 0) {
120 return d3alphar_dtau3;
121 } else if (idelta == 1) {
122 return d4alphar_ddelta_dtau3;
123 } else {
124 throw ValueError();
125 }
126 } else if (itau == 4) {
127 if (idelta == 0) {
128 return d4alphar_dtau4;
129 } else {
130 throw ValueError();
131 }
132 } else {
133 throw ValueError();
134 }
135 }
136};
137#undef LIST_OF_DERIVATIVE_VARIABLES
138
140
165class BaseHelmholtzTerm
166{
167 public:
168 BaseHelmholtzTerm(){};
169 virtual ~BaseHelmholtzTerm(){};
170
172
175 virtual CoolPropDbl base(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
176 HelmholtzDerivatives deriv;
177 all(tau, delta, deriv);
178 return deriv.alphar;
179 };
181
184 virtual CoolPropDbl dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
185 HelmholtzDerivatives deriv;
186 all(tau, delta, deriv);
187 return deriv.dalphar_dtau;
188 };
190
193 virtual CoolPropDbl dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
194 HelmholtzDerivatives deriv;
195 all(tau, delta, deriv);
196 return deriv.d2alphar_dtau2;
197 };
199
202 virtual CoolPropDbl dDelta_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
203 HelmholtzDerivatives deriv;
204 all(tau, delta, deriv);
205 return deriv.d2alphar_ddelta_dtau;
206 };
208
211 virtual CoolPropDbl dDelta(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
212 HelmholtzDerivatives deriv;
213 all(tau, delta, deriv);
214 return deriv.dalphar_ddelta;
215 };
217
220 virtual CoolPropDbl dDelta2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
221 HelmholtzDerivatives deriv;
222 all(tau, delta, deriv);
223 return deriv.d2alphar_ddelta2;
224 };
226
229 virtual CoolPropDbl dDelta2_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
230 HelmholtzDerivatives deriv;
231 all(tau, delta, deriv);
232 return deriv.d3alphar_ddelta2_dtau;
233 };
235
238 virtual CoolPropDbl dDelta_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
239 HelmholtzDerivatives deriv;
240 all(tau, delta, deriv);
241 return deriv.d3alphar_ddelta_dtau2;
242 };
244
247 virtual CoolPropDbl dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
248 HelmholtzDerivatives deriv;
249 all(tau, delta, deriv);
250 return deriv.d3alphar_dtau3;
251 };
253
256 virtual CoolPropDbl dDelta3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
257 HelmholtzDerivatives deriv;
258 all(tau, delta, deriv);
259 return deriv.d3alphar_ddelta3;
260 };
262
265 virtual CoolPropDbl dTau4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
266 HelmholtzDerivatives deriv;
267 all(tau, delta, deriv);
268 return deriv.d4alphar_dtau4;
269 };
270 virtual CoolPropDbl dDelta_dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
271 HelmholtzDerivatives deriv;
272 all(tau, delta, deriv);
273 return deriv.d4alphar_ddelta_dtau3;
274 };
275 virtual CoolPropDbl dDelta2_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
276 HelmholtzDerivatives deriv;
277 all(tau, delta, deriv);
278 return deriv.d4alphar_ddelta2_dtau2;
279 };
280 virtual CoolPropDbl dDelta3_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
281 HelmholtzDerivatives deriv;
282 all(tau, delta, deriv);
283 return deriv.d4alphar_ddelta3_dtau;
284 };
285 virtual CoolPropDbl dDelta4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
286 HelmholtzDerivatives deriv;
287 all(tau, delta, deriv);
288 return deriv.d4alphar_ddelta4;
289 };
290
291 virtual void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) = 0;
292 #if ENABLE_CATCH
293 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const {
294 throw CoolProp::NotImplementedError("The mcx derivative function was not implemented");
295 }
296 #endif
297};
298
299struct ResidualHelmholtzGeneralizedExponentialElement
300{
302 CoolPropDbl n, d, t;
305 CoolPropDbl c, l_double, omega, m_double, eta1, epsilon1, eta2, epsilon2, beta1, gamma1, beta2, gamma2;
307 int l_int, m_int;
309 bool l_is_int, m_is_int;
310
311 ResidualHelmholtzGeneralizedExponentialElement() {
312 n = 0;
313 d = 0;
314 t = 0;
315 c = 0;
316 l_double = 0;
317 omega = 0;
318 m_double = 0;
319 eta1 = 0;
320 epsilon1 = 0;
321 eta2 = 0;
322 epsilon2 = 0;
323 beta1 = 0;
324 gamma1 = 0;
325 beta2 = 0;
326 gamma2 = 0;
327 l_int = 0;
328 m_int = 0;
329 l_is_int = false;
330 m_is_int = true;
331 };
332};
341class ResidualHelmholtzGeneralizedExponential : public BaseHelmholtzTerm
342{
343
344 public:
345 bool delta_li_in_u, tau_mi_in_u, eta1_in_u, eta2_in_u, beta1_in_u, beta2_in_u, finished;
346 std::vector<CoolPropDbl> s;
347 std::size_t N;
348
349 // These variables are for the exp(u) part
350 // u is given by -c*delta^l_i-omega*tau^m_i-eta1*(delta-epsilon1)-eta2*(delta-epsilon2)^2-beta1*(tau-gamma1)-beta2*(tau-gamma2)^2
351 std::vector<double> n, d, t, c, l_double, omega, m_double, eta1, epsilon1, eta2, epsilon2, beta1, gamma1, beta2, gamma2;
352 // If l_i or m_i are integers, we will store them as integers in order to call pow(double, int) rather than pow(double, double)
353 std::vector<int> l_int, m_int;
354
355 //Eigen::ArrayXd uE, du_ddeltaE, du_dtauE, d2u_ddelta2E, d2u_dtau2E, d3u_ddelta3E, d3u_dtau3E;
356
357 std::vector<ResidualHelmholtzGeneralizedExponentialElement> elements;
358 // Default Constructor
359 ResidualHelmholtzGeneralizedExponential()
360 : delta_li_in_u(false), tau_mi_in_u(false), eta1_in_u(false), eta2_in_u(false), beta1_in_u(false), beta2_in_u(false), finished(false), N(0){};
366 void add_Power(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
367 const std::vector<CoolPropDbl>& l) {
368 for (std::size_t i = 0; i < n.size(); ++i) {
369 ResidualHelmholtzGeneralizedExponentialElement el;
370 el.n = n[i];
371 el.d = d[i];
372 el.t = t[i];
373 el.l_double = l[i];
374 el.l_int = (int)el.l_double;
375 if (el.l_double > 0)
376 el.c = 1.0;
377 else
378 el.c = 0.0;
379 elements.push_back(el);
380 }
381 delta_li_in_u = true;
382 };
388 void add_Exponential(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
389 const std::vector<CoolPropDbl>& g, const std::vector<CoolPropDbl>& l) {
390 for (std::size_t i = 0; i < n.size(); ++i) {
391 ResidualHelmholtzGeneralizedExponentialElement el;
392 el.n = n[i];
393 el.d = d[i];
394 el.t = t[i];
395 el.c = g[i];
396 el.l_double = l[i];
397 el.l_int = (int)el.l_double;
398 elements.push_back(el);
399 }
400 delta_li_in_u = true;
401 }
407 void add_Gaussian(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
408 const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& beta,
409 const std::vector<CoolPropDbl>& gamma) {
410 for (std::size_t i = 0; i < n.size(); ++i) {
411 ResidualHelmholtzGeneralizedExponentialElement el;
412 el.n = n[i];
413 el.d = d[i];
414 el.t = t[i];
415 el.eta2 = eta[i];
416 el.epsilon2 = epsilon[i];
417 el.beta2 = beta[i];
418 el.gamma2 = gamma[i];
419 elements.push_back(el);
420 }
421 eta2_in_u = true;
422 beta2_in_u = true;
423 };
429 void add_GERG2008Gaussian(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
430 const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& beta,
431 const std::vector<CoolPropDbl>& gamma) {
432 for (std::size_t i = 0; i < n.size(); ++i) {
433 ResidualHelmholtzGeneralizedExponentialElement el;
434 el.n = n[i];
435 el.d = d[i];
436 el.t = t[i];
437 el.eta2 = eta[i];
438 el.epsilon2 = epsilon[i];
439 el.eta1 = beta[i];
440 el.epsilon1 = gamma[i];
441 elements.push_back(el);
442 }
443 eta2_in_u = true;
444 eta1_in_u = true;
445 };
451 void add_Lemmon2005(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
452 const std::vector<CoolPropDbl>& l, const std::vector<CoolPropDbl>& m) {
453 for (std::size_t i = 0; i < n.size(); ++i) {
454 ResidualHelmholtzGeneralizedExponentialElement el;
455 el.n = n[i];
456 el.d = d[i];
457 el.t = t[i];
458 el.c = 1.0;
459 el.omega = 1.0;
460 el.l_double = l[i];
461 el.m_double = m[i];
462 el.l_int = (int)el.l_double;
463 el.m_int = (int)el.m_double;
464 elements.push_back(el);
465 }
466 delta_li_in_u = true;
467 tau_mi_in_u = true;
468 };
469
470 void finish() {
471 n.resize(elements.size());
472 d.resize(elements.size());
473 t.resize(elements.size());
474 c.resize(elements.size());
475 omega.resize(elements.size());
476 l_double.resize(elements.size());
477 l_int.resize(elements.size());
478 m_double.resize(elements.size());
479 m_int.resize(elements.size());
480 epsilon2.resize(elements.size());
481 eta2.resize(elements.size());
482 gamma2.resize(elements.size());
483 beta2.resize(elements.size());
484
485 for (std::size_t i = 0; i < elements.size(); ++i) {
486 n[i] = elements[i].n;
487 d[i] = elements[i].d;
488 t[i] = elements[i].t;
489 c[i] = elements[i].c;
490 omega[i] = elements[i].omega;
491 l_double[i] = elements[i].l_double;
492 l_int[i] = elements[i].l_int;
493 m_double[i] = elements[i].m_double;
494 m_int[i] = elements[i].m_int;
495 epsilon2[i] = elements[i].epsilon2;
496 eta2[i] = elements[i].eta2;
497 gamma2[i] = elements[i].gamma2;
498 beta2[i] = elements[i].beta2;
499
500 // See if l is an integer, and store a flag if it is
501 elements[i].l_is_int = (std::abs(static_cast<long>(elements[i].l_double) - elements[i].l_double) < 1e-14);
502 }
503 // uE.resize(elements.size());
504 // du_ddeltaE.resize(elements.size());
505 // du_dtauE.resize(elements.size());
506 // d2u_ddelta2E.resize(elements.size());
507 // d2u_dtau2E.resize(elements.size());
508 // d3u_ddelta3E.resize(elements.size());
509 // d3u_dtau3E.resize(elements.size());
510
511 finished = true;
512 };
513
514 void to_json(rapidjson::Value& el, rapidjson::Document& doc);
515
516 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
517 //void allEigen(const CoolPropDbl &tau, const CoolPropDbl &delta, HelmholtzDerivatives &derivs) throw();
518
519 #if ENABLE_CATCH
520 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
521 #endif
522};
523
524struct ResidualHelmholtzNonAnalyticElement
525{
526 CoolPropDbl n, a, b, beta, A, B, C, D;
527};
528class ResidualHelmholtzNonAnalytic : public BaseHelmholtzTerm
529{
530
531 public:
532 std::size_t N;
533 std::vector<CoolPropDbl> s;
534 std::vector<ResidualHelmholtzNonAnalyticElement> elements;
536 ResidualHelmholtzNonAnalytic() {
537 N = 0;
538 };
540 ~ResidualHelmholtzNonAnalytic(){};
542 ResidualHelmholtzNonAnalytic(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& a, const std::vector<CoolPropDbl>& b,
543 const std::vector<CoolPropDbl>& beta, const std::vector<CoolPropDbl>& A, const std::vector<CoolPropDbl>& B,
544 const std::vector<CoolPropDbl>& C, const std::vector<CoolPropDbl>& D) {
545 N = n.size();
546 s.resize(N);
547 for (std::size_t i = 0; i < n.size(); ++i) {
548 ResidualHelmholtzNonAnalyticElement el;
549 el.n = n[i];
550 el.a = a[i];
551 el.b = b[i];
552 el.beta = beta[i];
553 el.A = A[i];
554 el.B = B[i];
555 el.C = C[i];
556 el.D = D[i];
557 elements.push_back(el);
558 }
559 };
560 void to_json(rapidjson::Value& el, rapidjson::Document& doc);
561 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
562#if ENABLE_CATCH
563 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
564#endif
565};
566
567class ResidualHelmholtzGeneralizedCubic : public BaseHelmholtzTerm
568{
569 protected:
570 shared_ptr<AbstractCubic> m_abstractcubic;
571 std::vector<double> z;
572 public:
573 bool enabled;
574
576 ResidualHelmholtzGeneralizedCubic() {
577 enabled = false;
578 };
580 ResidualHelmholtzGeneralizedCubic(shared_ptr<AbstractCubic>& ac) : m_abstractcubic(ac) {
581 enabled = true;
582 z = std::vector<double>(1, 1); // Init the vector to [1.0]
583 };
584
585 void to_json(rapidjson::Value& el, rapidjson::Document& doc);
586 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
587};
588
589class ResidualHelmholtzGaoB : public BaseHelmholtzTerm
590{
591 protected:
592 std::vector<double> n, t, d, eta, beta, gamma, epsilon, b;
593
594 public:
595 bool enabled;
596
598 ResidualHelmholtzGaoB() {
599 enabled = false;
600 };
601
603 ResidualHelmholtzGaoB(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& t, const std::vector<CoolPropDbl>& d,
604 const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& beta, const std::vector<CoolPropDbl>& gamma,
605 const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& b)
606 : n(n), t(t), d(d), eta(eta), beta(beta), gamma(gamma), epsilon(epsilon), b(b) {
607 enabled = true;
608 };
609
610 void to_json(rapidjson::Value& el, rapidjson::Document& doc);
611 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
612
613 #if ENABLE_CATCH
614 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
615 #endif
616};
617
619class ResidualHelmholtzXiangDeiters : public BaseHelmholtzTerm
620{
621
622 public:
623 bool enabled;
624 ResidualHelmholtzGeneralizedExponential phi0, phi1, phi2;
625 CoolPropDbl Tc, pc, rhomolarc, acentric, R, theta;
627 ResidualHelmholtzXiangDeiters() : Tc(_HUGE), pc(_HUGE), rhomolarc(_HUGE), acentric(_HUGE), R(_HUGE), theta(_HUGE) {
628 enabled = false;
629 };
631 ResidualHelmholtzXiangDeiters(const CoolPropDbl Tc, const CoolPropDbl pc, const CoolPropDbl rhomolarc, const CoolPropDbl acentric,
632 const CoolPropDbl R);
633 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
634#if ENABLE_CATCH
635 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
636#endif
637};
638
639class ResidualHelmholtzSAFTAssociating : public BaseHelmholtzTerm
640{
641
642 protected:
643 double a, m, epsilonbar, vbarn, kappabar;
644
645 CoolPropDbl Deltabar(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
646 CoolPropDbl dDeltabar_ddelta__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
647 CoolPropDbl d2Deltabar_ddelta2__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
648 CoolPropDbl dDeltabar_dtau__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
649 CoolPropDbl d2Deltabar_dtau2__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
650 CoolPropDbl d2Deltabar_ddelta_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
651 CoolPropDbl d3Deltabar_dtau3__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
652 CoolPropDbl d3Deltabar_ddelta_dtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
653 CoolPropDbl d3Deltabar_ddelta3__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
654 CoolPropDbl d3Deltabar_ddelta2_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
655
656 CoolPropDbl X(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
657 CoolPropDbl dX_dDeltabar__constdelta(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
658 CoolPropDbl dX_ddelta__constDeltabar(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
659 CoolPropDbl dX_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
660 CoolPropDbl dX_ddelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
661 CoolPropDbl d2X_dtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
662 CoolPropDbl d2X_ddeltadtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
663 CoolPropDbl d2X_ddelta2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
664
665 CoolPropDbl d3X_dtau3(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
666 CoolPropDbl d3X_ddelta3(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
667 CoolPropDbl d3X_ddeltadtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
668 CoolPropDbl d3X_ddelta2dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
669
670 CoolPropDbl g(const CoolPropDbl& eta) const;
671 CoolPropDbl dg_deta(const CoolPropDbl& eta) const;
672 CoolPropDbl d2g_deta2(const CoolPropDbl& eta) const;
673 CoolPropDbl d3g_deta3(const CoolPropDbl& eta) const;
674 CoolPropDbl eta(const CoolPropDbl& delta) const;
675
676 public:
678 ResidualHelmholtzSAFTAssociating() : a(_HUGE), m(_HUGE), epsilonbar(_HUGE), vbarn(_HUGE), kappabar(_HUGE) {
679 disabled = true;
680 };
681
682 // Constructor
683 ResidualHelmholtzSAFTAssociating(double a, double m, double epsilonbar, double vbarn, double kappabar)
684 : a(a), m(m), epsilonbar(epsilonbar), vbarn(vbarn), kappabar(kappabar) {
685 disabled = false;
686 };
687
688 bool disabled;
689
690 //Destructor. No Implementation
691 ~ResidualHelmholtzSAFTAssociating(){};
692
693 void to_json(rapidjson::Value& el, rapidjson::Document& doc);
694
695 CoolPropDbl dTau4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
696 return 1e99;
697 };
698 CoolPropDbl dDelta_dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
699 return 1e99;
700 };
701 CoolPropDbl dDelta2_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
702 return 1e99;
703 };
704 CoolPropDbl dDelta3_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
705 return 1e99;
706 };
707 CoolPropDbl dDelta4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
708 return 1e99;
709 };
710
711 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& deriv) throw();
712};
713
714class BaseHelmholtzContainer
715{
716 protected:
717 CachedElement _base, _dDelta, _dTau, _dDelta2, _dTau2, _dDelta_dTau, _dDelta3, _dDelta2_dTau, _dDelta_dTau2, _dTau3;
718 CachedElement _dDelta4, _dDelta3_dTau, _dDelta2_dTau2, _dDelta_dTau3, _dTau4;
719
720 public:
721 void clear() {
722 _base.clear();
723 _dDelta.clear();
724 _dTau.clear();
725 _dDelta2.clear();
726 _dTau2.clear();
727 _dDelta_dTau.clear();
728 _dDelta3.clear();
729 _dTau3.clear();
730 _dDelta2_dTau.clear();
731 _dDelta_dTau2.clear();
732 _dDelta4.clear();
733 _dDelta3_dTau.clear();
734 _dDelta2_dTau2.clear();
735 _dDelta_dTau3.clear();
736 _dTau4.clear();
737 };
738
739 virtual void empty_the_EOS() = 0;
740 virtual HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values) = 0;
741
742 CoolPropDbl base(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
743 if (!_base || dont_use_cache)
744 return all(tau, delta, false).alphar;
745 else
746 return _base;
747 };
748 CoolPropDbl dDelta(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
749 if (!_dDelta || dont_use_cache)
750 return all(tau, delta, false).dalphar_ddelta;
751 else
752 return _dDelta;
753 };
754 CoolPropDbl dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
755 if (!_dTau || dont_use_cache)
756 return all(tau, delta, false).dalphar_dtau;
757 else
758 return _dTau;
759 };
760 CoolPropDbl dDelta2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
761 if (!_dDelta2 || dont_use_cache)
762 return all(tau, delta, false).d2alphar_ddelta2;
763 else
764 return _dDelta2;
765 };
766 CoolPropDbl dDelta_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
767 if (!_dDelta_dTau || dont_use_cache)
768 return all(tau, delta, false).d2alphar_ddelta_dtau;
769 else
770 return _dDelta_dTau;
771 };
772 CoolPropDbl dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
773 if (!_dTau2 || dont_use_cache)
774 return all(tau, delta, false).d2alphar_dtau2;
775 else
776 return _dTau2;
777 };
778 CoolPropDbl dDelta3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
779 if (!_dDelta3 || dont_use_cache)
780 return all(tau, delta, false).d3alphar_ddelta3;
781 else
782 return _dDelta3;
783 };
784 CoolPropDbl dDelta2_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
785 if (!_dDelta2_dTau || dont_use_cache)
786 return all(tau, delta, false).d3alphar_ddelta2_dtau;
787 else
788 return _dDelta2_dTau;
789 };
790 CoolPropDbl dDelta_dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
791 if (!_dDelta_dTau2 || dont_use_cache)
792 return all(tau, delta, false).d3alphar_ddelta_dtau2;
793 else
794 return _dDelta_dTau2;
795 };
796 CoolPropDbl dTau3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
797 if (!_dTau3 || dont_use_cache)
798 return all(tau, delta, false).d3alphar_dtau3;
799 else
800 return _dTau3;
801 };
802 CoolPropDbl dDelta4(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
803 return all(tau, delta, false).d4alphar_ddelta4;
804 };
805 CoolPropDbl dDelta3_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
806 return all(tau, delta, false).d4alphar_ddelta3_dtau;
807 };
808 CoolPropDbl dDelta2_dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
809 return all(tau, delta, false).d4alphar_ddelta2_dtau2;
810 };
811 CoolPropDbl dDelta_dTau3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
812 return all(tau, delta, false).d4alphar_ddelta_dtau3;
813 };
814 CoolPropDbl dTau4(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
815 return all(tau, delta, false).d4alphar_dtau4;
816 };
817};
818
819class ResidualHelmholtzContainer : public BaseHelmholtzContainer
820{
821 public:
822 ResidualHelmholtzNonAnalytic NonAnalytic;
823 ResidualHelmholtzSAFTAssociating SAFT;
824 ResidualHelmholtzGeneralizedExponential GenExp;
825 ResidualHelmholtzGeneralizedCubic cubic;
826 ResidualHelmholtzXiangDeiters XiangDeiters;
827 ResidualHelmholtzGaoB GaoB;
828
829 void empty_the_EOS() {
830 NonAnalytic = ResidualHelmholtzNonAnalytic();
831 SAFT = ResidualHelmholtzSAFTAssociating();
832 GenExp = ResidualHelmholtzGeneralizedExponential();
833 cubic = ResidualHelmholtzGeneralizedCubic();
834 XiangDeiters = ResidualHelmholtzXiangDeiters();
835 GaoB = ResidualHelmholtzGaoB();
836 };
837
838 HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values = false) {
839 HelmholtzDerivatives derivs; // zeros out the elements
840 GenExp.all(tau, delta, derivs);
841 NonAnalytic.all(tau, delta, derivs);
842 SAFT.all(tau, delta, derivs);
843 cubic.all(tau, delta, derivs);
844 XiangDeiters.all(tau, delta, derivs);
845 GaoB.all(tau, delta, derivs);
846 if (cache_values) {
847 _base = derivs.alphar;
848 _dDelta = derivs.dalphar_ddelta;
849 _dTau = derivs.dalphar_dtau;
850 _dDelta2 = derivs.d2alphar_ddelta2;
851 _dTau2 = derivs.d2alphar_dtau2;
852 _dDelta_dTau = derivs.d2alphar_ddelta_dtau;
853 _dDelta3 = derivs.d3alphar_ddelta3;
854 _dTau3 = derivs.d3alphar_dtau3;
855 _dDelta2_dTau = derivs.d3alphar_ddelta2_dtau;
856 _dDelta_dTau2 = derivs.d3alphar_ddelta_dtau2;
857 }
858 return derivs;
859 };
860};
861
862// #############################################################################
863// #############################################################################
864// #############################################################################
865// IDEAL GAS TERMS
866// #############################################################################
867// #############################################################################
868// #############################################################################
869
871
876class IdealHelmholtzLead : public BaseHelmholtzTerm
877{
878
879 private:
880 CoolPropDbl a1, a2;
881 bool enabled;
882
883 public:
884 // Default constructor
885 IdealHelmholtzLead() : a1(_HUGE), a2(_HUGE), enabled(false) {}
886
887 // Constructor
888 IdealHelmholtzLead(CoolPropDbl a1, CoolPropDbl a2) : a1(a1), a2(a2), enabled(true) {}
889
890 bool is_enabled() const {
891 return enabled;
892 }
893
894 void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
895 el.AddMember("type", "IdealHelmholtzLead", doc.GetAllocator());
896 el.AddMember("a1", static_cast<double>(a1), doc.GetAllocator());
897 el.AddMember("a2", static_cast<double>(a2), doc.GetAllocator());
898 };
899
900 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
901};
902
904
909class IdealHelmholtzEnthalpyEntropyOffset : public BaseHelmholtzTerm
910{
911 private:
912 CoolPropDbl a1, a2; // Use these variables internally
913 std::string reference;
914 bool enabled;
915
916 public:
917 IdealHelmholtzEnthalpyEntropyOffset() : a1(_HUGE), a2(_HUGE), enabled(false) {}
918
919 // Constructor
920 IdealHelmholtzEnthalpyEntropyOffset(CoolPropDbl a1, CoolPropDbl a2, const std::string& ref) : a1(a1), a2(a2), reference(ref), enabled(true) {}
921
922 // Set the values in the class
923 void set(CoolPropDbl a1, CoolPropDbl a2, const std::string& ref) {
924 // If it doesn't already exist, just set the values
925 if (enabled == false) {
926 this->a1 = a1;
927 this->a2 = a2;
928 enabled = true;
929 } else if (ref == "DEF") {
930 this->a1 = 0.0;
931 this->a2 = 0.0;
932 enabled = false;
933 } else {
934 // Otherwise, increment the values
935 this->a1 += a1;
936 this->a2 += a2;
937 enabled = true;
938 }
939 this->reference = ref;
940 }
941
942 bool is_enabled() const {
943 return enabled;
944 };
945
946 void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
947 el.AddMember("type", "IdealHelmholtzEnthalpyEntropyOffset", doc.GetAllocator());
948 el.AddMember("a1", static_cast<double>(a1), doc.GetAllocator());
949 el.AddMember("a2", static_cast<double>(a2), doc.GetAllocator());
950 };
951 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
952};
953
959class IdealHelmholtzLogTau : public BaseHelmholtzTerm
960{
961 private:
962 CoolPropDbl a1;
963 bool enabled;
964
965 public:
967 IdealHelmholtzLogTau() : a1(_HUGE), enabled(false) {}
968
969 // Constructor
970 IdealHelmholtzLogTau(CoolPropDbl a1) : a1(a1), enabled(true) {}
971
972 bool is_enabled() const {
973 return enabled;
974 };
975
976 void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
977 el.AddMember("type", "IdealHelmholtzLogTau", doc.GetAllocator());
978 el.AddMember("a1", static_cast<double>(a1), doc.GetAllocator());
979 };
980 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
981};
982
988class IdealHelmholtzPower : public BaseHelmholtzTerm
989{
990
991 private:
992 std::vector<CoolPropDbl> n, t; // Use these variables internally
993 std::size_t N;
994 bool enabled;
995
996 public:
997 IdealHelmholtzPower() : N(0), enabled(false){};
998 // Constructor
999 IdealHelmholtzPower(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& t) : n(n), t(t), N(n.size()), enabled(true){};
1000
1001 bool is_enabled() const {
1002 return enabled;
1003 };
1004
1005 void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
1006 el.AddMember("type", "IdealHelmholtzPower", doc.GetAllocator());
1007 cpjson::set_long_double_array("n", n, el, doc);
1008 cpjson::set_long_double_array("t", t, el, doc);
1009 };
1010 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1011};
1012
1059class IdealHelmholtzPlanckEinsteinGeneralized : public BaseHelmholtzTerm
1060{
1061
1062 private:
1063 std::vector<CoolPropDbl> n, theta, c, d; // Use these variables internally
1064 std::size_t N;
1065 bool enabled;
1066
1067 public:
1068 IdealHelmholtzPlanckEinsteinGeneralized() : N(0), enabled(false) {}
1069 // Constructor with std::vector instances
1070 IdealHelmholtzPlanckEinsteinGeneralized(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta,
1071 const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& d)
1072 : n(n), theta(theta), c(c), d(d), N(n.size()), enabled(true) {}
1073
1074 // Extend the vectors to allow for multiple instances feeding values to this function
1075 void extend(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, const std::vector<CoolPropDbl>& c,
1076 const std::vector<CoolPropDbl>& d) {
1077 this->n.insert(this->n.end(), n.begin(), n.end());
1078 this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1079 this->c.insert(this->c.end(), c.begin(), c.end());
1080 this->d.insert(this->d.end(), d.begin(), d.end());
1081 N += n.size();
1082 }
1083
1084 bool is_enabled() const {
1085 return enabled;
1086 };
1087
1088 void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
1089 el.AddMember("type", "IdealHelmholtzPlanckEinsteinGeneralized", doc.GetAllocator());
1090 cpjson::set_long_double_array("n", n, el, doc);
1091 cpjson::set_long_double_array("theta", theta, el, doc);
1092 };
1093 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1094};
1095
1096class IdealHelmholtzCP0Constant : public BaseHelmholtzTerm
1097{
1098
1099 private:
1100 double cp_over_R, Tc, T0, tau0; // Use these variables internally
1101 bool enabled;
1102
1103 public:
1105 IdealHelmholtzCP0Constant() : cp_over_R(_HUGE), Tc(_HUGE), T0(_HUGE), tau0(_HUGE) {
1106 enabled = false;
1107 };
1108
1110 IdealHelmholtzCP0Constant(CoolPropDbl cp_over_R, CoolPropDbl Tc, CoolPropDbl T0) : cp_over_R(cp_over_R), Tc(Tc), T0(T0) {
1111 enabled = true;
1112 tau0 = Tc / T0;
1113 };
1114
1116 ~IdealHelmholtzCP0Constant(){};
1117
1118 bool is_enabled() const {
1119 return enabled;
1120 };
1121
1122 void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
1123 el.AddMember("type", "IdealGasHelmholtzCP0Constant", doc.GetAllocator());
1124 el.AddMember("cp_over_R", cp_over_R, doc.GetAllocator());
1125 el.AddMember("Tc", Tc, doc.GetAllocator());
1126 el.AddMember("T0", T0, doc.GetAllocator());
1127 };
1128
1129 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1130};
1131
1132class IdealHelmholtzCP0PolyT : public BaseHelmholtzTerm
1133{
1134 private:
1135 std::vector<CoolPropDbl> c, t;
1136 CoolPropDbl Tc, T0, tau0; // Use these variables internally
1137 std::size_t N;
1138 bool enabled;
1139
1140 public:
1141 IdealHelmholtzCP0PolyT() : Tc(_HUGE), T0(_HUGE), tau0(_HUGE), N(0), enabled(false) {}
1142
1144 IdealHelmholtzCP0PolyT(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t, double Tc, double T0)
1145 : c(c), t(t), Tc(Tc), T0(T0), tau0(Tc / T0), N(c.size()), enabled(true) {
1146 assert(c.size() == t.size());
1147 }
1148
1149 void extend(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t) {
1150 this->c.insert(this->c.end(), c.begin(), c.end());
1151 this->t.insert(this->t.end(), t.begin(), t.end());
1152 N += c.size();
1153 }
1154
1155 bool is_enabled() const {
1156 return enabled;
1157 };
1158
1159 void to_json(rapidjson::Value& el, rapidjson::Document& doc);
1160 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
1161#if ENABLE_CATCH
1162 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
1163#endif
1164};
1168class IdealHelmholtzGERG2004Sinh : public BaseHelmholtzTerm
1169{
1170 private:
1171 std::vector<CoolPropDbl> n, theta;
1172 CoolPropDbl Tc, _Tr;
1173 std::size_t N;
1174 bool enabled;
1175
1176 public:
1177 IdealHelmholtzGERG2004Sinh() : Tc(_HUGE), _Tr(_HUGE), N(0), enabled(false) {}
1178
1180 IdealHelmholtzGERG2004Sinh(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, double Tc)
1181 : n(n), theta(theta), Tc(Tc), _Tr(_HUGE), N(n.size()), enabled(true) {
1182 assert(n.size() == theta.size());
1183 }
1184
1185 void extend(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t) {
1186 this->n.insert(this->n.end(), n.begin(), n.end());
1187 this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1188 N += c.size();
1189 }
1190 void set_Tred(CoolPropDbl Tr) {
1191 this->_Tr = Tr;
1192 }
1193
1194 bool is_enabled() const {
1195 return enabled;
1196 };
1197 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
1198
1199 #if ENABLE_CATCH
1200 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
1201 #endif
1202};
1203
1204class IdealHelmholtzGERG2004Cosh : public BaseHelmholtzTerm
1205{
1206 private:
1207 std::vector<CoolPropDbl> n, theta;
1208 CoolPropDbl Tc, _Tr;
1209 std::size_t N;
1210 bool enabled;
1211
1212 public:
1213 IdealHelmholtzGERG2004Cosh() : Tc(_HUGE), _Tr(_HUGE), N(0), enabled(false) {}
1214
1216 IdealHelmholtzGERG2004Cosh(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, double Tc)
1217 : n(n), theta(theta), Tc(Tc), _Tr(_HUGE), N(n.size()), enabled(true) {
1218 assert(n.size() == theta.size());
1219 }
1220
1221 void extend(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta) {
1222 this->n.insert(this->n.end(), n.begin(), n.end());
1223 this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1224 N += n.size();
1225 }
1226 void set_Tred(CoolPropDbl Tr) {
1227 this->_Tr = Tr;
1228 }
1229
1230 bool is_enabled() const {
1231 return enabled;
1232 };
1233 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
1234
1235 #if ENABLE_CATCH
1236 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
1237 #endif
1238
1239};
1240
1243//\f[
1244//\frac{c_p^0}{R_u} = A + B\left(\frac{C/T}{\sinh(C/T)}\right)^2 + D\left(\frac{E/T}{\cosh(E/T)}\right)^2
1245//\f]
1246//Second partial of ideal-gas Helmholtz energy given directly by specific heat (\f$\displaystyle\alpha_{\tau\tau}^0=-\frac{1}{\tau^2}\frac{c_p^0}{R_u} \f$) - this is obtained by real gas \f$c_p\f$ relationship, and killing off residual Helmholtz terms
1247//\f[
1248//\alpha^0_{\tau\tau} = -\frac{A}{\tau^2} - \frac{B}{\tau^2}\left(\frac{C/T}{\sinh(C/T)}\right)^2 - \frac{D}{\tau^2}\left(\frac{E/T}{\cosh(E/T)}\right)^2
1249//\f]
1250//or in terms of \f$ \tau \f$:
1251//\f[
1252//\alpha^0_{\tau\tau} = -\frac{A}{\tau^2} - \frac{BC^2}{T_c^2}\left(\frac{1}{\sinh(C\tau/T_c)}\right)^2 - \frac{DE^2}{T_c^2}\left(\frac{1}{\cosh(E\tau/T_c)}\right)^2
1253//\f]
1254//Third partial:
1255//\f[
1256//\alpha^0_{\tau\tau\tau} = 2\frac{A}{\tau^3} + 2\frac{BC^3}{T_c^3}\frac{\cosh(C\tau/T_c)}{\sinh^3(C\tau/T_c)} +2 \frac{DE^3}{T_c^3}\frac{\sinh(E\tau/T_c)}{\cosh^3(E\tau/T_c)}
1257//\f]
1258//Now coming back to the ideal gas Helmholtz energy definition:
1259//\f[
1260//\alpha^0 = -\tau\displaystyle\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau+\displaystyle\int_{\tau_0}^{\tau} \frac{1}{\tau}\frac{c_p^0}{R_u}d\tau
1261//\f]
1262//Applying derivative
1263//\f[
1264//\alpha^0_{\tau} = -\displaystyle\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau-\tau\frac{\partial}{\partial \tau}\left[\displaystyle\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau \right]+\frac{\partial}{\partial \tau}\left[\displaystyle\int_{\tau_0}^{\tau} \frac{1}{\tau}\frac{c_p^0}{R_u}d\tau \right]
1265//\f]
1266//Fundamental theorem of calculus
1267//\f[
1268//\alpha^0_{\tau} = -\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau-\tau \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau+\frac{1}{\tau}\frac{c_p^0}{R_u}
1269//\f]
1270//Last two terms cancel, leaving
1271//\f[
1272//\alpha^0_{\tau} = -\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau
1273//\f]
1274//Another derivative yields (from fundamental theorem of calculus)
1275//\f[
1276//\alpha^0_{\tau\tau} = - \frac{1}{\tau^2}\frac{c_p^0}{R_u}
1277//\f]
1278//
1279//see also Jaeschke and Schley, 1995, (http://link.springer.com/article/10.1007%2FBF02083547#page-1)
1280//*/
1282//class IdealHelmholtzCP0AlyLee : public BaseHelmholtzTerm{
1283//private:
1284// std::vector<CoolPropDbl> c;
1285// CoolPropDbl Tc, tau0, T0; // Use these variables internally
1286// bool enabled;
1287//public:
1288// IdealHelmholtzCP0AlyLee(){enabled = false;};
1289//
1290// /// Constructor with std::vectors
1291// IdealHelmholtzCP0AlyLee(const std::vector<CoolPropDbl> &c, double Tc, double T0)
1292// :c(c), Tc(Tc), T0(T0)
1293// {
1294// tau0=Tc/T0;
1295// enabled = true;
1296// };
1297//
1298// /// Destructor
1299// ~IdealHelmholtzCP0AlyLee(){};
1300//
1301// bool is_enabled() const {return enabled;};
1302//
1303// void to_json(rapidjson::Value &el, rapidjson::Document &doc);
1304//
1305//
1306// /// The antiderivative given by \f$ \displaystyle\int \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau \f$
1307// /**
1308// sympy code for this derivative:
1309//
1310// from sympy import *
1311// a1,a2,a3,a4,a5,Tc,tau = symbols('a1,a2,a3,a4,a5,Tc,tau', real = True)
1312// integrand = a1 + a2*(a3/Tc/sinh(a3*tau/Tc))**2 + a4*(a5/Tc/cosh(a5*tau/Tc))**2
1313// integrand = integrand.rewrite(exp)
1314// antideriv = trigsimp(integrate(integrand,tau))
1315// display(antideriv)
1316// print latex(antideriv)
1317// print ccode(antideriv)
1318//
1319// \f[
1320// \displaystyle\int \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau = -\frac{a_0}{\tau}+\frac{2a_1a_2}{T_c\left[\exp\left(-\frac{2a_2\tau}{T_c}\right)-1\right]}+\frac{2a_3a_4}{T_c\left[\exp\left(-\frac{2a_4\tau}{T_c}\right)+1\right]}
1321// \f]
1322// */
1323// CoolPropDbl anti_deriv_cp0_tau2(const CoolPropDbl &tau);
1324//
1325// /// The antiderivative given by \f$ \displaystyle\int \frac{1}{\tau}\frac{c_p^0}{R_u}d\tau \f$
1326// /**
1327// sympy code for this derivative:
1328//
1329// a_0,a_1,a_2,a_3,a_4,Tc,tau = symbols('a_0,a_1,a_2,a_3,a_4,Tc,tau', real = True)
1330// integrand = a_0/tau + a_1/tau*(a_2*tau/Tc/sinh(a_2*tau/Tc))**2 + a_3/tau*(a_4*tau/Tc/cosh(a_4*tau/Tc))**2
1331//
1332// term2 = a_1/tau*(a_2*tau/Tc/sinh(a_2*tau/Tc))**2
1333// term2 = term2.rewrite(exp) # Unpack the sinh to exp functions
1334// antideriv2 = trigsimp(integrate(term2,tau))
1335// display(antideriv2)
1336// print latex(antideriv2)
1337// print ccode(antideriv2)
1338//
1339// term3 = a_3/tau*(a_4*tau/Tc/cosh(a_4*tau/Tc))**2
1340// term3 = term3.rewrite(exp) # Unpack the cosh to exp functions
1341// antideriv3 = factor(trigsimp(integrate(term3,tau).rewrite(exp)))
1342// display(antideriv3)
1343// print latex(antideriv3)
1344// print ccode(antideriv3)
1345//
1346// Can be broken into three parts (trick is to express \f$sinh\f$ and \f$cosh\f$ in terms of \f$exp\f$ function)
1347//
1348// Term 2:
1349// \f[
1350// \displaystyle\int \frac{a_1a_2^2}{T_c^2}\frac{\tau}{\sinh\left(\displaystyle\frac{a_2\tau}{T_c}\right)^2} d\tau = \frac{2 a_{1} a_{2} \tau}{- Tc + Tc e^{- \frac{2 a_{2}}{Tc} \tau}} + a_{1} \log{\left (-1 + e^{- \frac{2 a_{2}}{Tc} \tau} \right )} + \frac{2 a_{1}}{Tc} a_{2} \tau
1351// \f]
1352//
1353// Term 3:
1354// \f[
1355// \displaystyle\int \frac{a_1a_2^2}{T_c^2}\frac{\tau}{\cosh\left(\displaystyle\frac{a_2\tau}{T_c}\right)^2} d\tau = - \frac{a_{3}}{Tc \left(e^{\frac{2 a_{4}}{Tc} \tau} + 1\right)} \left(Tc e^{\frac{2 a_{4}}{Tc} \tau} \log{\left (e^{\frac{2 a_{4}}{Tc} \tau} + 1 \right )} + Tc \log{\left (e^{\frac{2 a_{4}}{Tc} \tau} + 1 \right )} - 2 a_{4} \tau e^{\frac{2 a_{4}}{Tc} \tau}\right)
1356// \f]
1357// */
1358// CoolPropDbl anti_deriv_cp0_tau(const CoolPropDbl &tau);
1359//
1360// CoolPropDbl base(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1361// CoolPropDbl dDelta(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1362// CoolPropDbl dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1363// CoolPropDbl dDelta2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1364// CoolPropDbl dDelta_dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1365// CoolPropDbl dTau2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1366// CoolPropDbl dDelta3(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1367// CoolPropDbl dDelta2_dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1368// CoolPropDbl dDelta_dTau2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1369// CoolPropDbl dTau3(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1370// CoolPropDbl dTau4(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1371//
1372//};
1373
1374class IdealHelmholtzContainer : public BaseHelmholtzContainer
1375{
1376 private:
1377 double _prefactor;
1378
1379 public:
1380 IdealHelmholtzLead Lead;
1381 IdealHelmholtzEnthalpyEntropyOffset EnthalpyEntropyOffsetCore, EnthalpyEntropyOffset;
1382 IdealHelmholtzLogTau LogTau;
1383 IdealHelmholtzPower Power;
1384 IdealHelmholtzPlanckEinsteinGeneralized PlanckEinstein;
1385
1386 IdealHelmholtzCP0Constant CP0Constant;
1387 IdealHelmholtzCP0PolyT CP0PolyT;
1388 IdealHelmholtzGERG2004Cosh GERG2004Cosh;
1389 IdealHelmholtzGERG2004Sinh GERG2004Sinh;
1390
1391 IdealHelmholtzContainer() : _prefactor(1.0){};
1392
1393 void set_prefactor(double prefactor) {
1394 _prefactor = prefactor;
1395 }
1396
1397 void set_Tred(double T_red) {
1398 GERG2004Cosh.set_Tred(T_red);
1399 GERG2004Sinh.set_Tred(T_red);
1400 }
1401
1402 void empty_the_EOS() {
1403 Lead = IdealHelmholtzLead();
1404 EnthalpyEntropyOffsetCore = IdealHelmholtzEnthalpyEntropyOffset();
1405 EnthalpyEntropyOffset = IdealHelmholtzEnthalpyEntropyOffset();
1406 LogTau = IdealHelmholtzLogTau();
1407 Power = IdealHelmholtzPower();
1408 PlanckEinstein = IdealHelmholtzPlanckEinsteinGeneralized();
1409 CP0Constant = IdealHelmholtzCP0Constant();
1410 CP0PolyT = IdealHelmholtzCP0PolyT();
1411 GERG2004Cosh = IdealHelmholtzGERG2004Cosh();
1412 GERG2004Sinh = IdealHelmholtzGERG2004Sinh();
1413 };
1414
1415 HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values = false) {
1416 HelmholtzDerivatives derivs; // zeros out the elements
1417 Lead.all(tau, delta, derivs);
1418 EnthalpyEntropyOffsetCore.all(tau, delta, derivs);
1419 EnthalpyEntropyOffset.all(tau, delta, derivs);
1420 LogTau.all(tau, delta, derivs);
1421 Power.all(tau, delta, derivs);
1422 PlanckEinstein.all(tau, delta, derivs);
1423 CP0Constant.all(tau, delta, derivs);
1424 CP0PolyT.all(tau, delta, derivs);
1425 GERG2004Cosh.all(tau, delta, derivs);
1426 GERG2004Sinh.all(tau, delta, derivs);
1427
1428 if (cache_values) {
1429 _base = derivs.alphar * _prefactor;
1430 _dDelta = derivs.dalphar_ddelta * _prefactor;
1431 _dTau = derivs.dalphar_dtau * _prefactor;
1432 _dDelta2 = derivs.d2alphar_ddelta2 * _prefactor;
1433 _dTau2 = derivs.d2alphar_dtau2 * _prefactor;
1434 _dDelta_dTau = derivs.d2alphar_ddelta_dtau * _prefactor;
1435 _dDelta3 = derivs.d3alphar_ddelta3 * _prefactor;
1436 _dTau3 = derivs.d3alphar_dtau3 * _prefactor;
1437 _dDelta2_dTau = derivs.d3alphar_ddelta2_dtau * _prefactor;
1438 _dDelta_dTau2 = derivs.d3alphar_ddelta_dtau2 * _prefactor;
1439 }
1440 return derivs * _prefactor;
1441 };
1442};
1443}; /* namespace CoolProp */
1444
1445#endif
Generated on Fri Apr 25 2025 03:25:14 for CoolProp by doxygen 1.9.4