CoolProp 7.2.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 <vector>
8#include "rapidjson_include.h"
9//#include "Eigen/Core"
10#include "time.h"
11#include "CachedElement.h"
12#include "Backends/Cubics/GeneralizedCubic.h"
13#include "crossplatform_shared_ptr.h"
14#include "CPnumerics.h"
15
16#if ENABLE_CATCH
17# include "MultiComplex/MultiComplex.hpp"
18#endif
19
20namespace CoolProp {
21
22// #############################################################################
23// #############################################################################
24// #############################################################################
25// RESIDUAL TERMS
26// #############################################################################
27// #############################################################################
28// #############################################################################
29
30#define LIST_OF_DERIVATIVE_VARIABLES \
31 X(alphar) \
32 X(dalphar_ddelta) \
33 X(dalphar_dtau) \
34 X(d2alphar_ddelta2) \
35 X(d2alphar_dtau2) \
36 X(d2alphar_ddelta_dtau) \
37 X(d3alphar_ddelta3) \
38 X(d3alphar_ddelta_dtau2) \
39 X(d3alphar_ddelta2_dtau) \
40 X(d3alphar_dtau3) \
41 X(d4alphar_ddelta4) \
42 X(d4alphar_ddelta3_dtau) \
43 X(d4alphar_ddelta2_dtau2) \
44 X(d4alphar_ddelta_dtau3) \
45 X(d4alphar_dtau4) \
46 X(delta_x_dalphar_ddelta) \
47 X(tau_x_dalphar_dtau) \
48 X(delta2_x_d2alphar_ddelta2) \
49 X(deltatau_x_d2alphar_ddelta_dtau) \
50 X(tau2_x_d2alphar_dtau2)
51
52struct HelmholtzDerivatives
53{
54#define X(name) CoolPropDbl name;
55 LIST_OF_DERIVATIVE_VARIABLES
56#undef X
57 CoolPropDbl tau, delta, T_red, rhomolar_red;
58
59 void reset(CoolPropDbl v){
60#define X(name) name = v;
61 LIST_OF_DERIVATIVE_VARIABLES
62#undef X
63 } HelmholtzDerivatives operator+(const HelmholtzDerivatives& other) const {
64 HelmholtzDerivatives _new;
65#define X(name) _new.name = name + other.name;
66 LIST_OF_DERIVATIVE_VARIABLES
67#undef X
68 return _new;
69 }
70 HelmholtzDerivatives operator*(const CoolPropDbl& other) const {
71 HelmholtzDerivatives _new;
72#define X(name) _new.name = name * other;
73 LIST_OF_DERIVATIVE_VARIABLES
74#undef X
75 return _new;
76 }
77 HelmholtzDerivatives() : tau(0.0), delta(0.0), T_red(_HUGE), rhomolar_red(_HUGE) {
78 reset(0.0);
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 };
474 void add_DoubleExponential(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& d, const std::vector<CoolPropDbl>& t,
475 const std::vector<CoolPropDbl>& gd, const std::vector<CoolPropDbl>& ld, const std::vector<CoolPropDbl>& gt,
476 const std::vector<CoolPropDbl>& lt) {
477 for (std::size_t i = 0; i < n.size(); ++i) {
478 ResidualHelmholtzGeneralizedExponentialElement el;
479 el.n = n[i];
480 el.d = d[i];
481 el.t = t[i];
482 el.c = gd[i];
483 el.l_double = ld[i];
484 el.l_int = (int)el.l_double;
485 el.omega = gt[i];
486 el.m_double = lt[i];
487 elements.push_back(el);
488 }
489 delta_li_in_u = true;
490 tau_mi_in_u = true;
491 };
492
493 void finish() {
494 n.resize(elements.size());
495 d.resize(elements.size());
496 t.resize(elements.size());
497 c.resize(elements.size());
498 omega.resize(elements.size());
499 l_double.resize(elements.size());
500 l_int.resize(elements.size());
501 m_double.resize(elements.size());
502 m_int.resize(elements.size());
503 epsilon2.resize(elements.size());
504 eta2.resize(elements.size());
505 gamma2.resize(elements.size());
506 beta2.resize(elements.size());
507
508 for (std::size_t i = 0; i < elements.size(); ++i) {
509 n[i] = elements[i].n;
510 d[i] = elements[i].d;
511 t[i] = elements[i].t;
512 c[i] = elements[i].c;
513 omega[i] = elements[i].omega;
514 l_double[i] = elements[i].l_double;
515 l_int[i] = elements[i].l_int;
516 m_double[i] = elements[i].m_double;
517 m_int[i] = elements[i].m_int;
518 epsilon2[i] = elements[i].epsilon2;
519 eta2[i] = elements[i].eta2;
520 gamma2[i] = elements[i].gamma2;
521 beta2[i] = elements[i].beta2;
522
523 // See if l is an integer, and store a flag if it is
524 elements[i].l_is_int = (std::abs(static_cast<long>(elements[i].l_double) - elements[i].l_double) < 1e-14);
525 }
526 // uE.resize(elements.size());
527 // du_ddeltaE.resize(elements.size());
528 // du_dtauE.resize(elements.size());
529 // d2u_ddelta2E.resize(elements.size());
530 // d2u_dtau2E.resize(elements.size());
531 // d3u_ddelta3E.resize(elements.size());
532 // d3u_dtau3E.resize(elements.size());
533
534 finished = true;
535 };
536
537 void to_json(rapidjson::Value& el, rapidjson::Document& doc);
538
539 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
540 //void allEigen(const CoolPropDbl &tau, const CoolPropDbl &delta, HelmholtzDerivatives &derivs) throw();
541
542#if ENABLE_CATCH
543 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
544#endif
545};
546
547struct ResidualHelmholtzNonAnalyticElement
548{
549 CoolPropDbl n, a, b, beta, A, B, C, D;
550};
551class ResidualHelmholtzNonAnalytic : public BaseHelmholtzTerm
552{
553
554 public:
555 std::size_t N;
556 std::vector<CoolPropDbl> s;
557 std::vector<ResidualHelmholtzNonAnalyticElement> elements;
559 ResidualHelmholtzNonAnalytic() {
560 N = 0;
561 };
563 ~ResidualHelmholtzNonAnalytic() {};
565 ResidualHelmholtzNonAnalytic(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& a, const std::vector<CoolPropDbl>& b,
566 const std::vector<CoolPropDbl>& beta, const std::vector<CoolPropDbl>& A, const std::vector<CoolPropDbl>& B,
567 const std::vector<CoolPropDbl>& C, const std::vector<CoolPropDbl>& D) {
568 N = n.size();
569 s.resize(N);
570 for (std::size_t i = 0; i < n.size(); ++i) {
571 ResidualHelmholtzNonAnalyticElement el;
572 el.n = n[i];
573 el.a = a[i];
574 el.b = b[i];
575 el.beta = beta[i];
576 el.A = A[i];
577 el.B = B[i];
578 el.C = C[i];
579 el.D = D[i];
580 elements.push_back(el);
581 }
582 };
583 void to_json(rapidjson::Value& el, rapidjson::Document& doc);
584 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
585#if ENABLE_CATCH
586 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
587#endif
588};
589
590class ResidualHelmholtzGeneralizedCubic : public BaseHelmholtzTerm
591{
592 protected:
593 shared_ptr<AbstractCubic> m_abstractcubic;
594 std::vector<double> z;
595 public:
596 bool enabled;
597
599 ResidualHelmholtzGeneralizedCubic() {
600 enabled = false;
601 };
603 ResidualHelmholtzGeneralizedCubic(shared_ptr<AbstractCubic>& ac) : m_abstractcubic(ac) {
604 enabled = true;
605 z = std::vector<double>(1, 1); // Init the vector to [1.0]
606 };
607
608 void to_json(rapidjson::Value& el, rapidjson::Document& doc);
609 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
610};
611
612class ResidualHelmholtzGaoB : public BaseHelmholtzTerm
613{
614 protected:
615 std::vector<double> n, t, d, eta, beta, gamma, epsilon, b;
616
617 public:
618 bool enabled;
619
621 ResidualHelmholtzGaoB() {
622 enabled = false;
623 };
624
626 ResidualHelmholtzGaoB(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& t, const std::vector<CoolPropDbl>& d,
627 const std::vector<CoolPropDbl>& eta, const std::vector<CoolPropDbl>& beta, const std::vector<CoolPropDbl>& gamma,
628 const std::vector<CoolPropDbl>& epsilon, const std::vector<CoolPropDbl>& b)
629 : n(n), t(t), d(d), eta(eta), beta(beta), gamma(gamma), epsilon(epsilon), b(b) {
630 enabled = true;
631 };
632
633 void to_json(rapidjson::Value& el, rapidjson::Document& doc);
634 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
635
636#if ENABLE_CATCH
637 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
638#endif
639};
640
642class ResidualHelmholtzXiangDeiters : public BaseHelmholtzTerm
643{
644
645 public:
646 bool enabled;
647 ResidualHelmholtzGeneralizedExponential phi0, phi1, phi2;
648 CoolPropDbl Tc, pc, rhomolarc, acentric, R, theta;
650 ResidualHelmholtzXiangDeiters() : Tc(_HUGE), pc(_HUGE), rhomolarc(_HUGE), acentric(_HUGE), R(_HUGE), theta(_HUGE) {
651 enabled = false;
652 };
654 ResidualHelmholtzXiangDeiters(const CoolPropDbl Tc, const CoolPropDbl pc, const CoolPropDbl rhomolarc, const CoolPropDbl acentric,
655 const CoolPropDbl R);
656 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
657#if ENABLE_CATCH
658 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
659#endif
660};
661
662class ResidualHelmholtzSAFTAssociating : public BaseHelmholtzTerm
663{
664
665 protected:
666 double a, m, epsilonbar, vbarn, kappabar;
667
668 CoolPropDbl Deltabar(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
669 CoolPropDbl dDeltabar_ddelta__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
670 CoolPropDbl d2Deltabar_ddelta2__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
671 CoolPropDbl dDeltabar_dtau__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
672 CoolPropDbl d2Deltabar_dtau2__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
673 CoolPropDbl d2Deltabar_ddelta_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
674 CoolPropDbl d3Deltabar_dtau3__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
675 CoolPropDbl d3Deltabar_ddelta_dtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
676 CoolPropDbl d3Deltabar_ddelta3__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
677 CoolPropDbl d3Deltabar_ddelta2_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
678
679 CoolPropDbl X(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
680 CoolPropDbl dX_dDeltabar__constdelta(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
681 CoolPropDbl dX_ddelta__constDeltabar(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const;
682 CoolPropDbl dX_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
683 CoolPropDbl dX_ddelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
684 CoolPropDbl d2X_dtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
685 CoolPropDbl d2X_ddeltadtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
686 CoolPropDbl d2X_ddelta2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
687
688 CoolPropDbl d3X_dtau3(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
689 CoolPropDbl d3X_ddelta3(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
690 CoolPropDbl d3X_ddeltadtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
691 CoolPropDbl d3X_ddelta2dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const;
692
693 CoolPropDbl g(const CoolPropDbl& eta) const;
694 CoolPropDbl dg_deta(const CoolPropDbl& eta) const;
695 CoolPropDbl d2g_deta2(const CoolPropDbl& eta) const;
696 CoolPropDbl d3g_deta3(const CoolPropDbl& eta) const;
697 CoolPropDbl eta(const CoolPropDbl& delta) const;
698
699 public:
701 ResidualHelmholtzSAFTAssociating() : a(_HUGE), m(_HUGE), epsilonbar(_HUGE), vbarn(_HUGE), kappabar(_HUGE) {
702 disabled = true;
703 };
704
705 // Constructor
706 ResidualHelmholtzSAFTAssociating(double a, double m, double epsilonbar, double vbarn, double kappabar)
707 : a(a), m(m), epsilonbar(epsilonbar), vbarn(vbarn), kappabar(kappabar) {
708 disabled = false;
709 };
710
711 bool disabled;
712
713 //Destructor. No Implementation
714 ~ResidualHelmholtzSAFTAssociating() {};
715
716 void to_json(rapidjson::Value& el, rapidjson::Document& doc);
717
718 CoolPropDbl dTau4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
719 return 1e99;
720 };
721 CoolPropDbl dDelta_dTau3(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
722 return 1e99;
723 };
724 CoolPropDbl dDelta2_dTau2(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
725 return 1e99;
726 };
727 CoolPropDbl dDelta3_dTau(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
728 return 1e99;
729 };
730 CoolPropDbl dDelta4(const CoolPropDbl& tau, const CoolPropDbl& delta) throw() {
731 return 1e99;
732 };
733
734 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& deriv) throw();
735};
736
737class BaseHelmholtzContainer
738{
739 protected:
740 std::array<double, 16> cache = create_filled_array<double, 16>(_HUGE);
741 std::array<bool, 16> is_cached = create_filled_array<bool, 16>(false);
742 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,
743 i30 = 12, i31 = 13, i40 = 14;
744
745 bool cache_valid(std::size_t i) const {
746 return is_cached[i];
747 }
748
749 public:
750 void clear() {
751 memset(cache.data(), 0, sizeof(cache));
752 memset(is_cached.data(), false, sizeof(is_cached));
753 };
754
755 virtual void empty_the_EOS() = 0;
756 virtual HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values) = 0;
757
758 CoolPropDbl base(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
759 if (!cache_valid(i00) || dont_use_cache)
760 return all(tau, delta, false).alphar;
761 else
762 return cache[i00];
763 };
764 CoolPropDbl dDelta(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
765 if (!cache_valid(i10) || dont_use_cache)
766 return all(tau, delta, false).dalphar_ddelta;
767 else
768 return cache[i10];
769 };
770 CoolPropDbl dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
771 if (!cache_valid(i01) || dont_use_cache)
772 return all(tau, delta, false).dalphar_dtau;
773 else
774 return cache[i01];
775 };
776 CoolPropDbl dDelta2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
777 if (!cache_valid(i20) || dont_use_cache)
778 return all(tau, delta, false).d2alphar_ddelta2;
779 else
780 return cache[i20];
781 };
782 CoolPropDbl dDelta_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
783 if (!cache_valid(i11) || dont_use_cache)
784 return all(tau, delta, false).d2alphar_ddelta_dtau;
785 else
786 return cache[i11];
787 };
788 CoolPropDbl dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
789 if (!cache_valid(i02) || dont_use_cache)
790 return all(tau, delta, false).d2alphar_dtau2;
791 else
792 return cache[i02];
793 };
794 CoolPropDbl dDelta3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
795 if (!cache_valid(i30) || dont_use_cache)
796 return all(tau, delta, false).d3alphar_ddelta3;
797 else
798 return cache[i30];
799 };
800 CoolPropDbl dDelta2_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
801 if (!cache_valid(i21) || dont_use_cache)
802 return all(tau, delta, false).d3alphar_ddelta2_dtau;
803 else
804 return cache[i21];
805 };
806 CoolPropDbl dDelta_dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
807 if (!cache_valid(i12) || dont_use_cache)
808 return all(tau, delta, false).d3alphar_ddelta_dtau2;
809 else
810 return cache[i12];
811 };
812 CoolPropDbl dTau3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
813 if (!cache_valid(i03) || dont_use_cache)
814 return all(tau, delta, false).d3alphar_dtau3;
815 else
816 return cache[i03];
817 };
818 CoolPropDbl dDelta4(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
819 return all(tau, delta, false).d4alphar_ddelta4;
820 };
821 CoolPropDbl dDelta3_dTau(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
822 return all(tau, delta, false).d4alphar_ddelta3_dtau;
823 };
824 CoolPropDbl dDelta2_dTau2(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
825 return all(tau, delta, false).d4alphar_ddelta2_dtau2;
826 };
827 CoolPropDbl dDelta_dTau3(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
828 return all(tau, delta, false).d4alphar_ddelta_dtau3;
829 };
830 CoolPropDbl dTau4(CoolPropDbl tau, CoolPropDbl delta, const bool dont_use_cache = false) {
831 return all(tau, delta, false).d4alphar_dtau4;
832 };
833};
834
835class ResidualHelmholtzContainer : public BaseHelmholtzContainer
836{
837 public:
838 ResidualHelmholtzNonAnalytic NonAnalytic;
839 ResidualHelmholtzSAFTAssociating SAFT;
840 ResidualHelmholtzGeneralizedExponential GenExp;
841 ResidualHelmholtzGeneralizedCubic cubic;
842 ResidualHelmholtzXiangDeiters XiangDeiters;
843 ResidualHelmholtzGaoB GaoB;
844
845 void empty_the_EOS() {
846 NonAnalytic = ResidualHelmholtzNonAnalytic();
847 SAFT = ResidualHelmholtzSAFTAssociating();
848 GenExp = ResidualHelmholtzGeneralizedExponential();
849 cubic = ResidualHelmholtzGeneralizedCubic();
850 XiangDeiters = ResidualHelmholtzXiangDeiters();
851 GaoB = ResidualHelmholtzGaoB();
852 };
853
854 HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values = false) {
855 HelmholtzDerivatives derivs; // zeros out the elements
856 GenExp.all(tau, delta, derivs);
857 NonAnalytic.all(tau, delta, derivs);
858 SAFT.all(tau, delta, derivs);
859 cubic.all(tau, delta, derivs);
860 XiangDeiters.all(tau, delta, derivs);
861 GaoB.all(tau, delta, derivs);
862 if (cache_values) {
863 cache[i00] = derivs.alphar;
864 cache[i10] = derivs.dalphar_ddelta;
865 cache[i01] = derivs.dalphar_dtau;
866 cache[i20] = derivs.d2alphar_ddelta2;
867 cache[i02] = derivs.d2alphar_dtau2;
868 cache[i11] = derivs.d2alphar_ddelta_dtau;
869 cache[i30] = derivs.d3alphar_ddelta3;
870 cache[i03] = derivs.d3alphar_dtau3;
871 cache[i21] = derivs.d3alphar_ddelta2_dtau;
872 cache[i12] = derivs.d3alphar_ddelta_dtau2;
873 memset(is_cached.data(), true, sizeof(is_cached));
874 }
875 return derivs;
876 };
877};
878
879// #############################################################################
880// #############################################################################
881// #############################################################################
882// IDEAL GAS TERMS
883// #############################################################################
884// #############################################################################
885// #############################################################################
886
888
893class IdealHelmholtzLead : public BaseHelmholtzTerm
894{
895
896 private:
897 CoolPropDbl a1, a2;
898 bool enabled;
899
900 public:
901 // Default constructor
902 IdealHelmholtzLead() : a1(_HUGE), a2(_HUGE), enabled(false) {}
903
904 // Constructor
905 IdealHelmholtzLead(CoolPropDbl a1, CoolPropDbl a2) : a1(a1), a2(a2), enabled(true) {}
906
907 bool is_enabled() const {
908 return enabled;
909 }
910
911 void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
912 el.AddMember("type", "IdealHelmholtzLead", doc.GetAllocator());
913 el.AddMember("a1", static_cast<double>(a1), doc.GetAllocator());
914 el.AddMember("a2", static_cast<double>(a2), doc.GetAllocator());
915 };
916
917 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
918};
919
921
926class IdealHelmholtzEnthalpyEntropyOffset : public BaseHelmholtzTerm
927{
928 private:
929 CoolPropDbl a1, a2; // Use these variables internally
930 std::string reference;
931 bool enabled;
932
933 public:
934 IdealHelmholtzEnthalpyEntropyOffset() : a1(_HUGE), a2(_HUGE), enabled(false) {}
935
936 // Constructor
937 IdealHelmholtzEnthalpyEntropyOffset(CoolPropDbl a1, CoolPropDbl a2, const std::string& ref) : a1(a1), a2(a2), reference(ref), enabled(true) {}
938
939 // Set the values in the class
940 void set(CoolPropDbl a1, CoolPropDbl a2, const std::string& ref) {
941 // If it doesn't already exist, just set the values
942 if (enabled == false) {
943 this->a1 = a1;
944 this->a2 = a2;
945 enabled = true;
946 } else if (ref == "DEF") {
947 this->a1 = 0.0;
948 this->a2 = 0.0;
949 enabled = false;
950 } else {
951 // Otherwise, increment the values
952 this->a1 += a1;
953 this->a2 += a2;
954 enabled = true;
955 }
956 this->reference = ref;
957 }
958
959 bool is_enabled() const {
960 return enabled;
961 };
962
963 void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
964 el.AddMember("type", "IdealHelmholtzEnthalpyEntropyOffset", doc.GetAllocator());
965 el.AddMember("a1", static_cast<double>(a1), doc.GetAllocator());
966 el.AddMember("a2", static_cast<double>(a2), doc.GetAllocator());
967 };
968 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
969};
970
976class IdealHelmholtzLogTau : public BaseHelmholtzTerm
977{
978 private:
979 CoolPropDbl a1;
980 bool enabled;
981
982 public:
984 IdealHelmholtzLogTau() : a1(_HUGE), enabled(false) {}
985
986 // Constructor
987 IdealHelmholtzLogTau(CoolPropDbl a1) : a1(a1), enabled(true) {}
988
989 bool is_enabled() const {
990 return enabled;
991 };
992
993 void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
994 el.AddMember("type", "IdealHelmholtzLogTau", doc.GetAllocator());
995 el.AddMember("a1", static_cast<double>(a1), doc.GetAllocator());
996 };
997 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
998};
999
1005class IdealHelmholtzPower : public BaseHelmholtzTerm
1006{
1007
1008 private:
1009 std::vector<CoolPropDbl> n, t; // Use these variables internally
1010 std::size_t N;
1011 bool enabled;
1012
1013 public:
1014 IdealHelmholtzPower() : N(0), enabled(false) {};
1015 // Constructor
1016 IdealHelmholtzPower(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& t) : n(n), t(t), N(n.size()), enabled(true) {};
1017
1018 bool is_enabled() const {
1019 return enabled;
1020 };
1021
1022 void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
1023 el.AddMember("type", "IdealHelmholtzPower", doc.GetAllocator());
1024 cpjson::set_long_double_array("n", n, el, doc);
1025 cpjson::set_long_double_array("t", t, el, doc);
1026 };
1027 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1028};
1029
1076class IdealHelmholtzPlanckEinsteinGeneralized : public BaseHelmholtzTerm
1077{
1078
1079 private:
1080 std::vector<CoolPropDbl> n, theta, c, d; // Use these variables internally
1081 std::size_t N;
1082 bool enabled;
1083
1084 public:
1085 IdealHelmholtzPlanckEinsteinGeneralized() : N(0), enabled(false) {}
1086 // Constructor with std::vector instances
1087 IdealHelmholtzPlanckEinsteinGeneralized(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta,
1088 const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& d)
1089 : n(n), theta(theta), c(c), d(d), N(n.size()), enabled(true) {}
1090
1091 // Extend the vectors to allow for multiple instances feeding values to this function
1092 void extend(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, const std::vector<CoolPropDbl>& c,
1093 const std::vector<CoolPropDbl>& d) {
1094 this->n.insert(this->n.end(), n.begin(), n.end());
1095 this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1096 this->c.insert(this->c.end(), c.begin(), c.end());
1097 this->d.insert(this->d.end(), d.begin(), d.end());
1098 N += n.size();
1099 }
1100
1101 bool is_enabled() const {
1102 return enabled;
1103 };
1104
1105 void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
1106 el.AddMember("type", "IdealHelmholtzPlanckEinsteinGeneralized", doc.GetAllocator());
1107 cpjson::set_long_double_array("n", n, el, doc);
1108 cpjson::set_long_double_array("theta", theta, el, doc);
1109 };
1110 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1111};
1112
1113class IdealHelmholtzCP0Constant : public BaseHelmholtzTerm
1114{
1115
1116 private:
1117 double cp_over_R, Tc, T0, tau0; // Use these variables internally
1118 bool enabled;
1119
1120 public:
1122 IdealHelmholtzCP0Constant() : cp_over_R(_HUGE), Tc(_HUGE), T0(_HUGE), tau0(_HUGE) {
1123 enabled = false;
1124 };
1125
1127 IdealHelmholtzCP0Constant(CoolPropDbl cp_over_R, CoolPropDbl Tc, CoolPropDbl T0) : cp_over_R(cp_over_R), Tc(Tc), T0(T0) {
1128 enabled = true;
1129 tau0 = Tc / T0;
1130 };
1131
1133 ~IdealHelmholtzCP0Constant() {};
1134
1135 bool is_enabled() const {
1136 return enabled;
1137 };
1138
1139 void to_json(rapidjson::Value& el, rapidjson::Document& doc) {
1140 el.AddMember("type", "IdealGasHelmholtzCP0Constant", doc.GetAllocator());
1141 el.AddMember("cp_over_R", cp_over_R, doc.GetAllocator());
1142 el.AddMember("Tc", Tc, doc.GetAllocator());
1143 el.AddMember("T0", T0, doc.GetAllocator());
1144 };
1145
1146 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw();
1147};
1148
1149class IdealHelmholtzCP0PolyT : public BaseHelmholtzTerm
1150{
1151 private:
1152 std::vector<CoolPropDbl> c, t;
1153 CoolPropDbl Tc, T0, tau0; // Use these variables internally
1154 std::size_t N;
1155 bool enabled;
1156
1157 public:
1158 IdealHelmholtzCP0PolyT() : Tc(_HUGE), T0(_HUGE), tau0(_HUGE), N(0), enabled(false) {}
1159
1161 IdealHelmholtzCP0PolyT(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t, double Tc, double T0)
1162 : c(c), t(t), Tc(Tc), T0(T0), tau0(Tc / T0), N(c.size()), enabled(true) {
1163 assert(c.size() == t.size());
1164 }
1165
1166 void extend(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t) {
1167 this->c.insert(this->c.end(), c.begin(), c.end());
1168 this->t.insert(this->t.end(), t.begin(), t.end());
1169 N += c.size();
1170 }
1171
1172 bool is_enabled() const {
1173 return enabled;
1174 };
1175
1176 void to_json(rapidjson::Value& el, rapidjson::Document& doc);
1177 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
1178#if ENABLE_CATCH
1179 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
1180#endif
1181};
1185class IdealHelmholtzGERG2004Sinh : public BaseHelmholtzTerm
1186{
1187 private:
1188 std::vector<CoolPropDbl> n, theta;
1189 CoolPropDbl Tc, _Tr;
1190 std::size_t N;
1191 bool enabled;
1192
1193 public:
1194 IdealHelmholtzGERG2004Sinh() : Tc(_HUGE), _Tr(_HUGE), N(0), enabled(false) {}
1195
1197 IdealHelmholtzGERG2004Sinh(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, double Tc)
1198 : n(n), theta(theta), Tc(Tc), _Tr(_HUGE), N(n.size()), enabled(true) {
1199 assert(n.size() == theta.size());
1200 }
1201
1202 void extend(const std::vector<CoolPropDbl>& c, const std::vector<CoolPropDbl>& t) {
1203 this->n.insert(this->n.end(), n.begin(), n.end());
1204 this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1205 N += c.size();
1206 }
1207 void set_Tred(CoolPropDbl Tr) {
1208 this->_Tr = Tr;
1209 }
1210
1211 bool is_enabled() const {
1212 return enabled;
1213 };
1214 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
1215
1216#if ENABLE_CATCH
1217 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
1218#endif
1219};
1220
1221class IdealHelmholtzGERG2004Cosh : public BaseHelmholtzTerm
1222{
1223 private:
1224 std::vector<CoolPropDbl> n, theta;
1225 CoolPropDbl Tc, _Tr;
1226 std::size_t N;
1227 bool enabled;
1228
1229 public:
1230 IdealHelmholtzGERG2004Cosh() : Tc(_HUGE), _Tr(_HUGE), N(0), enabled(false) {}
1231
1233 IdealHelmholtzGERG2004Cosh(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta, double Tc)
1234 : n(n), theta(theta), Tc(Tc), _Tr(_HUGE), N(n.size()), enabled(true) {
1235 assert(n.size() == theta.size());
1236 }
1237
1238 void extend(const std::vector<CoolPropDbl>& n, const std::vector<CoolPropDbl>& theta) {
1239 this->n.insert(this->n.end(), n.begin(), n.end());
1240 this->theta.insert(this->theta.end(), theta.begin(), theta.end());
1241 N += n.size();
1242 }
1243 void set_Tred(CoolPropDbl Tr) {
1244 this->_Tr = Tr;
1245 }
1246
1247 bool is_enabled() const {
1248 return enabled;
1249 };
1250 void all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) override;
1251
1252#if ENABLE_CATCH
1253 virtual mcx::MultiComplex<double> one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const override;
1254#endif
1255};
1256
1259//\f[
1260//\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
1261//\f]
1262//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
1263//\f[
1264//\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
1265//\f]
1266//or in terms of \f$ \tau \f$:
1267//\f[
1268//\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
1269//\f]
1270//Third partial:
1271//\f[
1272//\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)}
1273//\f]
1274//Now coming back to the ideal gas Helmholtz energy definition:
1275//\f[
1276//\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
1277//\f]
1278//Applying derivative
1279//\f[
1280//\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]
1281//\f]
1282//Fundamental theorem of calculus
1283//\f[
1284//\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}
1285//\f]
1286//Last two terms cancel, leaving
1287//\f[
1288//\alpha^0_{\tau} = -\int_{\tau_0}^{\tau} \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau
1289//\f]
1290//Another derivative yields (from fundamental theorem of calculus)
1291//\f[
1292//\alpha^0_{\tau\tau} = - \frac{1}{\tau^2}\frac{c_p^0}{R_u}
1293//\f]
1294//
1295//see also Jaeschke and Schley, 1995, (http://link.springer.com/article/10.1007%2FBF02083547#page-1)
1296//*/
1298//class IdealHelmholtzCP0AlyLee : public BaseHelmholtzTerm{
1299//private:
1300// std::vector<CoolPropDbl> c;
1301// CoolPropDbl Tc, tau0, T0; // Use these variables internally
1302// bool enabled;
1303//public:
1304// IdealHelmholtzCP0AlyLee(){enabled = false;};
1305//
1306// /// Constructor with std::vectors
1307// IdealHelmholtzCP0AlyLee(const std::vector<CoolPropDbl> &c, double Tc, double T0)
1308// :c(c), Tc(Tc), T0(T0)
1309// {
1310// tau0=Tc/T0;
1311// enabled = true;
1312// };
1313//
1314// /// Destructor
1315// ~IdealHelmholtzCP0AlyLee(){};
1316//
1317// bool is_enabled() const {return enabled;};
1318//
1319// void to_json(rapidjson::Value &el, rapidjson::Document &doc);
1320//
1321//
1322// /// The antiderivative given by \f$ \displaystyle\int \frac{1}{\tau^2}\frac{c_p^0}{R_u}d\tau \f$
1323// /**
1324// sympy code for this derivative:
1325//
1326// from sympy import *
1327// a1,a2,a3,a4,a5,Tc,tau = symbols('a1,a2,a3,a4,a5,Tc,tau', real = True)
1328// integrand = a1 + a2*(a3/Tc/sinh(a3*tau/Tc))**2 + a4*(a5/Tc/cosh(a5*tau/Tc))**2
1329// integrand = integrand.rewrite(exp)
1330// antideriv = trigsimp(integrate(integrand,tau))
1331// display(antideriv)
1332// print latex(antideriv)
1333// print ccode(antideriv)
1334//
1335// \f[
1336// \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]}
1337// \f]
1338// */
1339// CoolPropDbl anti_deriv_cp0_tau2(const CoolPropDbl &tau);
1340//
1341// /// The antiderivative given by \f$ \displaystyle\int \frac{1}{\tau}\frac{c_p^0}{R_u}d\tau \f$
1342// /**
1343// sympy code for this derivative:
1344//
1345// 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)
1346// 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
1347//
1348// term2 = a_1/tau*(a_2*tau/Tc/sinh(a_2*tau/Tc))**2
1349// term2 = term2.rewrite(exp) # Unpack the sinh to exp functions
1350// antideriv2 = trigsimp(integrate(term2,tau))
1351// display(antideriv2)
1352// print latex(antideriv2)
1353// print ccode(antideriv2)
1354//
1355// term3 = a_3/tau*(a_4*tau/Tc/cosh(a_4*tau/Tc))**2
1356// term3 = term3.rewrite(exp) # Unpack the cosh to exp functions
1357// antideriv3 = factor(trigsimp(integrate(term3,tau).rewrite(exp)))
1358// display(antideriv3)
1359// print latex(antideriv3)
1360// print ccode(antideriv3)
1361//
1362// Can be broken into three parts (trick is to express \f$sinh\f$ and \f$cosh\f$ in terms of \f$exp\f$ function)
1363//
1364// Term 2:
1365// \f[
1366// \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
1367// \f]
1368//
1369// Term 3:
1370// \f[
1371// \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)
1372// \f]
1373// */
1374// CoolPropDbl anti_deriv_cp0_tau(const CoolPropDbl &tau);
1375//
1376// CoolPropDbl base(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1377// CoolPropDbl dDelta(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1378// CoolPropDbl dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1379// CoolPropDbl dDelta2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1380// CoolPropDbl dDelta_dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1381// CoolPropDbl dTau2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1382// CoolPropDbl dDelta3(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1383// CoolPropDbl dDelta2_dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1384// CoolPropDbl dDelta_dTau2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw(){return 0.0;};
1385// CoolPropDbl dTau3(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1386// CoolPropDbl dTau4(const CoolPropDbl &tau, const CoolPropDbl &delta) throw();
1387//
1388//};
1389
1390class IdealHelmholtzContainer : public BaseHelmholtzContainer
1391{
1392 private:
1393 double _prefactor;
1394
1395 public:
1396 IdealHelmholtzLead Lead;
1397 IdealHelmholtzEnthalpyEntropyOffset EnthalpyEntropyOffsetCore, EnthalpyEntropyOffset;
1398 IdealHelmholtzLogTau LogTau;
1399 IdealHelmholtzPower Power;
1400 IdealHelmholtzPlanckEinsteinGeneralized PlanckEinstein;
1401
1402 IdealHelmholtzCP0Constant CP0Constant;
1403 IdealHelmholtzCP0PolyT CP0PolyT;
1404 IdealHelmholtzGERG2004Cosh GERG2004Cosh;
1405 IdealHelmholtzGERG2004Sinh GERG2004Sinh;
1406
1407 IdealHelmholtzContainer() : _prefactor(1.0) {};
1408
1409 void set_prefactor(double prefactor) {
1410 _prefactor = prefactor;
1411 }
1412
1413 void set_Tred(double T_red) {
1414 GERG2004Cosh.set_Tred(T_red);
1415 GERG2004Sinh.set_Tred(T_red);
1416 }
1417
1418 void empty_the_EOS() {
1419 Lead = IdealHelmholtzLead();
1420 EnthalpyEntropyOffsetCore = IdealHelmholtzEnthalpyEntropyOffset();
1421 EnthalpyEntropyOffset = IdealHelmholtzEnthalpyEntropyOffset();
1422 LogTau = IdealHelmholtzLogTau();
1423 Power = IdealHelmholtzPower();
1424 PlanckEinstein = IdealHelmholtzPlanckEinsteinGeneralized();
1425 CP0Constant = IdealHelmholtzCP0Constant();
1426 CP0PolyT = IdealHelmholtzCP0PolyT();
1427 GERG2004Cosh = IdealHelmholtzGERG2004Cosh();
1428 GERG2004Sinh = IdealHelmholtzGERG2004Sinh();
1429 };
1430
1431 HelmholtzDerivatives all(const CoolPropDbl tau, const CoolPropDbl delta, bool cache_values = false) {
1432 HelmholtzDerivatives derivs; // zeros out the elements
1433 Lead.all(tau, delta, derivs);
1434 EnthalpyEntropyOffsetCore.all(tau, delta, derivs);
1435 EnthalpyEntropyOffset.all(tau, delta, derivs);
1436 LogTau.all(tau, delta, derivs);
1437 Power.all(tau, delta, derivs);
1438 PlanckEinstein.all(tau, delta, derivs);
1439 CP0Constant.all(tau, delta, derivs);
1440 CP0PolyT.all(tau, delta, derivs);
1441 GERG2004Cosh.all(tau, delta, derivs);
1442 GERG2004Sinh.all(tau, delta, derivs);
1443
1444 if (cache_values) {
1445 cache[i00] = derivs.alphar * _prefactor;
1446 cache[i10] = derivs.dalphar_ddelta * _prefactor;
1447 cache[i01] = derivs.dalphar_dtau * _prefactor;
1448 cache[i20] = derivs.d2alphar_ddelta2 * _prefactor;
1449 cache[i02] = derivs.d2alphar_dtau2 * _prefactor;
1450 cache[i11] = derivs.d2alphar_ddelta_dtau * _prefactor;
1451 cache[i30] = derivs.d3alphar_ddelta3 * _prefactor;
1452 cache[i03] = derivs.d3alphar_dtau3 * _prefactor;
1453 cache[i21] = derivs.d3alphar_ddelta2_dtau * _prefactor;
1454 cache[i12] = derivs.d3alphar_ddelta_dtau2 * _prefactor;
1455 memset(is_cached.data(), true, sizeof(is_cached));
1456 }
1457 return derivs * _prefactor;
1458 };
1459};
1460}; /* namespace CoolProp */
1461
1462#endif
Generated on Sat Nov 1 2025 19:25:21 for CoolProp by doxygen 1.9.4