CoolProp  6.6.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 "GeneralizedCubic.h"
10 #include "Exceptions.h"
11 
12 #ifndef VTPRCubic_h
13 # define VTPRCubic_h
14 
15 class 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 Wed Nov 29 2023 14:26:23 for CoolProp by doxygen 1.9.1