CoolProp 8.0.0
An open-source fluid property and humid air property database
Helmholtz.h
Go to the documentation of this file.
1
2#ifndef HELMHOLTZ_H
3#define HELMHOLTZ_H
4
5#include <algorithm>
6#include <array>
7#include <cassert>
8#include <cstring> // for std::memset, used below -- libstdc++ does not pull it in transitively
9#include <vector>
10#include "CoolProp/detail/tools.h" // for CoolPropDbl
11//#include "Eigen/Core"
12#include "time.h"
13#include "CoolProp/detail/CachedElement.h"
14#include <memory>
15using std::shared_ptr;
16#include "CoolProp/numerics/numerics.h"
17
18// Forward declaration of AbstractCubic — only used via shared_ptr<AbstractCubic>
19// in this header, so the full type from Backends/Cubics/GeneralizedCubic.h is not
20// required at parse time. Note: AbstractCubic is declared in the global namespace
21// (not CoolProp::) by GeneralizedCubic.h, so this forward decl must match. TUs
22// that actually call methods on the pointer must include
23// "Backends/Cubics/GeneralizedCubic.h" directly.
24class AbstractCubic;
25
26#if ENABLE_CATCH
27# include "MultiComplex/MultiComplex.hpp"
28#endif
29
30namespace CoolProp {
31
32// #############################################################################
33// #############################################################################
34// #############################################################################
35// RESIDUAL TERMS
36// #############################################################################
37// #############################################################################
38// #############################################################################
39
40#define LIST_OF_DERIVATIVE_VARIABLES \
41 X(alphar) \
42 X(dalphar_ddelta) \
43 X(dalphar_dtau) \
44 X(d2alphar_ddelta2) \
45 X(d2alphar_dtau2) \
46 X(d2alphar_ddelta_dtau) \
47 X(d3alphar_ddelta3) \
48 X(d3alphar_ddelta_dtau2) \
49 X(d3alphar_ddelta2_dtau) \
50 X(d3alphar_dtau3) \
51 X(d4alphar_ddelta4) \
52 X(d4alphar_ddelta3_dtau) \
53 X(d4alphar_ddelta2_dtau2) \
54 X(d4alphar_ddelta_dtau3) \
55 X(d4alphar_dtau4) \
56 X(delta_x_dalphar_ddelta) \
57 X(tau_x_dalphar_dtau) \
58 X(delta2_x_d2alphar_ddelta2) \
59 X(deltatau_x_d2alphar_ddelta_dtau) \
60 X(tau2_x_d2alphar_dtau2)
61
62struct HelmholtzDerivatives
63{
64#define X(name) CoolPropDbl name;
65 LIST_OF_DERIVATIVE_VARIABLES
66#undef X
67 CoolPropDbl tau, delta, T_red, rhomolar_red;
68
69 void reset(CoolPropDbl v){
70#define X(name) name = v;
71 LIST_OF_DERIVATIVE_VARIABLES
72#undef X
73 } HelmholtzDerivatives operator+(const HelmholtzDerivatives& other) const {
74 HelmholtzDerivatives _new;
75#define X(name) _new.name = name + other.name;
76 LIST_OF_DERIVATIVE_VARIABLES
77#undef X
78 return _new;
79 }
80 HelmholtzDerivatives operator*(const CoolPropDbl& other) const {
81 HelmholtzDerivatives _new;
82#define X(name) _new.name = name * other;
83 LIST_OF_DERIVATIVE_VARIABLES
84#undef X
85 return _new;
86 }
87 HelmholtzDerivatives() : tau(0.0), delta(0.0), T_red(_HUGE), rhomolar_red(_HUGE) {
88 reset(0.0);
89 };
91 double get(std::size_t itau, std::size_t idelta) {
92 if (itau == 0) {
93 if (idelta == 0) {
94 return alphar;
95 } else if (idelta == 1) {
96 return dalphar_ddelta;
97 } else if (idelta == 2) {
98 return d2alphar_ddelta2;
99 } else if (idelta == 3) {
100 return d3alphar_ddelta3;
101 } else if (idelta == 4) {
102 return d4alphar_ddelta4;
103 } else {
104 throw ValueError();
105 }
106 } else if (itau == 1) {
107 if (idelta == 0) {
108 return dalphar_dtau;
109 } else if (idelta == 1) {
110 return d2alphar_ddelta_dtau;
111 } else if (idelta == 2) {
112 return d3alphar_ddelta2_dtau;
113 } else if (idelta == 3) {
114 return d4alphar_ddelta3_dtau;
115 } else {
116 throw ValueError();
117 }
118 } else if (itau == 2) {
119 if (idelta == 0) {
120 return d2alphar_dtau2;
121 } else if (idelta == 1) {
122 return d3alphar_ddelta_dtau2;
123 } else if (idelta == 2) {
124 return d4alphar_ddelta2_dtau2;
125 } else {
126 throw ValueError();
127 }
128 } else if (itau == 3) {
129 if (idelta == 0) {
130 return d3alphar_dtau3;
131 } else if (idelta == 1) {
132 return d4alphar_ddelta_dtau3;
133 } else {
134 throw ValueError();
135 }
136 } else if (itau == 4) {
137 if (idelta == 0) {
138 return d4alphar_dtau4;
139 } else {
140 throw ValueError();
141 }
142 } else {
143 throw ValueError();
144 }
145 }
146};
147#undef LIST_OF_DERIVATIVE_VARIABLES
148
150
175class BaseHelmholtzTerm
176{
177 public:
178 BaseHelmholtzTerm() = default;
179 virtual ~BaseHelmholtzTerm() = default;
180
182
185 virtual CoolPropDbl base(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
186 HelmholtzDerivatives deriv;
187 all(tau, delta, deriv);
188 return deriv.alphar;
189 };
191
194 virtual CoolPropDbl dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
195 HelmholtzDerivatives deriv;
196 all(tau, delta, deriv);
197 return deriv.dalphar_dtau;
198 };
200
203 virtual CoolPropDbl dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
204 HelmholtzDerivatives deriv;
205 all(tau, delta, deriv);
206 return deriv.d2alphar_dtau2;
207 };
209
212 virtual CoolPropDbl dDelta_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
213 HelmholtzDerivatives deriv;
214 all(tau, delta, deriv);
215 return deriv.d2alphar_ddelta_dtau;
216 };
218
221 virtual CoolPropDbl dDelta(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
222 HelmholtzDerivatives deriv;
223 all(tau, delta, deriv);
224 return deriv.dalphar_ddelta;
225 };
227
230 virtual CoolPropDbl dDelta2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
231 HelmholtzDerivatives deriv;
232 all(tau, delta, deriv);
233 return deriv.d2alphar_ddelta2;
234 };
236
239 virtual CoolPropDbl dDelta2_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
240 HelmholtzDerivatives deriv;
241 all(tau, delta, deriv);
242 return deriv.d3alphar_ddelta2_dtau;
243 };
245
248 virtual CoolPropDbl dDelta_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
249 HelmholtzDerivatives deriv;
250 all(tau, delta, deriv);
251 return deriv.d3alphar_ddelta_dtau2;
252 };
254
257 virtual CoolPropDbl dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
258 HelmholtzDerivatives deriv;
259 all(tau, delta, deriv);
260 return deriv.d3alphar_dtau3;
261 };
263
266 virtual CoolPropDbl dDelta3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
267 HelmholtzDerivatives deriv;
268 all(tau, delta, deriv);
269 return deriv.d3alphar_ddelta3;
270 };
272
275 virtual CoolPropDbl dTau4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
276 HelmholtzDerivatives deriv;
277 all(tau, delta, deriv);
278 return deriv.d4alphar_dtau4;
279 };
280 virtual CoolPropDbl dDelta_dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
281 HelmholtzDerivatives deriv;
282 all(tau, delta, deriv);
283 return deriv.d4alphar_ddelta_dtau3;
284 };
285 virtual CoolPropDbl dDelta2_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
286 HelmholtzDerivatives deriv;
287 all(tau, delta, deriv);
288 return deriv.d4alphar_ddelta2_dtau2;
289 };
290 virtual CoolPropDbl dDelta3_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
291 HelmholtzDerivatives deriv;
292 all(tau, delta, deriv);
293 return deriv.d4alphar_ddelta3_dtau;
294 };
295 virtual CoolPropDbl dDelta4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
296 HelmholtzDerivatives deriv;
297 all(tau, delta, deriv);
298 return deriv.d4alphar_ddelta4;
299 };
300
301 virtual void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) = 0;
302#if ENABLE_CATCH
303 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const {
304 throw CoolProp::NotImplementedError("The mcx derivative function was not implemented");
305 }
306#endif
307};
308
309struct ResidualHelmholtzGeneralizedExponentialElement
310{
312 CoolPropDbl n, d, t;
315 CoolPropDbl c, l_double, omega, m_double, eta1, epsilon1, eta2, epsilon2, beta1, gamma1, beta2, gamma2;
317 int l_int, m_int;
319 bool l_is_int, m_is_int;
320
321 ResidualHelmholtzGeneralizedExponentialElement()
322 : n(0),
323 d(0),
324 t(0),
325 c(0),
326 l_double(0),
327 omega(0),
328 m_double(0),
329 eta1(0),
330 epsilon1(0),
331 eta2(0),
332 epsilon2(0),
333 beta1(0),
334 gamma1(0),
335 beta2(0),
336 gamma2(0),
337 l_int(0),
338 m_int(0),
339 l_is_int(false),
340 m_is_int(true) {
341
342 };
343};
352class ResidualHelmholtzGeneralizedExponential : public BaseHelmholtzTerm
353{
354
355 public:
356 bool delta_li_in_u, tau_mi_in_u, eta1_in_u, eta2_in_u, beta1_in_u, beta2_in_u, finished;
357 std::vector<CoolPropDbl> s;
358 std::size_t N;
359
360 // These variables are for the exp(u) part
361 // 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
362 std::vector<double> n, d, t, c, l_double, omega, m_double, eta1, epsilon1, eta2, epsilon2, beta1, gamma1, beta2, gamma2;
363 // 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)
364 std::vector<int> l_int, m_int;
365
366 //Eigen::ArrayXd uE, du_ddeltaE, du_dtauE, d2u_ddelta2E, d2u_dtau2E, d3u_ddelta3E, d3u_dtau3E;
367
368 std::vector<ResidualHelmholtzGeneralizedExponentialElement> elements;
369 // Default Constructor
370 ResidualHelmholtzGeneralizedExponential()
371 : 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) {};
377 void add_Power(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
378 const std::vector<CoolPropDbl>& l) {
379 for (std::size_t i = 0; i < n.size(); ++i) {
380 ResidualHelmholtzGeneralizedExponentialElement el;
381 el.n = n[i];
382 el.d = d[i];
383 el.t = t[i];
384 el.l_double = l[i];
385 el.l_int = (int)el.l_double;
386 if (el.l_double > 0)
387 el.c = 1.0;
388 else
389 el.c = 0.0;
390 elements.push_back(el);
391 }
392 delta_li_in_u = true;
393 };
399 void add_Exponential(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
400 const std::vector<CoolPropDbl>& g, const std::vector<CoolPropDbl>& l) {
401 for (std::size_t i = 0; i < n.size(); ++i) {
402 ResidualHelmholtzGeneralizedExponentialElement el;
403 el.n = n[i];
404 el.d = d[i];
405 el.t = t[i];
406 el.c = g[i];
407 el.l_double = l[i];
408 el.l_int = (int)el.l_double;
409 elements.push_back(el);
410 }
411 delta_li_in_u = true;
412 }
418 void add_Gaussian(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
419 const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& beta,
420 const std::vector<CoolPropDbl>& gamma) {
421 for (std::size_t i = 0; i < n.size(); ++i) {
422 ResidualHelmholtzGeneralizedExponentialElement el;
423 el.n = n[i];
424 el.d = d[i];
425 el.t = t[i];
426 el.eta2 = eta[i];
427 el.epsilon2 = epsilon[i];
428 el.beta2 = beta[i];
429 el.gamma2 = gamma[i];
430 elements.push_back(el);
431 }
432 eta2_in_u = true;
433 beta2_in_u = true;
434 };
440 void add_GERG2008Gaussian(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
441 const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& beta,
442 const std::vector<CoolPropDbl>& gamma) {
443 for (std::size_t i = 0; i < n.size(); ++i) {
444 ResidualHelmholtzGeneralizedExponentialElement el;
445 el.n = n[i];
446 el.d = d[i];
447 el.t = t[i];
448 el.eta2 = eta[i];
449 el.epsilon2 = epsilon[i];
450 el.eta1 = beta[i];
451 el.epsilon1 = gamma[i];
452 elements.push_back(el);
453 }
454 eta2_in_u = true;
455 eta1_in_u = true;
456 };
462 void add_Lemmon2005(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
463 const std::vector<CoolPropDbl>& l, const std::vector<CoolPropDbl>& m) {
464 for (std::size_t i = 0; i < n.size(); ++i) {
465 ResidualHelmholtzGeneralizedExponentialElement el;
466 el.n = n[i];
467 el.d = d[i];
468 el.t = t[i];
469 el.c = 1.0;
470 el.omega = 1.0;
471 el.l_double = l[i];
472 el.m_double = m[i];
473 el.l_int = (int)el.l_double;
474 el.m_int = (int)el.m_double;
475 elements.push_back(el);
476 }
477 delta_li_in_u = true;
478 tau_mi_in_u = true;
479 };
485 void add_DoubleExponential(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
486 const std::vector<CoolPropDbl>& gd, const std::vector<CoolPropDbl>& ld, const std::vector<CoolPropDbl>& gt,
487 const std::vector<CoolPropDbl>& lt) {
488 for (std::size_t i = 0; i < n.size(); ++i) {
489 ResidualHelmholtzGeneralizedExponentialElement el;
490 el.n = n[i];
491 el.d = d[i];
492 el.t = t[i];
493 el.c = gd[i];
494 el.l_double = ld[i];
495 el.l_int = (int)el.l_double;
496 el.omega = gt[i];
497 el.m_double = lt[i];
498 elements.push_back(el);
499 }
500 delta_li_in_u = true;
501 tau_mi_in_u = true;
502 };
503
504 void finish() {
505 n.resize(elements.size());
506 d.resize(elements.size());
507 t.resize(elements.size());
508 c.resize(elements.size());
509 omega.resize(elements.size());
510 l_double.resize(elements.size());
511 l_int.resize(elements.size());
512 m_double.resize(elements.size());
513 m_int.resize(elements.size());
514 epsilon2.resize(elements.size());
515 eta2.resize(elements.size());
516 gamma2.resize(elements.size());
517 beta2.resize(elements.size());
518
519 for (std::size_t i = 0; i < elements.size(); ++i) {
520 n[i] = elements[i].n;
521 d[i] = elements[i].d;
522 t[i] = elements[i].t;
523 c[i] = elements[i].c;
524 omega[i] = elements[i].omega;
525 l_double[i] = elements[i].l_double;
526 l_int[i] = elements[i].l_int;
527 m_double[i] = elements[i].m_double;
528 m_int[i] = elements[i].m_int;
529 epsilon2[i] = elements[i].epsilon2;
530 eta2[i] = elements[i].eta2;
531 gamma2[i] = elements[i].gamma2;
532 beta2[i] = elements[i].beta2;
533
534 // See if l is an integer, and store a flag if it is
535 elements[i].l_is_int = (std::abs(static_cast<long>(elements[i].l_double) - elements[i].l_double) < 1e-14);
536 }
537 // uE.resize(elements.size());
538 // du_ddeltaE.resize(elements.size());
539 // du_dtauE.resize(elements.size());
540 // d2u_ddelta2E.resize(elements.size());
541 // d2u_dtau2E.resize(elements.size());
542 // d3u_ddelta3E.resize(elements.size());
543 // d3u_dtau3E.resize(elements.size());
544
545 finished = true;
546 };
547
548 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
549 //void allEigen(const CoolPropDbl &tau, const CoolPropDbl &delta, HelmholtzDerivatives &derivs) throw();
550
551#if ENABLE_CATCH
552 mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
553#endif
554};
555
556struct ResidualHelmholtzNonAnalyticElement
557{
558 CoolPropDbl n, a, b, beta, A, B, C, D;
559};
560class ResidualHelmholtzNonAnalytic : public BaseHelmholtzTerm
561{
562
563 public:
564 std::size_t N;
565 std::vector<CoolPropDbl> s;
566 std::vector<ResidualHelmholtzNonAnalyticElement> elements;
568 ResidualHelmholtzNonAnalytic()
569 : N(0) {
570
571 };
573 ~ResidualHelmholtzNonAnalytic() = default;
575 ResidualHelmholtzNonAnalytic(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& a, const std::vector<CoolPropDbl>& b,
576 const std::vector<CoolPropDbl>& beta, const std::vector<CoolPropDbl>& A, const std::vector<CoolPropDbl>& B,
577 const std::vector<CoolPropDbl>& C, const std::vector<CoolPropDbl>& D)
578 : N(n.size()) {
579
580 s.resize(N);
581 for (std::size_t i = 0; i < n.size(); ++i) {
582 ResidualHelmholtzNonAnalyticElement el;
583 el.n = n[i];
584 el.a = a[i];
585 el.b = b[i];
586 el.beta = beta[i];
587 el.A = A[i];
588 el.B = B[i];
589 el.C = C[i];
590 el.D = D[i];
591 elements.push_back(el);
592 }
593 };
594 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
595#if ENABLE_CATCH
596 mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
597#endif
598};
599
600class ResidualHelmholtzGeneralizedCubic : public BaseHelmholtzTerm
601{
602 protected:
603 shared_ptr<AbstractCubic> m_abstractcubic;
604 std::vector<double> z;
605 public:
606 bool enabled;
607
609 ResidualHelmholtzGeneralizedCubic()
610 : enabled(false) {
611
612 };
614 ResidualHelmholtzGeneralizedCubic(shared_ptr<AbstractCubic>& ac) : m_abstractcubic(ac), enabled(true) {
615
616 z = std::vector<double>(1, 1); // Init the vector to [1.0]
617 };
618
619 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() override;
620};
621
622class ResidualHelmholtzGaoB : public BaseHelmholtzTerm
623{
624 protected:
625 std::vector<double> n, t, d, eta, beta, gamma, epsilon, b;
626
627 public:
628 bool enabled;
629
631 ResidualHelmholtzGaoB()
632 : enabled(false) {
633
634 };
635
637 ResidualHelmholtzGaoB(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& t, const std::vector<CoolPropDbl>& d,
638 const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& beta, const std::vector<CoolPropDbl>& gamma,
639 const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& b)
640 : n(n), t(t), d(d), eta(eta), beta(beta), gamma(gamma), epsilon(epsilon), b(b), enabled(true) {
641
642 };
643
644 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
645
646#if ENABLE_CATCH
647 mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
648#endif
649};
650
652class ResidualHelmholtzXiangDeiters : public BaseHelmholtzTerm
653{
654
655 public:
656 bool enabled;
657 ResidualHelmholtzGeneralizedExponential phi0, phi1, phi2;
658 CoolPropDbl Tc, pc, rhomolarc, acentric, R, theta;
660 ResidualHelmholtzXiangDeiters()
661 : enabled(false), Tc(_HUGE), pc(_HUGE), rhomolarc(_HUGE), acentric(_HUGE), R(_HUGE), theta(_HUGE) {
662
663 };
665 ResidualHelmholtzXiangDeiters(const CoolPropDbl Tc, const CoolPropDbl pc, const CoolPropDbl rhomolarc, const CoolPropDbl acentric,
666 const CoolPropDbl R);
667 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
668#if ENABLE_CATCH
669 mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
670#endif
671};
672
673class ResidualHelmholtzSAFTAssociating : public BaseHelmholtzTerm
674{
675
676 protected:
677 double a, m, epsilonbar, vbarn, kappabar;
678
679 CoolPropDbl Deltabar(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
680 CoolPropDbl dDeltabar_ddelta__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
681 CoolPropDbl d2Deltabar_ddelta2__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
682 CoolPropDbl dDeltabar_dtau__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
683 CoolPropDbl d2Deltabar_dtau2__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
684 CoolPropDbl d2Deltabar_ddelta_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
685 CoolPropDbl d3Deltabar_dtau3__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
686 CoolPropDbl d3Deltabar_ddelta_dtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
687 CoolPropDbl d3Deltabar_ddelta3__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
688 CoolPropDbl d3Deltabar_ddelta2_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
689
690 CoolPropDbl X(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
691 CoolPropDbl dX_dDeltabar__constdelta(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
692 CoolPropDbl dX_ddelta__constDeltabar(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
693 CoolPropDbl dX_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
694 CoolPropDbl dX_ddelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
695 CoolPropDbl d2X_dtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
696 CoolPropDbl d2X_ddeltadtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
697 CoolPropDbl d2X_ddelta2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
698
699 CoolPropDbl d3X_dtau3(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
700 CoolPropDbl d3X_ddelta3(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
701 CoolPropDbl d3X_ddeltadtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
702 CoolPropDbl d3X_ddelta2dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
703
704 CoolPropDbl g(const CoolPropDbl& eta) const;
705 CoolPropDbl dg_deta(const CoolPropDbl& eta) const;
706 CoolPropDbl d2g_deta2(const CoolPropDbl& eta) const;
707 CoolPropDbl d3g_deta3(const CoolPropDbl& eta) const;
708 CoolPropDbl eta(const CoolPropDbl& delta) const;
709
710 public:
712 ResidualHelmholtzSAFTAssociating()
713 : a(_HUGE), m(_HUGE), epsilonbar(_HUGE), vbarn(_HUGE), kappabar(_HUGE), disabled(true) {
714
715 };
716
717 // Constructor
718 ResidualHelmholtzSAFTAssociating(double a, double m, double epsilonbar, double vbarn, double kappabar)
719 : a(a), m(m), epsilonbar(epsilonbar), vbarn(vbarn), kappabar(kappabar), disabled(false) {
720
721 };
722
723 bool disabled;
724
725 //Destructor. No Implementation
726 ~ResidualHelmholtzSAFTAssociating() = default;
727
728 CoolPropDbl dTau4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() override {
729 return 1e99;
730 };
731 CoolPropDbl dDelta_dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() override {
732 return 1e99;
733 };
734 CoolPropDbl dDelta2_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() override {
735 return 1e99;
736 };
737 CoolPropDbl dDelta3_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() override {
738 return 1e99;
739 };
740 CoolPropDbl dDelta4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() override {
741 return 1e99;
742 };
743
744 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& deriv) throw() override;
745};
746
747class BaseHelmholtzContainer
748{
749 protected:
750 std::array<double, 16> cache = create_filled_array<double, 16>(_HUGE);
751 std::array<bool, 16> is_cached = create_filled_array<bool, 16>(false);
752 constexpr static std::size_t i00 = 0, i01 = 1, i02 = 2, i03 = 3, i04 = 4, i10 = 5, i11 = 6, i12 = 7, i13 = 8, i20 = 9, i21 = 10, i22 = 11,
753 i30 = 12, i31 = 13, i40 = 14;
754
755 bool cache_valid(std::size_t i) const {
756 return is_cached[i];
757 }
758
759 public:
760 void clear() {
761 std::memset(cache.data(), 0, sizeof(cache));
762 std::memset(is_cached.data(), false, sizeof(is_cached));
763 };
764
765 virtual void empty_the_EOS() = 0;
766 virtual HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values) = 0;
767
768 CoolPropDbl base(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
769 if (!cache_valid(i00) || dont_use_cache)
770 return all(tau, delta, false).alphar;
771 else
772 return cache[i00];
773 };
774 CoolPropDbl dDelta(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
775 if (!cache_valid(i10) || dont_use_cache)
776 return all(tau, delta, false).dalphar_ddelta;
777 else
778 return cache[i10];
779 };
780 CoolPropDbl dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
781 if (!cache_valid(i01) || dont_use_cache)
782 return all(tau, delta, false).dalphar_dtau;
783 else
784 return cache[i01];
785 };
786 CoolPropDbl dDelta2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
787 if (!cache_valid(i20) || dont_use_cache)
788 return all(tau, delta, false).d2alphar_ddelta2;
789 else
790 return cache[i20];
791 };
792 CoolPropDbl dDelta_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
793 if (!cache_valid(i11) || dont_use_cache)
794 return all(tau, delta, false).d2alphar_ddelta_dtau;
795 else
796 return cache[i11];
797 };
798 CoolPropDbl dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
799 if (!cache_valid(i02) || dont_use_cache)
800 return all(tau, delta, false).d2alphar_dtau2;
801 else
802 return cache[i02];
803 };
804 CoolPropDbl dDelta3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
805 if (!cache_valid(i30) || dont_use_cache)
806 return all(tau, delta, false).d3alphar_ddelta3;
807 else
808 return cache[i30];
809 };
810 CoolPropDbl dDelta2_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
811 if (!cache_valid(i21) || dont_use_cache)
812 return all(tau, delta, false).d3alphar_ddelta2_dtau;
813 else
814 return cache[i21];
815 };
816 CoolPropDbl dDelta_dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
817 if (!cache_valid(i12) || dont_use_cache)
818 return all(tau, delta, false).d3alphar_ddelta_dtau2;
819 else
820 return cache[i12];
821 };
822 CoolPropDbl dTau3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
823 if (!cache_valid(i03) || dont_use_cache)
824 return all(tau, delta, false).d3alphar_dtau3;
825 else
826 return cache[i03];
827 };
828 CoolPropDbl dDelta4(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
829 return all(tau, delta, false).d4alphar_ddelta4;
830 };
831 CoolPropDbl dDelta3_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
832 return all(tau, delta, false).d4alphar_ddelta3_dtau;
833 };
834 CoolPropDbl dDelta2_dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
835 return all(tau, delta, false).d4alphar_ddelta2_dtau2;
836 };
837 CoolPropDbl dDelta_dTau3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
838 return all(tau, delta, false).d4alphar_ddelta_dtau3;
839 };
840 CoolPropDbl dTau4(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
841 return all(tau, delta, false).d4alphar_dtau4;
842 };
843};
844
845class ResidualHelmholtzContainer : public BaseHelmholtzContainer
846{
847 public:
848 ResidualHelmholtzNonAnalytic NonAnalytic;
849 ResidualHelmholtzSAFTAssociating SAFT;
850 ResidualHelmholtzGeneralizedExponential GenExp;
851 ResidualHelmholtzGeneralizedCubic cubic;
852 ResidualHelmholtzXiangDeiters XiangDeiters;
853 ResidualHelmholtzGaoB GaoB;
854
855 void empty_the_EOS() override {
856 NonAnalytic = ResidualHelmholtzNonAnalytic();
857 SAFT = ResidualHelmholtzSAFTAssociating();
858 GenExp = ResidualHelmholtzGeneralizedExponential();
859 cubic = ResidualHelmholtzGeneralizedCubic();
860 XiangDeiters = ResidualHelmholtzXiangDeiters();
861 GaoB = ResidualHelmholtzGaoB();
862 };
863
864 HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values = false) override {
865 HelmholtzDerivatives derivs; // zeros out the elements
866 GenExp.all(tau, delta, derivs);
867 NonAnalytic.all(tau, delta, derivs);
868 SAFT.all(tau, delta, derivs);
869 cubic.all(tau, delta, derivs);
870 XiangDeiters.all(tau, delta, derivs);
871 GaoB.all(tau, delta, derivs);
872 if (cache_values) {
873 cache[i00] = derivs.alphar;
874 cache[i10] = derivs.dalphar_ddelta;
875 cache[i01] = derivs.dalphar_dtau;
876 cache[i20] = derivs.d2alphar_ddelta2;
877 cache[i02] = derivs.d2alphar_dtau2;
878 cache[i11] = derivs.d2alphar_ddelta_dtau;
879 cache[i30] = derivs.d3alphar_ddelta3;
880 cache[i03] = derivs.d3alphar_dtau3;
881 cache[i21] = derivs.d3alphar_ddelta2_dtau;
882 cache[i12] = derivs.d3alphar_ddelta_dtau2;
883 std::memset(is_cached.data(), true, sizeof(is_cached));
884 }
885 return derivs;
886 };
887};
888
889// #############################################################################
890// #############################################################################
891// #############################################################################
892// IDEAL GAS TERMS
893// #############################################################################
894// #############################################################################
895// #############################################################################
896
898
903class IdealHelmholtzLead : public BaseHelmholtzTerm
904{
905
906 private:
907 CoolPropDbl a1, a2;
908 bool enabled;
909
910 public:
911 // Default constructor
912 IdealHelmholtzLead() : a1(_HUGE), a2(_HUGE), enabled(false) {}
913
914 // Constructor
915 IdealHelmholtzLead(CoolPropDbl a1, CoolPropDbl a2) : a1(a1), a2(a2), enabled(true) {}
916
917 bool is_enabled() const {
918 return enabled;
919 }
920
921 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() override;
922};
923
925
930class IdealHelmholtzEnthalpyEntropyOffset : public BaseHelmholtzTerm
931{
932 private:
933 CoolPropDbl a1, a2; // Use these variables internally
934 std::string reference;
935 bool enabled;
936
937 public:
938 IdealHelmholtzEnthalpyEntropyOffset() : a1(_HUGE), a2(_HUGE), enabled(false) {}
939
940 // Constructor
941 IdealHelmholtzEnthalpyEntropyOffset(CoolPropDbl a1, CoolPropDbl a2, const std::string& ref) : a1(a1), a2(a2), reference(ref), enabled(true) {}
942
943 // Set the values in the class
944 void set(CoolPropDbl a1, CoolPropDbl a2, const std::string& ref) {
945 // If it doesn't already exist, just set the values
946 if (enabled == false) {
947 this->a1 = a1;
948 this->a2 = a2;
949 enabled = true;
950 } else if (ref == "DEF") {
951 this->a1 = 0.0;
952 this->a2 = 0.0;
953 enabled = false;
954 } else {
955 // Otherwise, increment the values
956 this->a1 += a1;
957 this->a2 += a2;
958 enabled = true;
959 }
960 this->reference = ref;
961 }
962
963 bool is_enabled() const {
964 return enabled;
965 };
966
973 CoolPropDbl get_a1() const {
974 return a1;
975 }
976 CoolPropDbl get_a2() const {
977 return a2;
978 }
979
980 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() override;
981};
982
988class IdealHelmholtzLogTau : public BaseHelmholtzTerm
989{
990 private:
991 CoolPropDbl a1;
992 bool enabled;
993
994 public:
996 IdealHelmholtzLogTau() : a1(_HUGE), enabled(false) {}
997
998 // Constructor
999 IdealHelmholtzLogTau(CoolPropDbl a1) : a1(a1), enabled(true) {}
1000
1001 bool is_enabled() const {
1002 return enabled;
1003 };
1004
1005 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() override;
1006};
1007
1013class IdealHelmholtzPower : public BaseHelmholtzTerm
1014{
1015
1016 private:
1017 std::vector<CoolPropDbl> n, t; // Use these variables internally
1018 std::size_t N;
1019 bool enabled;
1020
1021 public:
1022 IdealHelmholtzPower() : N(0), enabled(false) {};
1023 // Constructor
1024 IdealHelmholtzPower(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& t) : n(n), t(t), N(n.size()), enabled(true) {};
1025
1026 bool is_enabled() const {
1027 return enabled;
1028 };
1029
1030 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() override;
1031};
1032
1079class IdealHelmholtzPlanckEinsteinGeneralized : public BaseHelmholtzTerm
1080{
1081
1082 private:
1083 std::vector<CoolPropDbl> n, theta, c, d; // Use these variables internally
1084 std::size_t N;
1085 bool enabled;
1086
1087 public:
1088 IdealHelmholtzPlanckEinsteinGeneralized() : N(0), enabled(false) {}
1089 // Constructor with std::vector instances
1090 IdealHelmholtzPlanckEinsteinGeneralized(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta,
1091 const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& d)
1092 : n(n), theta(theta), c(c), d(d), N(n.size()), enabled(true) {}
1093
1094 // Extend the vectors to allow for multiple instances feeding values to this function
1095 void extend(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, const std::vector<CoolPropDbl>& c,
1096 const std::vector<CoolPropDbl>& d) {
1097 this->n.insert(this->n.end(), n.begin(), n.end());
1098 this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1099 this->c.insert(this->c.end(), c.begin(), c.end());
1100 this->d.insert(this->d.end(), d.begin(), d.end());
1101 N += n.size();
1102 }
1103
1104 bool is_enabled() const {
1105 return enabled;
1106 };
1107
1108 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() override;
1109};
1110
1111class IdealHelmholtzCP0Constant : public BaseHelmholtzTerm
1112{
1113
1114 private:
1115 double cp_over_R, Tc, T0, tau0; // Use these variables internally
1116 bool enabled;
1117
1118 public:
1120 IdealHelmholtzCP0Constant()
1121 : cp_over_R(_HUGE), Tc(_HUGE), T0(_HUGE), tau0(_HUGE), enabled(false) {
1122
1123 };
1124
1126 IdealHelmholtzCP0Constant(CoolPropDbl cp_over_R, CoolPropDbl Tc, CoolPropDbl T0)
1127 : cp_over_R(cp_over_R), Tc(Tc), T0(T0), enabled(true), tau0(Tc / T0) {
1128
1129 };
1130
1132 ~IdealHelmholtzCP0Constant() = default;
1133
1134 bool is_enabled() const {
1135 return enabled;
1136 };
1137
1138 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() override;
1139};
1140
1141class IdealHelmholtzCP0PolyT : public BaseHelmholtzTerm
1142{
1143 private:
1144 std::vector<CoolPropDbl> c, t;
1145 CoolPropDbl Tc, T0, tau0; // Use these variables internally
1146 std::size_t N;
1147 bool enabled;
1148
1149 public:
1150 IdealHelmholtzCP0PolyT() : Tc(_HUGE), T0(_HUGE), tau0(_HUGE), N(0), enabled(false) {}
1151
1153 IdealHelmholtzCP0PolyT(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t, double Tc, double T0)
1154 : c(c), t(t), Tc(Tc), T0(T0), tau0(Tc / T0), N(c.size()), enabled(true) {
1155 assert(c.size() == t.size());
1156 }
1157
1158 void extend(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t) {
1159 this->c.insert(this->c.end(), c.begin(), c.end());
1160 this->t.insert(this->t.end(), t.begin(), t.end());
1161 N += c.size();
1162 }
1163
1164 bool is_enabled() const {
1165 return enabled;
1166 };
1167
1168 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
1169#if ENABLE_CATCH
1170 mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
1171#endif
1172};
1176class IdealHelmholtzGERG2004Sinh : public BaseHelmholtzTerm
1177{
1178 private:
1179 std::vector<CoolPropDbl> n, theta;
1180 CoolPropDbl Tc, _Tr;
1181 std::size_t N;
1182 bool enabled;
1183
1184 public:
1185 IdealHelmholtzGERG2004Sinh() : Tc(_HUGE), _Tr(_HUGE), N(0), enabled(false) {}
1186
1188 IdealHelmholtzGERG2004Sinh(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, double Tc)
1189 : n(n), theta(theta), Tc(Tc), _Tr(_HUGE), N(n.size()), enabled(true) {
1190 assert(n.size() == theta.size());
1191 }
1192
1193 void extend(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t) {
1194 this->n.insert(this->n.end(), n.begin(), n.end());
1195 this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1196 N += c.size();
1197 }
1198 void set_Tred(CoolPropDbl Tr) {
1199 this->_Tr = Tr;
1200 }
1201
1202 bool is_enabled() const {
1203 return enabled;
1204 };
1205 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
1206
1207#if ENABLE_CATCH
1208 mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
1209#endif
1210};
1211
1212class IdealHelmholtzGERG2004Cosh : public BaseHelmholtzTerm
1213{
1214 private:
1215 std::vector<CoolPropDbl> n, theta;
1216 CoolPropDbl Tc, _Tr;
1217 std::size_t N;
1218 bool enabled;
1219
1220 public:
1221 IdealHelmholtzGERG2004Cosh() : Tc(_HUGE), _Tr(_HUGE), N(0), enabled(false) {}
1222
1224 IdealHelmholtzGERG2004Cosh(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, double Tc)
1225 : n(n), theta(theta), Tc(Tc), _Tr(_HUGE), N(n.size()), enabled(true) {
1226 assert(n.size() == theta.size());
1227 }
1228
1229 void extend(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta) {
1230 this->n.insert(this->n.end(), n.begin(), n.end());
1231 this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1232 N += n.size();
1233 }
1234 void set_Tred(CoolPropDbl Tr) {
1235 this->_Tr = Tr;
1236 }
1237
1238 bool is_enabled() const {
1239 return enabled;
1240 };
1241 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
1242
1243#if ENABLE_CATCH
1244 mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
1245#endif
1246};
1247
1250//\f[
1251//\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
1252//\f]
1253//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
1254//\f[
1255//\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
1256//\f]
1257//or in terms of \f$ \tau \f$:
1258//\f[
1259//\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
1260//\f]
1261//Third partial:
1262//\f[
1263//\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)}
1264//\f]
1265//Now coming back to the ideal gas Helmholtz energy definition:
1266//\f[
1267//\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
1268//\f]
1269//Applying derivative
1270//\f[
1271//\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]
1272//\f]
1273//Fundamental theorem of calculus
1274//\f[
1275//\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}
1276//\f]
1277//Last two terms cancel, leaving
1278//\f[
1279//\alpha^0_{\tau} = -\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau
1280//\f]
1281//Another derivative yields (from fundamental theorem of calculus)
1282//\f[
1283//\alpha^0_{\tau\tau} = - \frac{1}{\tau^2}\frac{c_p^0}{R_u}
1284//\f]
1285//
1286//see also Jaeschke and Schley, 1995, (http://link.springer.com/article/10.1007%2FBF02083547#page-1)
1287//*/
1289//class IdealHelmholtzCP0AlyLee : public BaseHelmholtzTerm{
1290//private:
1291// std::vector<CoolPropDbl> c;
1292// CoolPropDbl Tc, tau0, T0; // Use these variables internally
1293// bool enabled;
1294//public:
1295// IdealHelmholtzCP0AlyLee(){enabled = false;};
1296//
1297// /// Constructor with std::vectors
1298// IdealHelmholtzCP0AlyLee(const std::vector<CoolPropDbl> &c, double Tc, double T0)
1299// :c(c), Tc(Tc), T0(T0)
1300// {
1301// tau0=Tc/T0;
1302// enabled = true;
1303// };
1304//
1305// /// Destructor
1306// ~IdealHelmholtzCP0AlyLee(){};
1307//
1308// bool is_enabled() const {return enabled;};
1309//
1310//
1311// /// The antiderivative given by \f$ \displaystyle\int \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau \f$
1312// /**
1313// sympy code for this derivative:
1314//
1315// from sympy import *
1316// a1,a2,a3,a4,a5,Tc,tau = symbols('a1,a2,a3,a4,a5,Tc,tau', real = True)
1317// integrand = a1 + a2*(a3/Tc/sinh(a3*tau/Tc))**2 + a4*(a5/Tc/cosh(a5*tau/Tc))**2
1318// integrand = integrand.rewrite(exp)
1319// antideriv = trigsimp(integrate(integrand,tau))
1320// display(antideriv)
1321// print latex(antideriv)
1322// print ccode(antideriv)
1323//
1324// \f[
1325// \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]}
1326// \f]
1327// */
1328// CoolPropDbl anti_deriv_cp0_tau2(const CoolPropDbl &tau);
1329//
1330// /// The antiderivative given by \f$ \displaystyle\int \frac{1}{\tau}\frac{c_p^0}{R_u}d\tau \f$
1331// /**
1332// sympy code for this derivative:
1333//
1334// 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)
1335// 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
1336//
1337// term2 = a_1/tau*(a_2*tau/Tc/sinh(a_2*tau/Tc))**2
1338// term2 = term2.rewrite(exp) # Unpack the sinh to exp functions
1339// antideriv2 = trigsimp(integrate(term2,tau))
1340// display(antideriv2)
1341// print latex(antideriv2)
1342// print ccode(antideriv2)
1343//
1344// term3 = a_3/tau*(a_4*tau/Tc/cosh(a_4*tau/Tc))**2
1345// term3 = term3.rewrite(exp) # Unpack the cosh to exp functions
1346// antideriv3 = factor(trigsimp(integrate(term3,tau).rewrite(exp)))
1347// display(antideriv3)
1348// print latex(antideriv3)
1349// print ccode(antideriv3)
1350//
1351// Can be broken into three parts (trick is to express \f$sinh\f$ and \f$cosh\f$ in terms of \f$exp\f$ function)
1352//
1353// Term 2:
1354// \f[
1355// \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
1356// \f]
1357//
1358// Term 3:
1359// \f[
1360// \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)
1361// \f]
1362// */
1363// CoolPropDbl anti_deriv_cp0_tau(const CoolPropDbl &tau);
1364//
1365// CoolPropDbl base(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1366// CoolPropDbl dDelta(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1367// CoolPropDbl dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1368// CoolPropDbl dDelta2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1369// CoolPropDbl dDelta_dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1370// CoolPropDbl dTau2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1371// CoolPropDbl dDelta3(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1372// CoolPropDbl dDelta2_dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1373// CoolPropDbl dDelta_dTau2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1374// CoolPropDbl dTau3(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1375// CoolPropDbl dTau4(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1376//
1377//};
1378
1379class IdealHelmholtzContainer : public BaseHelmholtzContainer
1380{
1381 private:
1382 double _prefactor;
1383
1384 public:
1385 IdealHelmholtzLead Lead;
1386 IdealHelmholtzEnthalpyEntropyOffset EnthalpyEntropyOffsetCore, EnthalpyEntropyOffset;
1387 IdealHelmholtzLogTau LogTau;
1388 IdealHelmholtzPower Power;
1389 IdealHelmholtzPlanckEinsteinGeneralized PlanckEinstein;
1390
1391 IdealHelmholtzCP0Constant CP0Constant;
1392 IdealHelmholtzCP0PolyT CP0PolyT;
1393 IdealHelmholtzGERG2004Cosh GERG2004Cosh;
1394 IdealHelmholtzGERG2004Sinh GERG2004Sinh;
1395
1396 IdealHelmholtzContainer() : _prefactor(1.0) {};
1397
1398 void set_prefactor(double prefactor) {
1399 _prefactor = prefactor;
1400 }
1401
1405 double get_prefactor() const {
1406 return _prefactor;
1407 }
1408
1409 void set_Tred(double T_red) {
1410 GERG2004Cosh.set_Tred(T_red);
1411 GERG2004Sinh.set_Tred(T_red);
1412 }
1413
1414 void empty_the_EOS() override {
1415 Lead = IdealHelmholtzLead();
1416 EnthalpyEntropyOffsetCore = IdealHelmholtzEnthalpyEntropyOffset();
1417 EnthalpyEntropyOffset = IdealHelmholtzEnthalpyEntropyOffset();
1418 LogTau = IdealHelmholtzLogTau();
1419 Power = IdealHelmholtzPower();
1420 PlanckEinstein = IdealHelmholtzPlanckEinsteinGeneralized();
1421 CP0Constant = IdealHelmholtzCP0Constant();
1422 CP0PolyT = IdealHelmholtzCP0PolyT();
1423 GERG2004Cosh = IdealHelmholtzGERG2004Cosh();
1424 GERG2004Sinh = IdealHelmholtzGERG2004Sinh();
1425 };
1426
1427 HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values = false) override {
1428 HelmholtzDerivatives derivs; // zeros out the elements
1429 Lead.all(tau, delta, derivs);
1430 EnthalpyEntropyOffsetCore.all(tau, delta, derivs);
1431 EnthalpyEntropyOffset.all(tau, delta, derivs);
1432 LogTau.all(tau, delta, derivs);
1433 Power.all(tau, delta, derivs);
1434 PlanckEinstein.all(tau, delta, derivs);
1435 CP0Constant.all(tau, delta, derivs);
1436 CP0PolyT.all(tau, delta, derivs);
1437 GERG2004Cosh.all(tau, delta, derivs);
1438 GERG2004Sinh.all(tau, delta, derivs);
1439
1440 if (cache_values) {
1441 cache[i00] = derivs.alphar * _prefactor;
1442 cache[i10] = derivs.dalphar_ddelta * _prefactor;
1443 cache[i01] = derivs.dalphar_dtau * _prefactor;
1444 cache[i20] = derivs.d2alphar_ddelta2 * _prefactor;
1445 cache[i02] = derivs.d2alphar_dtau2 * _prefactor;
1446 cache[i11] = derivs.d2alphar_ddelta_dtau * _prefactor;
1447 cache[i30] = derivs.d3alphar_ddelta3 * _prefactor;
1448 cache[i03] = derivs.d3alphar_dtau3 * _prefactor;
1449 cache[i21] = derivs.d3alphar_ddelta2_dtau * _prefactor;
1450 cache[i12] = derivs.d3alphar_ddelta_dtau2 * _prefactor;
1451 std::memset(is_cached.data(), true, sizeof(is_cached));
1452 }
1453 return derivs * _prefactor;
1454 };
1455};
1456}; /* namespace CoolProp */
1457
1458#endif
Generated by doxygen 1.9.4