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