CoolProp 6.8.1dev
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 <vector>
15#include <cmath>
16#include "crossplatform_shared_ptr.h"
17#include "Exceptions.h"
18
20class AbstractCubicAlphaFunction
21{
22 protected:
23 double a0;
24 double Tr_over_Tci,
25 sqrt_Tr_Tci;
26 std::vector<double> c;
27 public:
28 virtual ~AbstractCubicAlphaFunction(){};
29 virtual double term(double tau, std::size_t itau) = 0;
30 void set_Tr_over_Tci(double Tr_over_Tci) {
31 this->Tr_over_Tci = Tr_over_Tci;
32 this->sqrt_Tr_Tci = sqrt(Tr_over_Tci);
33 };
34 AbstractCubicAlphaFunction(double a0, double Tr_over_Tci) : a0(a0), Tr_over_Tci(Tr_over_Tci), sqrt_Tr_Tci(sqrt(Tr_over_Tci)){};
35};
36
38class BasicMathiasCopemanAlphaFunction : public AbstractCubicAlphaFunction
39{
40 double m;
41 public:
42 BasicMathiasCopemanAlphaFunction(double a0, double m_ii, double Tr_over_Tci) : AbstractCubicAlphaFunction(a0, Tr_over_Tci) {
43 this->m = m_ii;
44 };
45 double term(double tau, std::size_t itau);
46};
47
49class TwuAlphaFunction : public AbstractCubicAlphaFunction
50{
51 public:
52 TwuAlphaFunction(double a0, double L, double M, double N, double Tr_over_Tci) : AbstractCubicAlphaFunction(a0, Tr_over_Tci) {
53 c.resize(3);
54 c[0] = L;
55 c[1] = M;
56 c[2] = N;
57 };
58 double term(double tau, std::size_t itau);
59};
60
62class MathiasCopemanAlphaFunction : public AbstractCubicAlphaFunction
63{
64 public:
65 MathiasCopemanAlphaFunction(double a0, double c1, double c2, double c3, double Tr_over_Tci) : AbstractCubicAlphaFunction(a0, Tr_over_Tci) {
66 c.resize(3);
67 c[0] = c1;
68 c[1] = c2;
69 c[2] = c3;
70 };
71 double term(double tau, std::size_t itau);
72};
73
74class AbstractCubic
75{
76
77 protected:
78 double rho_r,
79 T_r;
80 std::vector<double> Tc,
81 pc,
82 acentric;
83 double R_u;
84 double Delta_1,
85 Delta_2;
86 int N;
87 std::vector<std::vector<double>> k;
88 double cm;
89 std::vector<shared_ptr<AbstractCubicAlphaFunction>> alpha;
90 public:
100 AbstractCubic(std::vector<double> Tc, std::vector<double> pc, std::vector<double> acentric, double R_u, double Delta_1, double Delta_2,
101 std::vector<double> C1 = std::vector<double>(), std::vector<double> C2 = std::vector<double>(),
102 std::vector<double> C3 = std::vector<double>());
103 virtual ~AbstractCubic(){};
105 void set_alpha(const std::vector<double>& C1, const std::vector<double>& C2, const std::vector<double>& C3);
107 void set_alpha_function(std::size_t i, shared_ptr<AbstractCubicAlphaFunction>& acaf) {
108 alpha[i] = acaf;
109 };
111 shared_ptr<AbstractCubicAlphaFunction> get_alpha_function(std::size_t i) {
112 return alpha[i];
113 };
115 const std::vector<shared_ptr<AbstractCubicAlphaFunction>>& get_all_alpha_functions() {
116 return this->alpha;
117 };
119 void set_all_alpha_functions(const std::vector<shared_ptr<AbstractCubicAlphaFunction>>& alpha) {
120 this->alpha = alpha;
121 };
122
124 const std::vector<std::vector<double>>& get_kmat() {
125 return k;
126 };
128 void set_kmat(const std::vector<std::vector<double>>& k) {
129 this->k = k;
130 };
132 void set_kij(std::size_t i, std::size_t j, double val) {
133 k[i][j] = val;
134 k[j][i] = val;
135 }
137 double get_kij(std::size_t i, std::size_t j) {
138 return k[i][j];
139 }
141 std::vector<double>& get_Tc() {
142 return Tc;
143 }
145 std::vector<double>& get_pc() {
146 return pc;
147 }
149 std::vector<double>& get_acentric() {
150 return acentric;
151 }
153 double get_Delta_1() {
154 return Delta_1;
155 }
157 double get_Delta_2() {
158 return Delta_2;
159 }
161 double get_R_u() {
162 return R_u;
163 }
165 void set_Tr(double Tr) {
166 T_r = Tr;
167 for (std::size_t i = 0; i < alpha.size(); ++i) {
168 alpha[i]->set_Tr_over_Tci(T_r / Tc[i]);
169 }
170 }
172 void set_rhor(double rhor) {
173 rho_r = rhor;
174 }
176 double get_Tr() {
177 return T_r;
178 }
180 double get_rhor() {
181 return rho_r;
182 }
183
185 void set_C_MC(std::size_t i, double c1, double c2, double c3) {
186 alpha[i].reset(new MathiasCopemanAlphaFunction(a0_ii(i), c1, c2, c3, T_r / Tc[i]));
187 }
189 void set_C_Twu(std::size_t i, double L, double M, double N) {
190 alpha[i].reset(new TwuAlphaFunction(a0_ii(i), L, M, N, T_r / Tc[i]));
191 }
194 virtual double a0_ii(std::size_t i) = 0;
197 virtual double b0_ii(std::size_t i) = 0;
199 virtual double m_ii(std::size_t i) = 0;
200
202 virtual double alphar(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta);
204 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,
205 bool xN_independent);
207 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,
208 std::size_t j, bool xN_independent);
210 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,
211 std::size_t j, std::size_t k, bool xN_independent);
212
219 virtual double am_term(double tau, const std::vector<double>& x, std::size_t itau);
228 virtual double d_am_term_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent);
238 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);
249 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,
250 bool xN_independent);
251
256 virtual double bm_term(const std::vector<double>& x);
262 virtual double d_bm_term_dxi(const std::vector<double>& x, std::size_t i, bool xN_independent);
269 virtual double d2_bm_term_dxidxj(const std::vector<double>& x, std::size_t i, std::size_t j, bool xN_independent);
277 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);
281 virtual double cm_term();
283 void set_cm(double val) {
284 cm = val;
285 }
287 double get_cm() {
288 return cm;
289 }
290
292 virtual void set_Q_k(const size_t sgi, const double value) {
293 throw CoolProp::ValueError("set_Q_k not defined for AbstractCubic");
294 };
296 virtual double get_Q_k(const size_t sgi) const {
297 throw CoolProp::ValueError("get_Q_k not defined for AbstractCubic");
298 };
299
307 double aij_term(double tau, std::size_t i, std::size_t j, std::size_t itau);
314 double u_term(double tau, std::size_t i, std::size_t j, std::size_t itau);
320 double aii_term(double tau, std::size_t i, std::size_t itau);
321
329 double psi_minus(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta);
339 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);
350 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,
351 bool xN_independent);
363 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,
364 std::size_t k, bool xN_independent);
365
374 double PI_12(double delta, const std::vector<double>& x, std::size_t idelta);
383 double d_PI_12_dxi(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, bool xN_independent);
393 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);
404 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,
405 bool xN_independent);
406
413 double tau_times_a(double tau, const std::vector<double>& x, std::size_t itau);
422 double d_tau_times_a_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent);
432 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);
443 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,
444 bool xN_independent);
445
454 double psi_plus(double delta, const std::vector<double>& x, std::size_t idelta);
463 double d_psi_plus_dxi(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, bool xN_independent);
473 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);
484 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,
485 bool xN_independent);
486
495 double c_term(const std::vector<double>& x) {
496 return 1 / bm_term(x);
497 };
504 double d_c_term_dxi(const std::vector<double>& x, std::size_t i, bool xN_independent) {
505 return -d_bm_term_dxi(x, i, xN_independent) / pow(bm_term(x), 2);
506 };
514 double d2_c_term_dxidxj(const std::vector<double>& x, std::size_t i, std::size_t j, bool xN_independent) {
515 double bm = bm_term(x);
516 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))
517 / pow(bm, 3);
518 };
527 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) {
528 double bm = bm_term(x);
529 return 1 / pow(bm, 4)
530 * (2 * bm
531 * (d_bm_term_dxi(x, i, xN_independent) * d2_bm_term_dxidxj(x, j, k, xN_independent)
532 + d_bm_term_dxi(x, j, xN_independent) * d2_bm_term_dxidxj(x, i, k, xN_independent)
533 + d_bm_term_dxi(x, k, xN_independent) * d2_bm_term_dxidxj(x, i, j, xN_independent))
534 - pow(bm, 2) * d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent)
535 - 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));
536 };
537
548 double A_term(double delta, const std::vector<double>& x) {
549 double bm = bm_term(x);
550 double cm = cm_term();
551 return log((delta * rho_r * (Delta_1 * bm + cm) + 1) / (delta * rho_r * (Delta_2 * bm + cm) + 1));
552 };
560 double d_A_term_dxi(double delta, const std::vector<double>& x, std::size_t i, bool xN_independent) {
561 std::size_t idelta = 0;
562 return delta * rho_r * d_bm_term_dxi(x, i, xN_independent) * (Delta_1 - Delta_2) / PI_12(delta, x, idelta);
563 };
572 double d2_A_term_dxidxj(double delta, const std::vector<double>& x, std::size_t i, std::size_t j, bool xN_independent) {
573 std::size_t idelta = 0;
574 double PI12 = PI_12(delta, x, idelta);
575 return delta * rho_r * (Delta_1 - Delta_2) / pow(PI12, 2)
576 * (PI12 * d2_bm_term_dxidxj(x, i, j, xN_independent)
577 - d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_bm_term_dxi(x, i, xN_independent));
578 };
588 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) {
589 std::size_t idelta = 0;
590 double PI12 = PI_12(delta, x, idelta);
591 // The leading factor
592 double lead = delta * rho_r * (Delta_1 - Delta_2) / pow(PI12, 3);
593 return lead
594 * (-PI12
595 * (d_PI_12_dxi(delta, x, idelta, j, xN_independent) * d2_bm_term_dxidxj(x, i, k, xN_independent)
596 + d_PI_12_dxi(delta, x, idelta, k, xN_independent) * d2_bm_term_dxidxj(x, i, j, xN_independent)
597 + d_bm_term_dxi(x, i, xN_independent) * d2_PI_12_dxidxj(delta, x, idelta, j, k, xN_independent))
598 + pow(PI12, 2) * d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent)
599 + 2 * d_PI_12_dxi(delta, x, idelta, j, xN_independent) * d_PI_12_dxi(delta, x, idelta, k, xN_independent)
600 * d_bm_term_dxi(x, i, xN_independent));
601 };
602 // Allows to modify the unifac interaction parameters aij, bij and cij. Only for use with VTPR backend.
603 virtual void set_interaction_parameter(const std::size_t mgi1, const std::size_t mgi2, const std::string& parameter, const double value) {
604 throw CoolProp::NotImplementedError("set_interaction_parameter is not implemented for this backend");
605 }
606 // Allows to get the unifac interaction parameters aij, bij and cij. Only for use with VTPR backend.
607 virtual double get_interaction_parameter(const std::size_t mgi1, const std::size_t mgi2, const std::string& parameter) {
608 throw CoolProp::NotImplementedError("get_interaction_parameter is not implemented for this backend");
609 }
610};
611
612class PengRobinson : public AbstractCubic
613{
614 public:
615 PengRobinson(std::vector<double> Tc, std::vector<double> pc, std::vector<double> acentric, double R_u,
616 std::vector<double> C1 = std::vector<double>(), std::vector<double> C2 = std::vector<double>(),
617 std::vector<double> C3 = std::vector<double>())
618 : AbstractCubic(Tc, pc, acentric, R_u, 1 + sqrt(2.0), 1 - sqrt(2.0), C1, C2, C3) {
619 set_alpha(C1, C2, C3);
620 };
621
622 PengRobinson(double Tc, double pc, double acentric, double R_u)
623 : 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)) {
624 set_alpha(std::vector<double>(), std::vector<double>(), std::vector<double>());
625 };
626
627 double a0_ii(std::size_t i);
628 double b0_ii(std::size_t i);
629 double m_ii(std::size_t i);
630};
631
632class SRK : public AbstractCubic
633{
634 public:
635 SRK(std::vector<double> Tc, std::vector<double> pc, std::vector<double> acentric, double R_u, std::vector<double> C1 = std::vector<double>(),
636 std::vector<double> C2 = std::vector<double>(), std::vector<double> C3 = std::vector<double>())
637 : AbstractCubic(Tc, pc, acentric, R_u, 1, 0, C1, C2, C3) {
638 set_alpha(C1, C2, C3);
639 };
640 SRK(double Tc, double pc, double acentric, double R_u)
641 : AbstractCubic(std::vector<double>(1, Tc), std::vector<double>(1, pc), std::vector<double>(1, acentric), R_u, 1, 0) {
642 set_alpha(std::vector<double>(), std::vector<double>(), std::vector<double>());
643 };
644
645 double a0_ii(std::size_t i);
646 double b0_ii(std::size_t i);
647 double m_ii(std::size_t i);
648};
649
650#endif
Generated on Wed Jun 4 2025 04:12:26 for CoolProp by doxygen 1.9.4