CoolProp  6.6.0
An open-source fluid property and humid air property database
ReducingFunctions.cpp
Go to the documentation of this file.
1 #include "ReducingFunctions.h"
2 
3 namespace CoolProp {
4 
5 CoolPropDbl ReducingFunction::d_ndTrdni_dxj__constxi(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
6  x_N_dependency_flag xN_flag) const {
7  if (xN_flag == XN_INDEPENDENT) {
8  CoolPropDbl s = 0;
9  for (std::size_t k = 0; k < N; k++) {
10  s += x[k] * d2Trdxidxj(x, j, k, xN_flag);
11  }
12  return d2Trdxidxj(x, i, j, xN_flag) - dTrdxi__constxj(x, j, xN_flag) - s;
13  } else if (xN_flag == XN_DEPENDENT) {
14  if (j == N - 1) {
15  return 0;
16  }
17  if (N == 0) {
18  return 0;
19  }
20  CoolPropDbl s = 0;
21  for (std::size_t k = 0; k < N - 1; k++) {
22  s += x[k] * d2Trdxidxj(x, k, j, xN_flag);
23  }
24  CoolPropDbl val = d2Trdxidxj(x, j, i, xN_flag) - dTrdxi__constxj(x, j, xN_flag) - s;
25  return val;
26  } else {
27  throw ValueError(format("xN dependency flag invalid"));
28  }
29 }
30 CoolPropDbl ReducingFunction::d2_ndTrdni_dxj_dxk__constxi(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, std::size_t k,
31  x_N_dependency_flag xN_flag) const {
32  if (xN_flag == XN_INDEPENDENT) {
33  CoolPropDbl s = 0;
34  for (std::size_t m = 0; m < N; m++) {
35  s += x[m] * d3Trdxidxjdxk(x, j, k, m, xN_flag);
36  }
37  return d3Trdxidxjdxk(x, i, j, k, xN_flag) - 2 * d2Trdxidxj(x, j, k, xN_flag) - s;
38  } else if (xN_flag == XN_DEPENDENT) {
39  if (N == 0) {
40  return 0;
41  }
42  if (j == N - 1) {
43  return 0;
44  }
45  CoolPropDbl s = 0;
46  for (std::size_t m = 0; m < N - 1; m++) {
47  s += x[m] * d3Trdxidxjdxk(x, k, j, m, xN_flag);
48  }
49  CoolPropDbl val = d3Trdxidxjdxk(x, i, j, k, xN_flag) - d2Trdxidxj(x, j, k, xN_flag) - s;
50  return val;
51  } else {
52  throw ValueError(format("xN dependency flag invalid"));
53  }
54 }
55 CoolPropDbl ReducingFunction::d_ndrhorbardni_dxj__constxi(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
56  x_N_dependency_flag xN_flag) const {
57  CoolPropDbl s = 0;
58  for (std::size_t k = 0; k < N; k++) {
59  s += x[k] * d2rhormolardxidxj(x, j, k, xN_flag);
60  }
61  return d2rhormolardxidxj(x, j, i, xN_flag) - drhormolardxi__constxj(x, j, xN_flag) - s;
62 }
63 CoolPropDbl ReducingFunction::d2_ndrhorbardni_dxj_dxk__constxi(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, std::size_t k,
64  x_N_dependency_flag xN_flag) const {
65  CoolPropDbl s = 0;
66  for (std::size_t m = 0; m < N; m++) {
67  s += x[m] * d3rhormolardxidxjdxk(x, j, k, m, xN_flag);
68  }
69  return d3rhormolardxidxjdxk(x, i, j, k, xN_flag) - 2 * d2rhormolardxidxj(x, j, k, xN_flag) - s;
70 }
71 CoolPropDbl ReducingFunction::ndrhorbardni__constnj(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
72  if (xN_flag == XN_INDEPENDENT) {
73  CoolPropDbl summer_term1 = 0;
74  for (std::size_t j = 0; j < N; j++) {
75  summer_term1 += x[j] * drhormolardxi__constxj(x, j, xN_flag);
76  }
77  return drhormolardxi__constxj(x, i, xN_flag) - summer_term1;
78  } else if (xN_flag == XN_DEPENDENT) {
79  CoolPropDbl summer_term1 = 0;
80  if (N == 0) {
81  return 0;
82  }
83  for (std::size_t k = 0; k < N - 1; ++k) {
84  summer_term1 += x[k] * drhormolardxi__constxj(x, k, xN_flag);
85  }
86  return drhormolardxi__constxj(x, i, xN_flag) - summer_term1;
87  } else {
88  throw ValueError(format("xN dependency flag invalid"));
89  }
90 }
91 CoolPropDbl ReducingFunction::ndTrdni__constnj(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
92  if (xN_flag == XN_INDEPENDENT) {
93  // GERG Equation 7.54
94  CoolPropDbl summer_term1 = 0;
95  for (std::size_t j = 0; j < N; j++) {
96  summer_term1 += x[j] * dTrdxi__constxj(x, j, xN_flag);
97  }
98  return dTrdxi__constxj(x, i, xN_flag) - summer_term1;
99  } else if (xN_flag == XN_DEPENDENT) {
100  CoolPropDbl summer_term1 = 0;
101  if (N == 0) {
102  return 0;
103  }
104  for (std::size_t k = 0; k < N - 1; ++k) {
105  summer_term1 += x[k] * dTrdxi__constxj(x, k, xN_flag);
106  }
107  return dTrdxi__constxj(x, i, xN_flag) - summer_term1;
108  } else {
109  throw ValueError(format("xN dependency flag invalid"));
110  }
111 }
112 CoolPropDbl ReducingFunction::PSI_rho(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
113  return 1 - 1 / rhormolar(x) * ndrhorbardni__constnj(x, i, xN_flag);
114 }
115 CoolPropDbl ReducingFunction::d_PSI_rho_dxj(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) const {
116  return -1 / rhormolar(x) * (d_ndrhorbardni_dxj__constxi(x, i, j, xN_flag) - drhormolardxi__constxj(x, j, xN_flag) * (1 - PSI_rho(x, i, xN_flag)));
117 }
118 CoolPropDbl ReducingFunction::d2_PSI_rho_dxj_dxk(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, std::size_t k,
119  x_N_dependency_flag xN_flag) const {
120  double line1 = d2_ndrhorbardni_dxj_dxk__constxi(x, i, j, k, xN_flag);
121  double line2 = -1 / rhormolar(x) * drhormolardxi__constxj(x, k, xN_flag) * d_ndrhorbardni_dxj__constxi(x, i, j, xN_flag);
122  double line3 = drhormolardxi__constxj(x, j, xN_flag) * d_PSI_rho_dxj(x, i, k, xN_flag);
123  double line4 =
124  -(d2rhormolardxidxj(x, j, k, xN_flag) - 1 / rhormolar(x) * drhormolardxi__constxj(x, k, xN_flag) * drhormolardxi__constxj(x, j, xN_flag))
125  * (1 - PSI_rho(x, i, xN_flag));
126  return -1 / rhormolar(x) * (line1 + line2 + line3 + line4);
127 }
128 CoolPropDbl ReducingFunction::PSI_T(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
129  return 1 / Tr(x) * ndTrdni__constnj(x, i, xN_flag);
130 }
131 CoolPropDbl ReducingFunction::d_PSI_T_dxj(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) const {
132  return 1 / Tr(x) * (d_ndTrdni_dxj__constxi(x, i, j, xN_flag) - dTrdxi__constxj(x, j, xN_flag) * PSI_T(x, i, xN_flag));
133 }
134 CoolPropDbl ReducingFunction::d2_PSI_T_dxj_dxk(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, std::size_t k,
135  x_N_dependency_flag xN_flag) const {
136  double line1 = d2_ndTrdni_dxj_dxk__constxi(x, i, j, k, xN_flag);
137  double line2 = -1 / Tr(x) * dTrdxi__constxj(x, k, xN_flag) * d_ndTrdni_dxj__constxi(x, i, j, xN_flag);
138  double line3 = -dTrdxi__constxj(x, j, xN_flag) * d_PSI_T_dxj(x, i, k, xN_flag);
139  double line4 =
140  -(d2Trdxidxj(x, j, k, xN_flag) - 1 / Tr(x) * dTrdxi__constxj(x, k, xN_flag) * dTrdxi__constxj(x, j, xN_flag)) * PSI_T(x, i, xN_flag);
141  return 1 / Tr(x) * (line1 + line2 + line3 + line4);
142 }
143 
144 CoolPropDbl GERG2008ReducingFunction::Tr(const std::vector<CoolPropDbl>& x) const {
145  return Yr(x, beta_T, gamma_T, T_c, Yc_T);
146 }
147 CoolPropDbl GERG2008ReducingFunction::dTr_dbetaT(const std::vector<CoolPropDbl>& x) const {
148  return dYr_dbeta(x, beta_T, gamma_T, T_c, Yc_T);
149 }
150 CoolPropDbl GERG2008ReducingFunction::dTr_dgammaT(const std::vector<CoolPropDbl>& x) const {
151  return dYr_dgamma(x, beta_T, gamma_T, T_c, Yc_T);
152 }
153 CoolPropDbl GERG2008ReducingFunction::d2Tr_dxidgammaT(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
154  return d2Yrdxidgamma(x, i, beta_T, gamma_T, T_c, Yc_T, xN_flag);
155 };
156 CoolPropDbl GERG2008ReducingFunction::d2Tr_dxidbetaT(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
157  return d2Yrdxidbeta(x, i, beta_T, gamma_T, T_c, Yc_T, xN_flag);
158 };
159 
160 CoolPropDbl GERG2008ReducingFunction::dTrdxi__constxj(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
161  return dYrdxi__constxj(x, i, beta_T, gamma_T, T_c, Yc_T, xN_flag);
162 }
163 CoolPropDbl GERG2008ReducingFunction::d2Trdxi2__constxj(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
164  return d2Yrdxi2__constxj(x, i, beta_T, gamma_T, T_c, Yc_T, xN_flag);
165 }
166 CoolPropDbl GERG2008ReducingFunction::d2Trdxidxj(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) const {
167  return d2Yrdxidxj(x, i, j, beta_T, gamma_T, T_c, Yc_T, xN_flag);
168 }
169 CoolPropDbl GERG2008ReducingFunction::d3Trdxidxjdxk(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, std::size_t k,
170  x_N_dependency_flag xN_flag) const {
171  return d3Yrdxidxjdxk(x, i, j, k, beta_T, gamma_T, T_c, Yc_T, xN_flag);
172 }
173 CoolPropDbl GERG2008ReducingFunction::rhormolar(const std::vector<CoolPropDbl>& x) const {
174  return 1 / Yr(x, beta_v, gamma_v, v_c, Yc_v);
175 }
176 
177 CoolPropDbl GERG2008ReducingFunction::d2rhormolar_dxidgammaV(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
178  CoolPropDbl rhor = rhormolar(x);
179  return -rhor * rhor * d2vrmolar_dxidgammaV(x, i, xN_flag)
180  + 2 * POW3(rhor) * dvrmolardxi__constxj(x, i, xN_flag) * dYr_dgamma(x, beta_v, gamma_v, v_c, Yc_v);
181 }
182 CoolPropDbl GERG2008ReducingFunction::d2rhormolar_dxidbetaV(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
183  CoolPropDbl rhor = rhormolar(x);
184  return -rhor * rhor * d2vrmolar_dxidbetaV(x, i, xN_flag)
185  + 2 * POW3(rhor) * dvrmolardxi__constxj(x, i, xN_flag) * dYr_dbeta(x, beta_v, gamma_v, v_c, Yc_v);
186 }
187 CoolPropDbl GERG2008ReducingFunction::d2vrmolar_dxidgammaV(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
188  return d2Yrdxidgamma(x, i, beta_v, gamma_v, v_c, Yc_v, xN_flag);
189 };
190 CoolPropDbl GERG2008ReducingFunction::d2vrmolar_dxidbetaV(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
191  return d2Yrdxidbeta(x, i, beta_v, gamma_v, v_c, Yc_v, xN_flag);
192 };
193 CoolPropDbl GERG2008ReducingFunction::drhormolar_dgammaV(const std::vector<CoolPropDbl>& x) const {
194  CoolPropDbl rhor = rhormolar(x);
195  return -rhor * rhor * dYr_dgamma(x, beta_v, gamma_v, v_c, Yc_v);
196 }
197 CoolPropDbl GERG2008ReducingFunction::drhormolar_dbetaV(const std::vector<CoolPropDbl>& x) const {
198  CoolPropDbl rhor = rhormolar(x);
199  return -rhor * rhor * dYr_dbeta(x, beta_v, gamma_v, v_c, Yc_v);
200 }
201 CoolPropDbl GERG2008ReducingFunction::drhormolardxi__constxj(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
202  CoolPropDbl rhor = rhormolar(x);
203  return -rhor * rhor * dvrmolardxi__constxj(x, i, xN_flag);
204 }
205 CoolPropDbl GERG2008ReducingFunction::dvrmolardxi__constxj(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
206  return dYrdxi__constxj(x, i, beta_v, gamma_v, v_c, Yc_v, xN_flag);
207 }
208 CoolPropDbl GERG2008ReducingFunction::d2vrmolardxi2__constxj(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
209  return d2Yrdxi2__constxj(x, i, beta_v, gamma_v, v_c, Yc_v, xN_flag);
210 }
211 CoolPropDbl GERG2008ReducingFunction::d2vrmolardxidxj(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
212  x_N_dependency_flag xN_flag) const {
213  return d2Yrdxidxj(x, i, j, beta_v, gamma_v, v_c, Yc_v, xN_flag);
214 }
215 CoolPropDbl GERG2008ReducingFunction::d3vrmolardxidxjdxk(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, std::size_t k,
216  x_N_dependency_flag xN_flag) const {
217  return d3Yrdxidxjdxk(x, i, j, k, beta_v, gamma_v, v_c, Yc_v, xN_flag);
218 }
219 CoolPropDbl GERG2008ReducingFunction::d2rhormolardxi2__constxj(const std::vector<CoolPropDbl>& x, std::size_t i, x_N_dependency_flag xN_flag) const {
220  CoolPropDbl rhor = this->rhormolar(x);
221  CoolPropDbl dvrbardxi = this->dvrmolardxi__constxj(x, i, xN_flag);
222  return 2 * pow(rhor, 3) * pow(dvrbardxi, 2) - pow(rhor, 2) * this->d2vrmolardxi2__constxj(x, i, xN_flag);
223 }
224 CoolPropDbl GERG2008ReducingFunction::d2rhormolardxidxj(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j,
225  x_N_dependency_flag xN_flag) const {
226  double rhor = this->rhormolar(x);
227  double dvrbardxi = this->dvrmolardxi__constxj(x, i, xN_flag);
228  double dvrbardxj = this->dvrmolardxi__constxj(x, j, xN_flag);
229  return 2 * pow(rhor, 3) * dvrbardxi * dvrbardxj - pow(rhor, 2) * this->d2vrmolardxidxj(x, i, j, xN_flag);
230 }
231 CoolPropDbl GERG2008ReducingFunction::d3rhormolardxidxjdxk(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, std::size_t k,
232  x_N_dependency_flag xN_flag) const {
233  double rhor = this->rhormolar(x);
234  double line1 = -pow(rhor, 2) * d3vrmolardxidxjdxk(x, i, j, k, xN_flag);
235  double line2 = 2 * pow(rhor, 3)
236  * (dvrmolardxi__constxj(x, k, xN_flag) * d2vrmolardxidxj(x, i, j, xN_flag)
237  + dvrmolardxi__constxj(x, j, xN_flag) * d2vrmolardxidxj(x, i, k, xN_flag)
238  + dvrmolardxi__constxj(x, i, xN_flag) * d2vrmolardxidxj(x, j, k, xN_flag));
239  double line3 =
240  -6 * pow(rhor, 4) * dvrmolardxi__constxj(x, i, xN_flag) * dvrmolardxi__constxj(x, j, xN_flag) * dvrmolardxi__constxj(x, k, xN_flag);
241  return line1 + line2 + line3;
242 }
243 
244 CoolPropDbl GERG2008ReducingFunction::Yr(const std::vector<CoolPropDbl>& x, const STLMatrix& beta, const STLMatrix& gamma, const STLMatrix& Y_c_ij,
245  const std::vector<CoolPropDbl>& Yc) const {
246  CoolPropDbl Yr = 0;
247  for (std::size_t i = 0; i < N; i++) {
248  double xi = x[i];
249  Yr += xi * xi * Yc[i];
250 
251  // The last term is only used for the pure component, as it is sum_{i=1}^{N-1}sum_{j=1}^{N}
252  if (i == N - 1) {
253  break;
254  }
255 
256  for (std::size_t j = i + 1; j < N; j++) {
257  Yr += c_Y_ij(i, j, beta, gamma, Y_c_ij) * f_Y_ij(x, i, j, beta);
258  }
259  }
260  return Yr;
261 }
262 
263 CoolPropDbl GERG2008ReducingFunction::dYr_dgamma(const std::vector<CoolPropDbl>& x, const STLMatrix& beta, const STLMatrix& gamma,
264  const STLMatrix& Y_c_ij, const std::vector<CoolPropDbl>& Yc) const {
265  CoolPropDbl dYr_dgamma = 0;
266  for (std::size_t i = 0; i < N; i++) {
267  // The last term is only used for the pure component, as it is sum_{i=1}^{N-1}sum_{j=1}^{N}
268  if (i == N - 1) {
269  break;
270  }
271  for (std::size_t j = i + 1; j < N; j++) {
272  dYr_dgamma += 2 * beta[i][j] * Y_c_ij[i][j] * f_Y_ij(x, i, j, beta);
273  }
274  }
275  return dYr_dgamma;
276 }
277 CoolPropDbl GERG2008ReducingFunction::dYr_dbeta(const std::vector<CoolPropDbl>& x, const STLMatrix& beta, const STLMatrix& gamma,
278  const STLMatrix& Y_c_ij, const std::vector<CoolPropDbl>& Yc) const {
279  CoolPropDbl dYr_dbeta = 0;
280  for (std::size_t i = 0; i < N; i++) {
281  // The last term is only used for the pure component, as it is sum_{i=1}^{N-1}sum_{j=1}^{N}
282  if (i == N - 1) {
283  break;
284  }
285 
286  for (std::size_t j = i + 1; j < N; j++) {
287  double xj = x[j], xi = x[i], beta_Y = beta[i][j], beta_Y_squared = beta_Y * beta_Y;
288  if (std::abs(xi) < 10 * DBL_EPSILON && std::abs(xj) < 10 * DBL_EPSILON) {
289  return 0;
290  }
291  double dfYij_dbeta = xi * xj * (-(xi + xj) * (2 * beta_Y * xi)) / pow(beta_Y_squared * xi + xj, 2);
292  dYr_dbeta += c_Y_ij(i, j, beta, gamma, Y_c_ij) * dfYij_dbeta + f_Y_ij(x, i, j, beta) * 2 * gamma[i][j] * Y_c_ij[i][j];
293  }
294  }
295  return dYr_dbeta;
296 }
297 
298 CoolPropDbl GERG2008ReducingFunction::dYrdxi__constxj(const std::vector<CoolPropDbl>& x, std::size_t i, const STLMatrix& beta, const STLMatrix& gamma,
299  const STLMatrix& Y_c_ij, const std::vector<CoolPropDbl>& Yc,
300  x_N_dependency_flag xN_flag) const {
301  if (xN_flag == XN_INDEPENDENT) {
302  // See Table B9 from Kunz Wagner 2012 (GERG 2008)
303  CoolPropDbl xi = x[i];
304  CoolPropDbl dYr_dxi = 2 * xi * Yc[i];
305  for (std::size_t k = 0; k < i; k++) {
306  dYr_dxi += c_Y_ij(k, i, beta, gamma, Y_c_ij) * dfYkidxi__constxk(x, k, i, beta);
307  }
308  for (std::size_t k = i + 1; k < N; k++) {
309  dYr_dxi += c_Y_ij(i, k, beta, gamma, Y_c_ij) * dfYikdxi__constxk(x, i, k, beta);
310  }
311  return dYr_dxi;
312  } else if (xN_flag == XN_DEPENDENT) {
313  // Table S1 from Gernert, 2014, supplemental information
314  if (i == N - 1) {
315  return 0.0;
316  }
317  CoolPropDbl dYr_dxi = 2 * x[i] * Yc[i] - 2 * x[N - 1] * Yc[N - 1];
318  for (std::size_t k = 0; k < i; k++) {
319  dYr_dxi += c_Y_ij(k, i, beta, gamma, Y_c_ij) * dfYkidxi__constxk(x, k, i, beta);
320  }
321  for (std::size_t k = i + 1; k < N - 1; k++) {
322  dYr_dxi += c_Y_ij(i, k, beta, gamma, Y_c_ij) * dfYikdxi__constxk(x, i, k, beta);
323  }
324  double beta_Y_iN = beta[i][N - 1], xN = x[N - 1];
325  dYr_dxi += c_Y_ij(i, N - 1, beta, gamma, Y_c_ij)
326  * (xN * (x[i] + xN) / (pow(beta_Y_iN, 2) * x[i] + xN)
327  + (1 - beta_Y_iN * beta_Y_iN) * x[i] * xN * xN / POW2(beta_Y_iN * beta_Y_iN * x[i] + xN));
328  for (std::size_t k = 0; k < N - 1; ++k) {
329  double beta_Y_kN = beta[k][N - 1], xk = x[k], beta_Y_kN_squared = beta_Y_kN * beta_Y_kN;
330  dYr_dxi +=
331  c_Y_ij(k, N - 1, beta, gamma, Y_c_ij)
332  * (-xk * (xk + xN) / (beta_Y_kN_squared * xk + xN) + (1 - beta_Y_kN_squared) * xN * xk * xk / POW2(beta_Y_kN_squared * xk + xN));
333  }
334  return dYr_dxi;
335  } else {
336  throw ValueError(format("xN dependency flag invalid"));
337  }
338 }
339 CoolPropDbl GERG2008ReducingFunction::d2Yrdxidbeta(const std::vector<CoolPropDbl>& x, std::size_t i, const STLMatrix& beta, const STLMatrix& gamma,
340  const STLMatrix& Y_c_ij, const std::vector<CoolPropDbl>& Yc, x_N_dependency_flag xN_flag) const {
341  if (xN_flag == XN_INDEPENDENT) {
342  // See Table B9 from Kunz Wagner 2012 (GERG 2008)
343  CoolPropDbl xi = x[i];
344  CoolPropDbl deriv = 0;
345  for (std::size_t k = 0; k < i; k++) {
346  /*
347  Sympy code:
348  x_k,x_i,beta_Y = symbols('x_k,x_i,beta_Y')
349  dfkidxi = x_k*(x_k+x_i)/(beta_Y**2*x_k+x_i) + x_k*x_i/(beta_Y**2*x_k+x_i)*(1-(x_k+x_i)/(beta_Y**2*x_k+x_i))
350  simplify(diff(dfkidxi, beta_Y))
351  */
352  double xk = x[k], beta_Y = beta[k][i], beta_Y_squared = beta_Y * beta_Y;
353  double d2fYkidxidbeta = 2 * beta_Y * pow(xk, 2) * (xi * (xi + xk * (-beta_Y_squared + 1) + xk) - (xi + xk) * (beta_Y_squared * xk + xi))
354  / pow(beta_Y_squared * xk + xi, 3);
355  deriv += c_Y_ij(k, i, beta, gamma, Y_c_ij) * d2fYkidxidbeta + dfYkidxi__constxk(x, k, i, beta) * 2 * gamma[k][i] * Y_c_ij[k][i];
356  }
357  for (std::size_t k = i + 1; k < N; k++) {
358  /*
359  x_k,x_i,beta_Y = symbols('x_k,x_i,beta_Y')
360  dfikdxi = x_k*(x_i+x_k)/(beta_Y**2*x_i+x_k) + x_i*x_k/(beta_Y**2*x_i+x_k)*(1-beta_Y**2*(x_i+x_k)/(beta_Y**2*x_i+x_k))
361  print(ccode(simplify(diff(dfikdxi, beta_Y))))
362  */
363  double xk = x[k], beta_Y = beta[i][k], beta_Y_squared = beta_Y * beta_Y;
364  double d2fYikdxidbeta =
365  2 * beta_Y * xi * xk
366  * (xi * (-beta_Y_squared * xi + beta_Y_squared * (xi + xk) - xk) - xk * (xi + xk) - (xi + xk) * (beta_Y_squared * xi + xk))
367  / pow(beta_Y_squared * xi + xk, 3);
368  deriv += c_Y_ij(i, k, beta, gamma, Y_c_ij) * d2fYikdxidbeta + dfYikdxi__constxk(x, i, k, beta) * 2 * gamma[i][k] * Y_c_ij[i][k];
369  }
370  return deriv;
371  } else if (xN_flag == XN_DEPENDENT) {
372  throw NotImplementedError("Not yet implemented for xN_dependent");
373  } else {
374  throw ValueError(format("xN dependency flag invalid"));
375  }
376 }
377 CoolPropDbl GERG2008ReducingFunction::d2Yrdxidgamma(const std::vector<CoolPropDbl>& x, std::size_t i, const STLMatrix& beta, const STLMatrix& gamma,
378  const STLMatrix& Y_c_ij, const std::vector<CoolPropDbl>& Yc, x_N_dependency_flag xN_flag) const {
379  if (xN_flag == XN_INDEPENDENT) {
380  // See Table B9 from Kunz Wagner 2012 (GERG 2008)
381  CoolPropDbl deriv = 0;
382  for (std::size_t k = 0; k < i; k++) {
383  deriv += 2 * beta[k][i] * Y_c_ij[k][i] * dfYkidxi__constxk(x, k, i, beta);
384  }
385  for (std::size_t k = i + 1; k < N; k++) {
386  deriv += 2 * beta[i][k] * Y_c_ij[i][k] * dfYikdxi__constxk(x, i, k, beta);
387  }
388  return deriv;
389  } else if (xN_flag == XN_DEPENDENT) {
390  // Table S1 from Gernert, 2014, supplemental information
391  if (i == N - 1) {
392  return 0.0;
393  }
394  CoolPropDbl deriv = 0;
395  for (std::size_t k = 0; k < i; k++) {
396  deriv += 2 * beta[k][i] * Y_c_ij[k][i] * dfYkidxi__constxk(x, k, i, beta);
397  }
398  for (std::size_t k = i + 1; k < N - 1; k++) {
399  deriv += 2 * beta[i][k] * Y_c_ij[i][k] * dfYikdxi__constxk(x, i, k, beta);
400  }
401  double beta_Y_iN = beta[i][N - 1], xN = x[N - 1];
402  deriv += 2 * beta[i][N - 1] * Y_c_ij[i][N - 1]
403  * (xN * (x[i] + xN) / (pow(beta_Y_iN, 2) * x[i] + xN)
404  + (1 - beta_Y_iN * beta_Y_iN) * x[i] * xN * xN / POW2(beta_Y_iN * beta_Y_iN * x[i] + xN));
405  for (std::size_t k = 0; k < N - 1; ++k) {
406  double beta_Y_kN = beta[k][N - 1], xk = x[k], beta_Y_kN_squared = beta_Y_kN * beta_Y_kN;
407  deriv += 2 * beta[k][N - 1] * Y_c_ij[k][N - 1]
408  * (-xk * (xk + xN) / (beta_Y_kN_squared * xk + xN) + (1 - beta_Y_kN_squared) * xN * xk * xk / POW2(beta_Y_kN_squared * xk + xN));
409  }
410  return deriv;
411  } else {
412  throw ValueError(format("xN dependency flag invalid"));
413  }
414 }
415 CoolPropDbl GERG2008ReducingFunction::d2Yrdxi2__constxj(const std::vector<CoolPropDbl>& x, std::size_t i, const STLMatrix& beta,
416  const STLMatrix& gamma, const STLMatrix& Y_c_ij, const std::vector<CoolPropDbl>& Yc,
417  x_N_dependency_flag xN_flag) const {
418  if (xN_flag == XN_INDEPENDENT) {
419  // See Table B9 from Kunz Wagner 2012 (GERG 2008)
420  CoolPropDbl d2Yr_dxi2 = 2 * Yc[i];
421  for (std::size_t k = 0; k < i; k++) {
422  d2Yr_dxi2 += c_Y_ij(k, i, beta, gamma, Y_c_ij) * d2fYkidxi2__constxk(x, k, i, beta);
423  }
424  for (std::size_t k = i + 1; k < N; k++) {
425  d2Yr_dxi2 += c_Y_ij(i, k, beta, gamma, Y_c_ij) * d2fYikdxi2__constxk(x, i, k, beta);
426  }
427  return d2Yr_dxi2;
428  } else if (xN_flag == XN_DEPENDENT) {
429  // Table S1 from Gernert, 2014, supplemental information
430  if (i == N - 1) {
431  return 0.0;
432  }
433  CoolPropDbl d2Yr_dxi2 = 2 * Yc[i] + 2 * Yc[N - 1];
434  for (std::size_t k = 0; k < i; k++) {
435  d2Yr_dxi2 += c_Y_ij(k, i, beta, gamma, Y_c_ij) * d2fYkidxi2__constxk(x, k, i, beta);
436  }
437  for (std::size_t k = i + 1; k < N - 1; k++) {
438  d2Yr_dxi2 += c_Y_ij(i, k, beta, gamma, Y_c_ij) * d2fYikdxi2__constxk(x, i, k, beta);
439  }
440  double beta_Y_iN = beta[i][N - 1], xN = x[N - 1];
441  d2Yr_dxi2 += 2 * c_Y_ij(i, N - 1, beta, gamma, Y_c_ij)
442  * (-(x[i] + xN) / (pow(beta_Y_iN, 2) * x[i] + xN)
443  + (1 - beta_Y_iN * beta_Y_iN)
444  * (xN * xN / pow(beta_Y_iN * beta_Y_iN * x[i] + xN, 2)
445  + ((1 - beta_Y_iN * beta_Y_iN) * x[i] * xN * xN - beta_Y_iN * beta_Y_iN * x[i] * x[i] * xN)
446  / pow(beta_Y_iN * beta_Y_iN * x[i] + xN, 3)));
447  for (std::size_t k = 0; k < N - 1; ++k) {
448  double beta_Y_kN = beta[k][N - 1], xk = x[k], beta_Y_kN_squared = beta_Y_kN * beta_Y_kN;
449  d2Yr_dxi2 += 2 * c_Y_ij(k, N - 1, beta, gamma, Y_c_ij) * xk * xk * (1 - beta_Y_kN_squared) / pow(beta_Y_kN_squared * xk + xN, 2)
450  * (xN / (beta_Y_kN_squared * xk + xN) - 1);
451  }
452  return d2Yr_dxi2;
453  } else {
454  throw ValueError(format("xN dependency flag invalid"));
455  }
456 }
457 
458 CoolPropDbl GERG2008ReducingFunction::d2Yrdxidxj(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, const STLMatrix& beta,
459  const STLMatrix& gamma, const STLMatrix& Y_c_ij, const std::vector<CoolPropDbl>& Yc,
460  x_N_dependency_flag xN_flag) const {
461  if (xN_flag == XN_INDEPENDENT) {
462  if (i == j) {
463  return d2Yrdxi2__constxj(x, i, beta, gamma, Y_c_ij, Yc, xN_flag);
464  } else {
465  // See Table B9 from Kunz Wagner 2012 (GERG 2008)
466  return c_Y_ij(i, j, beta, gamma, Y_c_ij) * d2fYijdxidxj(x, i, j, beta);
467  }
468  } else if (xN_flag == XN_DEPENDENT) {
469  // Table S1 from Gernert, 2014, supplemental information
470  if (j == N - 1 || i == N - 1) {
471  return 0.0;
472  }
473  if (i == j) {
474  return d2Yrdxi2__constxj(x, i, beta, gamma, Y_c_ij, Yc, xN_flag);
475  }
476  CoolPropDbl d2Yr_dxidxj = 2 * Yc[N - 1];
477  d2Yr_dxidxj += c_Y_ij(i, j, beta, gamma, Y_c_ij) * d2fYijdxidxj(x, i, j, beta);
478 
479  for (std::size_t k = 0; k < N - 1; k++) {
480  d2Yr_dxidxj += c_Y_ij(k, N - 1, beta, gamma, Y_c_ij) * d2fYkidxi2__constxk(x, k, N - 1, beta);
481  }
482  d2Yr_dxidxj -= c_Y_ij(i, N - 1, beta, gamma, Y_c_ij) * d2fYijdxidxj(x, i, N - 1, beta);
483  d2Yr_dxidxj -= c_Y_ij(j, N - 1, beta, gamma, Y_c_ij) * d2fYijdxidxj(x, j, N - 1, beta);
484  return d2Yr_dxidxj;
485  } else {
486  throw ValueError(format("xN dependency flag invalid"));
487  }
488 }
489 
490 CoolPropDbl GERG2008ReducingFunction::d3Yrdxidxjdxk(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, std::size_t k,
491  const STLMatrix& beta, const STLMatrix& gamma, const STLMatrix& Y_c_ij,
492  const std::vector<CoolPropDbl>& Yc, x_N_dependency_flag xN_flag) const {
493  if (xN_flag == XN_INDEPENDENT) {
494  if (i != j && j != k && k != i) {
495  return 0;
496  } else if (k == i && i != j) {
497  return c_Y_ij(i, j, beta, gamma, Y_c_ij) * d3fYijdxi2dxj(x, i, j, beta);
498  } else if (k == j && i != j) {
499  return c_Y_ij(i, j, beta, gamma, Y_c_ij) * d3fYijdxidxj2(x, i, j, beta);
500  } else if (i == j && i != k) {
501  return c_Y_ij(i, k, beta, gamma, Y_c_ij) * d3fYijdxi2dxj(x, i, k, beta);
502  } else {
503  CoolPropDbl d3Yr_dxi3 = 0;
504  for (std::size_t m = 0; m < i; m++) {
505  d3Yr_dxi3 += c_Y_ij(m, i, beta, gamma, Y_c_ij) * d3fYkidxi3__constxk(x, m, i, beta);
506  }
507  for (std::size_t m = i + 1; m < N; m++) {
508  d3Yr_dxi3 += c_Y_ij(i, m, beta, gamma, Y_c_ij) * d3fYikdxi3__constxk(x, i, m, beta);
509  }
510  return d3Yr_dxi3;
511  }
512  } else if (xN_flag == XN_DEPENDENT) {
513  CoolPropDbl summer = 0;
514  // Needed for all third partials
515  for (std::size_t m = 0; m < N - 1; m++) {
516  summer -= c_Y_ij(m, N - 1, beta, gamma, Y_c_ij) * d3fYkidxi3__constxk(x, m, N - 1, beta);
517  }
518  if (i != j && j != k && k != i) {
519  summer += c_Y_ij(i, N - 1, beta, gamma, Y_c_ij) * d3fYijdxidxj2(x, i, N - 1, beta);
520  summer += c_Y_ij(j, N - 1, beta, gamma, Y_c_ij) * d3fYijdxidxj2(x, j, N - 1, beta);
521  summer += c_Y_ij(k, N - 1, beta, gamma, Y_c_ij) * d3fYijdxidxj2(x, k, N - 1, beta);
522  } else if (k == i && i != j) { // two i, one j
523  summer += c_Y_ij(i, j, beta, gamma, Y_c_ij) * d3fYijdxi2dxj(x, i, j, beta);
524  summer += c_Y_ij(j, N - 1, beta, gamma, Y_c_ij) * d3fYijdxidxj2(x, j, N - 1, beta);
525  summer += c_Y_ij(i, N - 1, beta, gamma, Y_c_ij) * (2 * d3fYijdxidxj2(x, i, N - 1, beta) - d3fYijdxi2dxj(x, i, N - 1, beta));
526  } else if (k == j && i != j) { // two j, one i
527  summer += c_Y_ij(i, j, beta, gamma, Y_c_ij) * d3fYijdxidxj2(x, i, j, beta);
528  summer += c_Y_ij(i, N - 1, beta, gamma, Y_c_ij) * d3fYijdxidxj2(x, i, N - 1, beta);
529  summer += c_Y_ij(j, N - 1, beta, gamma, Y_c_ij) * (2 * d3fYijdxidxj2(x, j, N - 1, beta) - d3fYijdxi2dxj(x, j, N - 1, beta));
530  } else if (i == j && i != k) { // two i, one k
531  summer += c_Y_ij(i, k, beta, gamma, Y_c_ij) * d3fYijdxi2dxj(x, i, k, beta);
532  summer += c_Y_ij(k, N - 1, beta, gamma, Y_c_ij) * d3fYijdxidxj2(x, k, N - 1, beta);
533  summer += c_Y_ij(i, N - 1, beta, gamma, Y_c_ij) * (2 * d3fYijdxidxj2(x, i, N - 1, beta) - d3fYijdxi2dxj(x, i, N - 1, beta));
534  } else {
535  for (std::size_t m = 0; m < i; m++) {
536  summer += c_Y_ij(m, i, beta, gamma, Y_c_ij) * d3fYkidxi3__constxk(x, m, i, beta);
537  }
538  for (std::size_t m = i + 1; m < N - 1; m++) {
539  summer += c_Y_ij(i, m, beta, gamma, Y_c_ij) * d3fYikdxi3__constxk(x, i, m, beta);
540  }
541  summer += c_Y_ij(i, N - 1, beta, gamma, Y_c_ij)
542  * (3 * d3fYijdxidxj2(x, i, N - 1, beta) - 3 * d3fYijdxi2dxj(x, i, N - 1, beta) + d3fYikdxi3__constxk(x, i, N - 1, beta));
543  }
544  return summer;
545  } else {
546  throw ValueError(format("xN dependency flag invalid"));
547  }
548 }
549 
550 CoolPropDbl GERG2008ReducingFunction::dfYkidxi__constxk(const std::vector<CoolPropDbl>& x, std::size_t k, std::size_t i,
551  const STLMatrix& beta) const {
552  double xk = x[k], xi = x[i], beta_Y = beta[k][i], beta_Y_squared = beta_Y * beta_Y;
553  return xk * (xk + xi) / (beta_Y_squared * xk + xi) + xk * xi / (beta_Y_squared * xk + xi) * (1 - (xk + xi) / (beta_Y_squared * xk + xi));
554 }
555 CoolPropDbl GERG2008ReducingFunction::dfYikdxi__constxk(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t k,
556  const STLMatrix& beta) const {
557  double xk = x[k], xi = x[i], beta_Y = beta[i][k], beta_Y_squared = beta_Y * beta_Y;
558  return xk * (xi + xk) / (beta_Y_squared * xi + xk)
559  + xi * xk / (beta_Y_squared * xi + xk) * (1 - beta_Y_squared * (xi + xk) / (beta_Y_squared * xi + xk));
560 }
561 const CoolPropDbl GERG2008ReducingFunction::c_Y_ij(const std::size_t i, const std::size_t j, const STLMatrix& beta, const STLMatrix& gamma,
562  const STLMatrix& Y_c) const {
563  return 2 * beta[i][j] * gamma[i][j] * Y_c[i][j];
564 }
565 CoolPropDbl GERG2008ReducingFunction::f_Y_ij(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, const STLMatrix& beta) const {
566  double xi = x[i], xj = x[j], beta_Y = beta[i][j];
567  return xi * xj * (xi + xj) / (beta_Y * beta_Y * xi + xj);
568 }
569 CoolPropDbl GERG2008ReducingFunction::d2fYikdxi2__constxk(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t k,
570  const STLMatrix& beta) const {
571  double xi = x[i], xk = x[k], beta_Y = beta[i][k];
572  return 1 / (beta_Y * beta_Y * xi + xk) * (1 - beta_Y * beta_Y * (xi + xk) / (beta_Y * beta_Y * xi + xk))
573  * (2 * xk - xi * xk * 2 * beta_Y * beta_Y / (beta_Y * beta_Y * xi + xk));
574 }
575 CoolPropDbl GERG2008ReducingFunction::d2fYkidxi2__constxk(const std::vector<CoolPropDbl>& x, std::size_t k, std::size_t i,
576  const STLMatrix& beta) const {
577  double xi = x[i], xk = x[k], beta_Y = beta[k][i];
578  return 1 / (beta_Y * beta_Y * xk + xi) * (1 - (xk + xi) / (beta_Y * beta_Y * xk + xi)) * (2 * xk - xk * xi * 2 / (beta_Y * beta_Y * xk + xi));
579 }
580 CoolPropDbl GERG2008ReducingFunction::d2fYijdxidxj(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, const STLMatrix& beta) const {
581  double xi = x[i], xj = x[j], beta_Y = beta[i][j], beta_Y2 = beta_Y * beta_Y;
582  return (xi + xj) / (beta_Y2 * xi + xj) + xj / (beta_Y2 * xi + xj) * (1 - (xi + xj) / (beta_Y2 * xi + xj))
583  + xi / (beta_Y2 * xi + xj) * (1 - beta_Y2 * (xi + xj) / (beta_Y2 * xi + xj))
584  - xi * xj / pow(beta_Y2 * xi + xj, 2) * (1 + beta_Y2 - 2 * beta_Y2 * (xi + xj) / (beta_Y2 * xi + xj));
585 }
586 CoolPropDbl GERG2008ReducingFunction::d3fYijdxi2dxj(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, const STLMatrix& beta) const {
587  double x_i = x[i], x_j = x[j], beta_Y = beta[i][j], beta_Y2 = beta_Y * beta_Y;
588  double den = pow(beta_Y, 8) * pow(x_i, 4) + 4 * pow(beta_Y, 6) * pow(x_i, 3) * x_j + 6 * pow(beta_Y, 4) * pow(x_i * x_j, 2)
589  + 4 * beta_Y2 * x_i * pow(x_j, 3) + pow(x_j, 4);
590  return -6 * beta_Y2 * x_i * x_j * x_j * (beta_Y2 - 1) / den;
591 }
592 CoolPropDbl GERG2008ReducingFunction::d3fYijdxidxj2(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t j, const STLMatrix& beta) const {
593  double x_i = x[i], x_j = x[j], beta_Y = beta[i][j], beta_Y2 = beta_Y * beta_Y;
594  double den = pow(beta_Y, 8) * pow(x_i, 4) + 4 * pow(beta_Y, 6) * pow(x_i, 3) * x_j + 6 * pow(beta_Y, 4) * pow(x_i * x_j, 2)
595  + 4 * beta_Y2 * x_i * pow(x_j, 3) + pow(x_j, 4);
596  return 6 * beta_Y2 * x_i * x_i * x_j * (beta_Y2 - 1) / den;
597 }
598 CoolPropDbl GERG2008ReducingFunction::d3fYikdxi3__constxk(const std::vector<CoolPropDbl>& x, std::size_t i, std::size_t k,
599  const STLMatrix& beta) const {
600  double x_i = x[i], x_k = x[k], beta_Y = beta[i][k], beta_Y2 = beta_Y * beta_Y;
601  double den = pow(beta_Y, 8) * pow(x_i, 4) + 4 * pow(beta_Y, 6) * pow(x_i, 3) * x_k + 6 * pow(beta_Y, 4) * pow(x_i * x_k, 2)
602  + 4 * beta_Y2 * x_i * pow(x_k, 3) + pow(x_k, 4);
603  return 6 * beta_Y2 * x_k * x_k * x_k * (beta_Y2 - 1) / den;
604 }
605 CoolPropDbl GERG2008ReducingFunction::d3fYkidxi3__constxk(const std::vector<CoolPropDbl>& x, std::size_t k, std::size_t i,
606  const STLMatrix& beta) const {
607  double x_i = x[i], x_k = x[k], beta_Y = beta[k][i], beta_Y2 = beta_Y * beta_Y;
608  return 6 * beta_Y2 * x_k * x_k * x_k * (1 - beta_Y2) / pow(beta_Y2 * x_k + x_i, 4);
609 }
610 
611 } /* namespace CoolProp */
Generated on Wed Nov 29 2023 14:26:24 for CoolProp by doxygen 1.9.1