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