CoolProp 6.8.1dev
An open-source fluid property and humid air property database
ReducingFunctions.cpp
Go to the documentation of this file.
1#include "ReducingFunctions.h"
2
3namespace CoolProp {
4
5CoolPropDbl 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}
30CoolPropDbl 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}
55CoolPropDbl 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}
63CoolPropDbl 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}
71CoolPropDbl 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}
91CoolPropDbl 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}
112CoolPropDbl 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}
115CoolPropDbl 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}
118CoolPropDbl 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}
128CoolPropDbl 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}
131CoolPropDbl 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}
134CoolPropDbl 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
144CoolPropDbl GERG2008ReducingFunction::Tr(const std::vector<CoolPropDbl>& x) const {
145 return Yr(x, beta_T, gamma_T, T_c, Yc_T);
146}
147CoolPropDbl GERG2008ReducingFunction::dTr_dbetaT(const std::vector<CoolPropDbl>& x) const {
148 return dYr_dbeta(x, beta_T, gamma_T, T_c, Yc_T);
149}
150CoolPropDbl GERG2008ReducingFunction::dTr_dgammaT(const std::vector<CoolPropDbl>& x) const {
151 return dYr_dgamma(x, beta_T, gamma_T, T_c, Yc_T);
152}
153CoolPropDbl 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};
156CoolPropDbl 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
160CoolPropDbl 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}
163CoolPropDbl 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}
166CoolPropDbl 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}
169CoolPropDbl 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}
173CoolPropDbl GERG2008ReducingFunction::rhormolar(const std::vector<CoolPropDbl>& x) const {
174 return 1 / Yr(x, beta_v, gamma_v, v_c, Yc_v);
175}
176
177CoolPropDbl 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}
182CoolPropDbl 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}
187CoolPropDbl 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};
190CoolPropDbl 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};
193CoolPropDbl 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}
197CoolPropDbl 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}
201CoolPropDbl 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}
205CoolPropDbl 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}
208CoolPropDbl 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}
211CoolPropDbl 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}
215CoolPropDbl 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}
219CoolPropDbl 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}
224CoolPropDbl 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}
231CoolPropDbl 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
244CoolPropDbl 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
263CoolPropDbl 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}
277CoolPropDbl 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
298CoolPropDbl 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}
339CoolPropDbl 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}
377CoolPropDbl 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}
415CoolPropDbl 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
458CoolPropDbl 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
490CoolPropDbl 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
550CoolPropDbl 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}
555CoolPropDbl 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}
561const 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}
565CoolPropDbl 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}
569CoolPropDbl 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}
575CoolPropDbl 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}
580CoolPropDbl 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}
586CoolPropDbl 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}
592CoolPropDbl 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}
598CoolPropDbl 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}
605CoolPropDbl 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 Jun 4 2025 04:12:26 for CoolProp by doxygen 1.9.4