CoolProp 8.0.0
An open-source fluid property and humid air property database
GeneralizedCubic.h
Go to the documentation of this file.
1
11#ifndef CUBIC_H
12#define CUBIC_H
13
14#include <utility>
15#include <vector>
16#include <cmath>
17#include <cstring>
18#include <limits>
19#include <array>
20#include <memory>
21using std::shared_ptr;
22#include "CoolProp/Exceptions.h"
23
25class AbstractCubicAlphaFunction
26{
27 protected:
28 double a0;
29 double Tr_over_Tci,
30 sqrt_Tr_Tci;
31 std::vector<double> c;
34 unsigned long m_version = 0;
35
36 public:
37 virtual ~AbstractCubicAlphaFunction() = default;
47 virtual double term(double tau, std::size_t itau) = 0;
49 virtual void calc_all_terms(double tau, std::array<double, 5>& terms) {
50 for (int k = 0; k < 5; ++k)
51 terms[k] = term(tau, k);
52 }
53 void set_Tr_over_Tci(double Tr_over_Tci) {
54 this->Tr_over_Tci = Tr_over_Tci;
55 this->sqrt_Tr_Tci = sqrt(Tr_over_Tci);
56 ++m_version;
57 };
61 unsigned long version() const {
62 return m_version;
63 }
64 AbstractCubicAlphaFunction(double a0, double Tr_over_Tci) : a0(a0), Tr_over_Tci(Tr_over_Tci), sqrt_Tr_Tci(sqrt(Tr_over_Tci)) {};
65};
66
68class BasicMathiasCopemanAlphaFunction : public AbstractCubicAlphaFunction
69{
70 double m;
71 public:
72 BasicMathiasCopemanAlphaFunction(double a0, double m_ii, double Tr_over_Tci)
73 : AbstractCubicAlphaFunction(a0, Tr_over_Tci), m(m_ii) {
74
75 };
76 double term(double tau, std::size_t itau) override;
77 void calc_all_terms(double tau, std::array<double, 5>& terms) override;
78};
79
81class TwuAlphaFunction : public AbstractCubicAlphaFunction
82{
83 public:
84 TwuAlphaFunction(double a0, double L, double M, double N, double Tr_over_Tci) : AbstractCubicAlphaFunction(a0, Tr_over_Tci) {
85 c.resize(3);
86 c[0] = L;
87 c[1] = M;
88 c[2] = N;
89 };
90 double term(double tau, std::size_t itau) override;
91 void calc_all_terms(double tau, std::array<double, 5>& terms) override;
92};
93
95class MathiasCopemanAlphaFunction : public AbstractCubicAlphaFunction
96{
97 public:
98 MathiasCopemanAlphaFunction(double a0, double c1, double c2, double c3, double Tr_over_Tci) : AbstractCubicAlphaFunction(a0, Tr_over_Tci) {
99 c.resize(3);
100 c[0] = c1;
101 c[1] = c2;
102 c[2] = c3;
103 };
104 double term(double tau, std::size_t itau) override;
105 void calc_all_terms(double tau, std::array<double, 5>& terms) override;
106};
107
108class AbstractCubic
109{
110
111 protected:
112 double rho_r,
113 T_r;
114 std::vector<double> Tc,
115 pc,
116 acentric;
117 double R_u;
118 double Delta_1,
119 Delta_2;
120 int N;
121 std::vector<std::vector<double>> k;
122 double cm;
123 std::vector<shared_ptr<AbstractCubicAlphaFunction>> alpha;
125 mutable double m_tau_cache;
126 mutable std::vector<std::array<double, 5>> m_aii_cache;
130 mutable std::vector<unsigned long> m_alpha_versions_cache;
134 mutable std::vector<double> m_b0_ii_cache;
135 void _ensure_aii_cache(double tau) const {
136 bool stale = (std::memcmp(&tau, &m_tau_cache, sizeof(double)) != 0) || (m_alpha_versions_cache.size() != alpha.size());
137 if (!stale) {
138 for (int i = 0; i < N; ++i) {
139 if (alpha[i]->version() != m_alpha_versions_cache[i]) {
140 stale = true;
141 break;
142 }
143 }
144 }
145 if (stale) {
146 m_aii_cache.resize(N);
147 m_alpha_versions_cache.resize(N);
148 for (int i = 0; i < N; ++i) {
149 alpha[i]->calc_all_terms(tau, m_aii_cache[i]);
150 m_alpha_versions_cache[i] = alpha[i]->version();
151 }
152 m_tau_cache = tau;
153 }
154 }
155
156 public:
166 AbstractCubic(const std::vector<double>& Tc, std::vector<double> pc, std::vector<double> acentric, double R_u, double Delta_1, double Delta_2,
167 const std::vector<double>& C1 = std::vector<double>(), const std::vector<double>& C2 = std::vector<double>(),
168 const std::vector<double>& C3 = std::vector<double>());
169 virtual ~AbstractCubic() = default;
171 void set_alpha(const std::vector<double>& C1, const std::vector<double>& C2, const std::vector<double>& C3);
173 void set_alpha_function(std::size_t i, shared_ptr<AbstractCubicAlphaFunction>& acaf) {
174 alpha[i] = acaf;
175 m_tau_cache = std::numeric_limits<double>::quiet_NaN(); // invalidate cache
176 };
178 shared_ptr<AbstractCubicAlphaFunction> get_alpha_function(std::size_t i) {
179 return alpha[i];
180 };
182 const std::vector<shared_ptr<AbstractCubicAlphaFunction>>& get_all_alpha_functions() {
183 return this->alpha;
184 };
186 void set_all_alpha_functions(const std::vector<shared_ptr<AbstractCubicAlphaFunction>>& alpha) {
187 this->alpha = alpha;
188 m_tau_cache = std::numeric_limits<double>::quiet_NaN(); // invalidate cache
189 };
190
192 const std::vector<std::vector<double>>& get_kmat() {
193 return k;
194 };
196 void set_kmat(const std::vector<std::vector<double>>& k) {
197 this->k = k;
198 };
200 void set_kij(std::size_t i, std::size_t j, double val) {
201 k[i][j] = val;
202 k[j][i] = val;
203 }
205 double get_kij(std::size_t i, std::size_t j) {
206 return k[i][j];
207 }
209 std::vector<double>& get_Tc() {
210 return Tc;
211 }
213 std::vector<double>& get_pc() {
214 return pc;
215 }
217 std::vector<double>& get_acentric() {
218 return acentric;
219 }
221 double get_Delta_1() {
222 return Delta_1;
223 }
225 double get_Delta_2() {
226 return Delta_2;
227 }
229 double get_R_u() {
230 return R_u;
231 }
233 void set_Tr(double Tr) {
234 T_r = Tr;
235 for (std::size_t i = 0; i < alpha.size(); ++i) {
236 alpha[i]->set_Tr_over_Tci(T_r / Tc[i]);
237 }
238 m_tau_cache = std::numeric_limits<double>::quiet_NaN(); // invalidate cache
239 }
241 void set_rhor(double rhor) {
242 rho_r = rhor;
243 }
245 double get_Tr() {
246 return T_r;
247 }
249 double get_rhor() {
250 return rho_r;
251 }
252
254 void set_C_MC(std::size_t i, double c1, double c2, double c3) {
255 alpha[i] = std::make_shared<MathiasCopemanAlphaFunction>(a0_ii(i), c1, c2, c3, T_r / Tc[i]);
256 m_tau_cache = std::numeric_limits<double>::quiet_NaN(); // invalidate cache
257 }
259 void set_C_Twu(std::size_t i, double L, double M, double N) {
260 alpha[i] = std::make_shared<TwuAlphaFunction>(a0_ii(i), L, M, N, T_r / Tc[i]);
261 m_tau_cache = std::numeric_limits<double>::quiet_NaN(); // invalidate cache
262 }
265 virtual double a0_ii(std::size_t i) = 0;
268 virtual double b0_ii(std::size_t i) = 0;
270 virtual double m_ii(std::size_t i) = 0;
271
273 virtual double alphar(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta);
275 virtual double d_alphar_dxi(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
276 bool xN_independent);
278 virtual double d2_alphar_dxidxj(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
279 std::size_t j, bool xN_independent);
281 virtual double d3_alphar_dxidxjdxk(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
282 std::size_t j, std::size_t k, bool xN_independent);
283
290 virtual double am_term(double tau, const std::vector<double>& x, std::size_t itau);
299 virtual double d_am_term_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent);
309 virtual double d2_am_term_dxidxj(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j, bool xN_independent);
320 virtual double d3_am_term_dxidxjdxk(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j, std::size_t k,
321 bool xN_independent);
322
327 virtual double bm_term(const std::vector<double>& x);
333 virtual double d_bm_term_dxi(const std::vector<double>& x, std::size_t i, bool xN_independent);
340 virtual double d2_bm_term_dxidxj(const std::vector<double>& x, std::size_t i, std::size_t j, bool xN_independent);
348 virtual double d3_bm_term_dxidxjdxk(const std::vector<double>& x, std::size_t i, std::size_t j, std::size_t k, bool xN_independent);
352 virtual double cm_term();
354 void set_cm(double val) {
355 cm = val;
356 }
358 double get_cm() {
359 return cm;
360 }
361
363 virtual void set_Q_k(const size_t sgi, const double value) {
364 throw CoolProp::ValueError("set_Q_k not defined for AbstractCubic");
365 };
367 virtual double get_Q_k(const size_t sgi) const {
368 throw CoolProp::ValueError("get_Q_k not defined for AbstractCubic");
369 };
370
378 double aij_term(double tau, std::size_t i, std::size_t j, std::size_t itau);
385 double u_term(double tau, std::size_t i, std::size_t j, std::size_t itau);
391 double aii_term(double tau, std::size_t i, std::size_t itau);
392
400 double psi_minus(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta);
410 double d_psi_minus_dxi(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i, bool xN_independent);
421 double d2_psi_minus_dxidxj(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i, std::size_t j,
422 bool xN_independent);
434 double d3_psi_minus_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i, std::size_t j,
435 std::size_t k, bool xN_independent);
436
445 double PI_12(double delta, const std::vector<double>& x, std::size_t idelta);
454 double d_PI_12_dxi(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, bool xN_independent);
464 double d2_PI_12_dxidxj(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j, bool xN_independent);
475 double d3_PI_12_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j, std::size_t k,
476 bool xN_independent);
477
484 double tau_times_a(double tau, const std::vector<double>& x, std::size_t itau);
493 double d_tau_times_a_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent);
503 double d2_tau_times_a_dxidxj(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j, bool xN_independent);
514 double d3_tau_times_a_dxidxjdxk(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j, std::size_t k,
515 bool xN_independent);
516
525 double psi_plus(double delta, const std::vector<double>& x, std::size_t idelta);
534 double d_psi_plus_dxi(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, bool xN_independent);
544 double d2_psi_plus_dxidxj(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j, bool xN_independent);
555 double d3_psi_plus_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j, std::size_t k,
556 bool xN_independent);
557
566 double c_term(const std::vector<double>& x) {
567 return 1 / bm_term(x);
568 };
575 double d_c_term_dxi(const std::vector<double>& x, std::size_t i, bool xN_independent) {
576 return -d_bm_term_dxi(x, i, xN_independent) / pow(bm_term(x), 2);
577 };
585 double d2_c_term_dxidxj(const std::vector<double>& x, std::size_t i, std::size_t j, bool xN_independent) {
586 double bm = bm_term(x);
587 return (2 * d_bm_term_dxi(x, i, xN_independent) * d_bm_term_dxi(x, j, xN_independent) - bm * d2_bm_term_dxidxj(x, i, j, xN_independent))
588 / pow(bm, 3);
589 };
598 double d3_c_term_dxidxjdxk(const std::vector<double>& x, std::size_t i, std::size_t j, std::size_t k, bool xN_independent) {
599 double bm = bm_term(x);
600 return 1 / pow(bm, 4)
601 * (2 * bm
602 * (d_bm_term_dxi(x, i, xN_independent) * d2_bm_term_dxidxj(x, j, k, xN_independent)
603 + d_bm_term_dxi(x, j, xN_independent) * d2_bm_term_dxidxj(x, i, k, xN_independent)
604 + d_bm_term_dxi(x, k, xN_independent) * d2_bm_term_dxidxj(x, i, j, xN_independent))
605 - pow(bm, 2) * d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent)
606 - 6 * d_bm_term_dxi(x, i, xN_independent) * d_bm_term_dxi(x, j, xN_independent) * d_bm_term_dxi(x, k, xN_independent));
607 };
608
619 double A_term(double delta, const std::vector<double>& x) {
620 double bm = bm_term(x);
621 double cm = cm_term();
622 return log((delta * rho_r * (Delta_1 * bm + cm) + 1) / (delta * rho_r * (Delta_2 * bm + cm) + 1));
623 };
631 double d_A_term_dxi(double delta, const std::vector<double>& x, std::size_t i, bool xN_independent) {
632 std::size_t idelta = 0;
633 return delta * rho_r * d_bm_term_dxi(x, i, xN_independent) * (Delta_1 - Delta_2) / PI_12(delta, x, idelta);
634 };
643 double d2_A_term_dxidxj(double delta, const std::vector<double>& x, std::size_t i, std::size_t j, bool xN_independent) {
644 std::size_t idelta = 0;
645 double PI12 = PI_12(delta, x, idelta);
646 return delta * rho_r * (Delta_1 - Delta_2) / pow(PI12, 2)
647 * (PI12 * d2_bm_term_dxidxj(x, i, j, xN_independent)
648 - d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_bm_term_dxi(x, i, xN_independent));
649 };
659 double d3_A_term_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t i, std::size_t j, std::size_t k, bool xN_independent) {
660 std::size_t idelta = 0;
661 double PI12 = PI_12(delta, x, idelta);
662 // The leading factor
663 double lead = delta * rho_r * (Delta_1 - Delta_2) / pow(PI12, 3);
664 return lead
665 * (-PI12
666 * (d_PI_12_dxi(delta, x, idelta, j, xN_independent) * d2_bm_term_dxidxj(x, i, k, xN_independent)
667 + d_PI_12_dxi(delta, x, idelta, k, xN_independent) * d2_bm_term_dxidxj(x, i, j, xN_independent)
668 + d_bm_term_dxi(x, i, xN_independent) * d2_PI_12_dxidxj(delta, x, idelta, j, k, xN_independent))
669 + pow(PI12, 2) * d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent)
670 + 2 * d_PI_12_dxi(delta, x, idelta, j, xN_independent) * d_PI_12_dxi(delta, x, idelta, k, xN_independent)
671 * d_bm_term_dxi(x, i, xN_independent));
672 };
673 // Allows to modify the unifac interaction parameters aij, bij and cij. Only for use with VTPR backend.
674 virtual void set_interaction_parameter(const std::size_t mgi1, const std::size_t mgi2, const std::string& parameter, const double value) {
675 throw CoolProp::NotImplementedError("set_interaction_parameter is not implemented for this backend");
676 }
677 // Allows to get the unifac interaction parameters aij, bij and cij. Only for use with VTPR backend.
678 virtual double get_interaction_parameter(const std::size_t mgi1, const std::size_t mgi2, const std::string& parameter) {
679 throw CoolProp::NotImplementedError("get_interaction_parameter is not implemented for this backend");
680 }
681};
682
683class PengRobinson : public AbstractCubic
684{
685 public:
686 PengRobinson(const std::vector<double>& Tc, std::vector<double> pc, std::vector<double> acentric, double R_u,
687 const std::vector<double>& C1 = std::vector<double>(), const std::vector<double>& C2 = std::vector<double>(),
688 const std::vector<double>& C3 = std::vector<double>())
689 : AbstractCubic(Tc, std::move(pc), std::move(acentric), R_u, 1 + sqrt(2.0), 1 - sqrt(2.0), C1, C2, C3) {
690 set_alpha(C1, C2, C3);
691 };
692
693 PengRobinson(double Tc, double pc, double acentric, double R_u)
694 : AbstractCubic(std::vector<double>(1, Tc), std::vector<double>(1, pc), std::vector<double>(1, acentric), R_u, 1 + sqrt(2.0), 1 - sqrt(2.0)) {
695 set_alpha(std::vector<double>(), std::vector<double>(), std::vector<double>());
696 };
697
698 double a0_ii(std::size_t i) override;
699 double b0_ii(std::size_t i) override;
700 double m_ii(std::size_t i) override;
701};
702
703class SRK : public AbstractCubic
704{
705 public:
706 SRK(const std::vector<double>& Tc, std::vector<double> pc, std::vector<double> acentric, double R_u,
707 const std::vector<double>& C1 = std::vector<double>(), const std::vector<double>& C2 = std::vector<double>(),
708 const std::vector<double>& C3 = std::vector<double>())
709 : AbstractCubic(Tc, std::move(pc), std::move(acentric), R_u, 1, 0, C1, C2, C3) {
710 set_alpha(C1, C2, C3);
711 };
712 SRK(double Tc, double pc, double acentric, double R_u)
713 : AbstractCubic(std::vector<double>(1, Tc), std::vector<double>(1, pc), std::vector<double>(1, acentric), R_u, 1, 0) {
714 set_alpha(std::vector<double>(), std::vector<double>(), std::vector<double>());
715 };
716
717 double a0_ii(std::size_t i) override;
718 double b0_ii(std::size_t i) override;
719 double m_ii(std::size_t i) override;
720};
721
722#endif
Generated by doxygen 1.9.4