CoolProp 6.8.1dev
An open-source fluid property and humid air property database
VTPRCubic.h
Go to the documentation of this file.
1//
2// VTPRCubic.h
3// CoolProp
4//
5// Created by Ian on 7/17/16.
6//
7//
8
9#include "GeneralizedCubic.h"
10#include "Exceptions.h"
11
12#ifndef VTPRCubic_h
13# define VTPRCubic_h
14
15class VTPRCubic : public PengRobinson
16{
17 private:
18 UNIFAC::UNIFACMixture unifaq;
19
20 public:
21 VTPRCubic(std::vector<double> Tc, std::vector<double> pc, std::vector<double> acentric, double R_u,
22 const UNIFACLibrary::UNIFACParameterLibrary& lib)
23 : PengRobinson(Tc, pc, acentric, R_u), unifaq(lib, T_r){};
24
25 VTPRCubic(double Tc, double pc, double acentric, double R_u, const UNIFACLibrary::UNIFACParameterLibrary& lib)
26 : PengRobinson(std::vector<double>(1, Tc), std::vector<double>(1, pc), std::vector<double>(1, acentric), R_u), unifaq(lib, T_r){};
27
29 UNIFAC::UNIFACMixture& get_unifaq() {
30 return unifaq;
31 }
32
34 double gE_R_RT(double tau, const std::vector<double>& x, std::size_t itau) {
35 double summer = 0;
36 for (std::size_t i = 0; i < x.size(); ++i) {
37 summer += x[i] * unifaq.ln_gamma_R(tau, i, itau);
38 }
39 return summer;
40 }
41 double d_gE_R_RT_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent) {
42 if (xN_independent) {
43 return unifaq.ln_gamma_R(tau, i, itau);
44 } else {
45 return unifaq.ln_gamma_R(tau, i, itau) - unifaq.ln_gamma_R(tau, N - 1, itau);
46 }
47 }
48 double gE_R(double tau, const std::vector<double>& x, std::size_t itau) {
49 if (x.size() == 1) {
50 return 0.;
51 } else {
52 switch (itau) {
53 case 0: {
54 return R_u * T_r / tau * gE_R_RT(tau, x, 0);
55 }
56 case 1: {
57 return R_u * T_r / tau * (-gE_R_RT(tau, x, 0) / tau + gE_R_RT(tau, x, 1));
58 }
59 case 2: {
60 return R_u * T_r / tau * (2 * gE_R_RT(tau, x, 0) / powInt(tau, 2) - 2 * gE_R_RT(tau, x, 1) / tau + gE_R_RT(tau, x, 2));
61 }
62 case 3: {
63 return R_u * T_r / tau
64 * (-6 * gE_R_RT(tau, x, 0) / powInt(tau, 3) + 6 * gE_R_RT(tau, x, 1) / powInt(tau, 2) - 3 * gE_R_RT(tau, x, 2) / tau
65 + gE_R_RT(tau, x, 3));
66 }
67 case 4: {
68 return R_u * T_r / tau
69 * (24 * gE_R_RT(tau, x, 0) / powInt(tau, 4) - 24 * gE_R_RT(tau, x, 1) / powInt(tau, 3)
70 + 12 * gE_R_RT(tau, x, 2) / powInt(tau, 2) - 4 * gE_R_RT(tau, x, 3) / tau + gE_R_RT(tau, x, 4));
71 }
72 default:
73 throw CoolProp::ValueError(format("itau (%d) is invalid", itau));
74 }
75 }
76 }
77 double d_gE_R_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent) {
78 if (x.size() == 1) {
79 return 0.;
80 } else {
81 switch (itau) {
82 case 0: {
83 return R_u * T_r / tau * d_gE_R_RT_dxi(tau, x, 0, i, xN_independent);
84 }
85 case 1: {
86 return R_u * T_r / tau * (-d_gE_R_RT_dxi(tau, x, 0, i, xN_independent) / tau + d_gE_R_RT_dxi(tau, x, 1, i, xN_independent));
87 }
88 case 2: {
89 return R_u * T_r / tau
90 * (2 * d_gE_R_RT_dxi(tau, x, 0, i, xN_independent) / powInt(tau, 2) - 2 * d_gE_R_RT_dxi(tau, x, 1, i, xN_independent) / tau
91 + d_gE_R_RT_dxi(tau, x, 2, i, xN_independent));
92 }
93 case 3: {
94 return R_u * T_r / tau
95 * (-6 * d_gE_R_RT_dxi(tau, x, 0, i, xN_independent) / powInt(tau, 3)
96 + 6 * d_gE_R_RT_dxi(tau, x, 1, i, xN_independent) / powInt(tau, 2)
97 - 3 * d_gE_R_RT_dxi(tau, x, 2, i, xN_independent) / tau + d_gE_R_RT_dxi(tau, x, 3, i, xN_independent));
98 }
99 case 4: {
100 return R_u * T_r / tau
101 * (24 * d_gE_R_RT_dxi(tau, x, 0, i, xN_independent) / powInt(tau, 4)
102 - 24 * d_gE_R_RT_dxi(tau, x, 1, i, xN_independent) / powInt(tau, 3)
103 + 12 * d_gE_R_RT_dxi(tau, x, 2, i, xN_independent) / powInt(tau, 2)
104 - 4 * d_gE_R_RT_dxi(tau, x, 3, i, xN_independent) / tau + d_gE_R_RT_dxi(tau, x, 4, i, xN_independent));
105 }
106 default:
107 throw CoolProp::ValueError(format("itau (%d) is invalid", itau));
108 }
109 }
110 }
111 double am_term(double tau, const std::vector<double>& x, std::size_t itau) {
112 return bm_term(x) * (sum_xi_aii_bii(tau, x, itau) + gE_R(tau, x, itau) / (-0.53087));
113 }
114 double sum_xi_aii_bii(double tau, const std::vector<double>& x, std::size_t itau) {
115 double summeram = 0;
116 for (int i = 0; i < N; ++i) {
117 summeram += x[i] * aii_term(tau, i, itau) / b0_ii(i);
118 }
119 return summeram;
120 }
121 double d_sum_xi_aii_bii_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent) {
122 if (xN_independent) {
123 return aii_term(tau, i, itau) / b0_ii(i);
124 } else {
125 return aii_term(tau, i, itau) / b0_ii(i) - aii_term(tau, N - 1, itau) / b0_ii(N - 1);
126 }
127 }
128 double d_am_term_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent) {
129 return d_bm_term_dxi(x, i, xN_independent) * (sum_xi_aii_bii(tau, x, itau) + gE_R(tau, x, itau) / (-0.53087))
130 + bm_term(x) * (d_sum_xi_aii_bii_dxi(tau, x, itau, i, xN_independent) + d_gE_R_dxi(tau, x, itau, i, xN_independent) / (-0.53087));
131 }
132 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) {
133 return d2_bm_term_dxidxj(x, i, j, xN_independent) * (sum_xi_aii_bii(tau, x, itau) + gE_R(tau, x, itau) / (-0.53087))
134 + d_bm_term_dxi(x, i, xN_independent)
135 * (d_sum_xi_aii_bii_dxi(tau, x, itau, i, xN_independent) + d_gE_R_dxi(tau, x, itau, i, xN_independent) / (-0.53087))
136 + d_bm_term_dxi(x, j, xN_independent)
137 * (d_sum_xi_aii_bii_dxi(tau, x, itau, i, xN_independent) + d_gE_R_dxi(tau, x, itau, i, xN_independent) / (-0.53087));
138 }
139 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,
140 bool xN_independent) {
141 return d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent) * (sum_xi_aii_bii(tau, x, itau) + gE_R(tau, x, itau) / (-0.53087))
142 + d2_bm_term_dxidxj(x, i, k, xN_independent)
143 * (d_sum_xi_aii_bii_dxi(tau, x, itau, i, xN_independent) + d_gE_R_dxi(tau, x, itau, i, xN_independent) / (-0.53087))
144 + d2_bm_term_dxidxj(x, j, k, xN_independent)
145 * (d_sum_xi_aii_bii_dxi(tau, x, itau, i, xN_independent) + d_gE_R_dxi(tau, x, itau, i, xN_independent) / (-0.53087));
146 }
147
148 double bm_term(const std::vector<double>& x) {
149 double summerbm = 0;
150 for (int i = 0; i < N; ++i) {
151 for (int j = 0; j < N; ++j) {
152 summerbm += x[i] * x[j] * bij_term(i, j);
153 }
154 }
155 return summerbm;
156 }
157 double bij_term(std::size_t i, std::size_t j) {
158 return pow((pow(b0_ii(i), 0.75) + pow(b0_ii(j), 0.75)) / 2.0, 4.0 / 3.0);
159 }
160 double d_bm_term_dxi(const std::vector<double>& x, std::size_t i, bool xN_independent) {
161 double summer = 0;
162 if (xN_independent) {
163 for (int j = N - 1; j >= 0; --j) {
164 summer += x[j] * bij_term(i, j);
165 }
166 return 2 * summer;
167 } else {
168 for (int k = N - 2; k >= 0; --k) {
169 summer += x[k] * (bij_term(i, k) - bij_term(k, N - 1));
170 }
171 return 2 * (summer + x[N - 1] * (bij_term(N - 1, i) - bij_term(N - 1, N - 1)));
172 }
173 }
174 double d2_bm_term_dxidxj(const std::vector<double>& x, std::size_t i, std::size_t j, bool xN_independent) {
175 if (xN_independent) {
176 return 2 * bij_term(i, j);
177 } else {
178 return 2 * (bij_term(i, j) - bij_term(j, N - 1) - bij_term(N - 1, i) + bij_term(N - 1, N - 1));
179 }
180 }
181 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) {
182 return 0;
183 }
184 // Allows to modify the unifac interaction parameters aij, bij and cij
185 void set_interaction_parameter(const std::size_t mgi1, const std::size_t mgi2, const std::string& parameter, const double value) {
186 unifaq.set_interaction_parameter(mgi1, mgi2, parameter, value);
187 }
188
189 // Allows to modify the surface parameter Q_k of the sub group sgi
190 void set_Q_k(const size_t sgi, const double value) {
191 unifaq.set_Q_k(sgi, value);
192 }
193
194 // Get the surface parameter Q_k of the sub group sgi
195 double get_Q_k(const size_t sgi) const {
196 return unifaq.get_Q_k(sgi);
197 }
198
199 // Allows to modify the unifac interaction parameters aij, bij and cij
200 double get_interaction_parameter(const std::size_t mgi1, const std::size_t mgi2, const std::string& parameter) {
201 return unifaq.get_interaction_parameter(mgi1, mgi2, parameter);
202 }
203 std::vector<double> ln_fugacity_coefficient(const std::vector<double>& z, double rhomolar, double p, double T) {
204 double v = 1 / rhomolar;
205 // Common terms for all components
206 double tau = get_Tr() / T;
207 double b = bm_term(z);
208 double c = cm_term();
209 double R = get_R_u();
210 std::vector<double> ln_phi;
211 double bracket = log((v + c + (1 + sqrt(2.0)) * b) / (v + c + (1 - sqrt(2.0)) * b));
212 for (std::size_t i = 0; i < z.size(); ++i) {
213 double summer1 = 0;
214 for (std::size_t j = 0; j < z.size(); ++j) {
215 summer1 += z[j] * bij_term(i, j);
216 }
217 double a_i_over_b_i = aii_term(tau, i, 0) / b0_ii(i);
218 double c_i = 0; // TODO: fix this, allow for volume translation
219 double _ln_phi = (2 / b * summer1 - 1) * (p * (v + c) / (R * T) - 1) - p * c_i / (R * T) - log(p * (v + c - b) / (R * T))
220 - 1.0 / (2.0 * sqrt(2.0) * R * T) * (a_i_over_b_i + R * T * unifaq.ln_gamma_R(tau, i, 0) / -0.53087) * bracket;
221 ln_phi.push_back(_ln_phi);
222 }
223 return ln_phi;
224 }
225};
226
227#endif /* VTPRCubic_h */
Generated on Thu Jun 5 2025 04:10:24 for CoolProp by doxygen 1.9.4