CoolProp 6.8.0
An open-source fluid property and humid air property database
MixtureDerivatives.cpp
Go to the documentation of this file.
3
4namespace CoolProp {
5
7 return dln_fugacity_coefficient_dT__constp_n(HEOS, i, xN_flag);
8}
10 return dln_fugacity_coefficient_dp__constT_n(HEOS, i, xN_flag) + 1 / HEOS.p();
11}
13 return HEOS.mole_fractions[i] * HEOS.rhomolar() * HEOS.gas_constant() * HEOS.T() * exp(dnalphar_dni__constT_V_nj(HEOS, i, xN_flag));
14}
16 return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i, xN_flag) - log(1 + HEOS._delta.pt() * HEOS.dalphar_dDelta());
17}
19 return 1 / HEOS.T() * (1 - HEOS.tau() * HEOS.dalphar_dTau() - HEOS.tau() * d_ndalphardni_dTau(HEOS, i, xN_flag));
20}
22 return 1 / HEOS.rhomolar() * (1 + HEOS.delta() * HEOS.dalphar_dDelta() + HEOS.delta() * d_ndalphardni_dDelta(HEOS, i, xN_flag));
23}
25 // GERG Equation 7.42
26 return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i, xN_flag);
27}
29 return -HEOS._tau.pt() / HEOS._T * (HEOS.dalphar_dTau() + d_ndalphardni_dTau(HEOS, i, xN_flag));
30}
32 double T = HEOS._reducing.T / HEOS._tau.pt();
33 CoolPropDbl R_u = HEOS.gas_constant();
34 return d2nalphar_dni_dT(HEOS, i, xN_flag) + 1 / T - partial_molar_volume(HEOS, i, xN_flag) / (R_u * T) * dpdT__constV_n(HEOS);
35}
37 return -ndpdni__constT_V_nj(HEOS, i, xN_flag) / ndpdV__constT_n(HEOS);
38}
39
41 // GERG equation 7.30
42 CoolPropDbl R_u = HEOS.gas_constant();
43 double partial_molar_volumeval = partial_molar_volume(HEOS, i, xN_flag); // [m^3/mol]
44 double term1 = partial_molar_volumeval / (R_u * HEOS._T); // m^3/mol/(N*m)*mol = m^2/N = 1/Pa
45 double term2 = 1.0 / HEOS.p();
46 return term1 - term2;
47}
49 return -1 / HEOS.tau() + HEOS.dalphar_dTau() + d_ndalphardni_dTau(HEOS, i, xN_flag);
50}
52 return 1 + HEOS.delta() * HEOS.dalphar_dDelta() + HEOS.delta() * d_ndalphardni_dDelta(HEOS, i, xN_flag);
53}
55 x_N_dependency_flag xN_flag) {
56 // This is a term to which some more might be added depending on i and j
58 const std::vector<CoolPropDbl>& x = HEOS.get_mole_fractions();
59 std::size_t N = x.size();
60 if (i == N - 1) {
61 val += -1 / x[N - 1];
62 } else if (i == j) {
63 val += 1 / x[j];
64 }
65 return val;
66}
68 x_N_dependency_flag xN_flag) {
69 if (xN_flag == XN_INDEPENDENT) {
70 throw ValueError("dln_fugacity_dxj__constT_rho_xi only valid for xN_DEPENDENT for now");
71 }
72 CoolPropDbl rhor = HEOS.Reducing->rhormolar(HEOS.get_mole_fractions());
73 CoolPropDbl Tr = HEOS.Reducing->Tr(HEOS.get_mole_fractions());
74 CoolPropDbl dTrdxj = HEOS.Reducing->dTrdxi__constxj(HEOS.get_mole_fractions(), j, xN_flag);
75 CoolPropDbl drhordxj = HEOS.Reducing->drhormolardxi__constxj(HEOS.get_mole_fractions(), j, xN_flag);
76
77 // These lines are all the same
78 CoolPropDbl line1 = dln_fugacity_i_dtau__constdelta_x(HEOS, i, xN_flag) * 1 / HEOS.T() * dTrdxj;
79 CoolPropDbl line2 = -dln_fugacity_i_ddelta__consttau_x(HEOS, i, xN_flag) * 1 / rhor * drhordxj;
80 CoolPropDbl line4 = HEOS.residual_helmholtz->dalphar_dxi(HEOS, j, xN_flag) + d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, j, xN_flag);
81
82 const std::vector<CoolPropDbl>& x = HEOS.get_mole_fractions();
83 std::size_t N = x.size();
84
85 CoolPropDbl line3 = 1 / rhor * HEOS.Reducing->drhormolardxi__constxj(x, j, xN_flag) + 1 / Tr * HEOS.Reducing->dTrdxi__constxj(x, j, xN_flag);
86 ;
87
88 // This is a term to which some more might be added depending on i and j
89 if (i == N - 1) {
90 line3 += -1 / x[N - 1];
91 } else if (i == j) {
92 line3 += 1 / x[j];
93 } else {
94 }
95
96 return line1 + line2 + line3 + line4;
97}
98
100 x_N_dependency_flag xN_flag) {
101 double s = (HEOS.mole_fractions[i] > DBL_EPSILON) ? Kronecker_delta(i, j) / HEOS.mole_fractions[i] : 0;
102 return s + nd2nalphardnidnj__constT_V(HEOS, i, j, xN_flag);
103}
105 x_N_dependency_flag xN_flag) {
106 return d_ndalphardni_dTau(HEOS, j, xN_flag) + d_nd_ndalphardni_dnj_dTau__constdelta_x(HEOS, i, j, xN_flag);
107}
109 x_N_dependency_flag xN_flag) {
110 return d2_ndalphardni_dTau2(HEOS, j, xN_flag) + d2_nd_ndalphardni_dnj_dTau2__constdelta_x(HEOS, i, j, xN_flag);
111}
113 x_N_dependency_flag xN_flag) {
114 return d2_ndalphardni_dDelta2(HEOS, j, xN_flag) + d2_nd_ndalphardni_dnj_dDelta2__consttau_x(HEOS, i, j, xN_flag);
115}
117 x_N_dependency_flag xN_flag) {
118 return d2_ndalphardni_dDelta_dTau(HEOS, j, xN_flag) + d2_nd_ndalphardni_dnj_dDelta_dTau__constx(HEOS, i, j, xN_flag);
119}
120
122 x_N_dependency_flag xN_flag) {
123 return d_ndalphardni_dDelta(HEOS, j, xN_flag) + d_nd_ndalphardni_dnj_dDelta__consttau_x(HEOS, i, j, xN_flag);
124}
126 std::size_t k, x_N_dependency_flag xN_flag) {
127 CoolPropDbl s = (HEOS.mole_fractions[i] > DBL_EPSILON) ? -Kronecker_delta(i, j) * Kronecker_delta(i, k) / pow(HEOS.mole_fractions[i], 2) : 0;
128 return s + d_ndalphardni_dxj__constdelta_tau_xi(HEOS, j, k, xN_flag) + d_nd_ndalphardni_dnj_dxk__consttau_delta(HEOS, i, j, k, xN_flag);
129}
131 std::size_t k, x_N_dependency_flag xN_flag) {
132 return d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, j, k, xN_flag) + d2_nd_ndalphardni_dnj_dxk_dTau__constdelta(HEOS, i, j, k, xN_flag);
133}
135 std::size_t k, x_N_dependency_flag xN_flag) {
136 return d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, j, k, xN_flag) + d2_nd_ndalphardni_dnj_dxk_dDelta__consttau(HEOS, i, j, k, xN_flag);
137}
139 x_N_dependency_flag xN_flag) {
140 double sum = d_ndln_fugacity_i_dnj_dtau__constdelta_x(HEOS, i, j, xN_flag) * ndtaudni__constT_V_nj(HEOS, k, xN_flag)
141 + d_ndln_fugacity_i_dnj_ddelta__consttau_x(HEOS, i, j, xN_flag) * nddeltadni__constT_V_nj(HEOS, k, xN_flag)
142 + d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HEOS, i, j, k, xN_flag);
143 std::size_t mmax = HEOS.mole_fractions.size();
144 if (xN_flag == XN_DEPENDENT) {
145 mmax--;
146 }
147 for (unsigned int m = 0; m < mmax; ++m) {
148 sum -= HEOS.mole_fractions[m] * d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HEOS, i, j, m, xN_flag);
149 }
150 return sum;
151}
152
154 x_N_dependency_flag xN_flag) {
155 // Gernert 3.115
156 CoolPropDbl R_u = HEOS.gas_constant();
157 // partial molar volume is -dpdn/dpdV, so need to flip the sign here
158 return d2nalphar_dxj_dni__constT_V(HEOS, j, i, xN_flag)
159 - partial_molar_volume(HEOS, i, xN_flag) / (R_u * HEOS._T) * dpdxj__constT_V_xi(HEOS, j, xN_flag);
160}
162 // Gernert 3.130
163 CoolPropDbl R_u = HEOS.gas_constant();
164 return HEOS._rhomolar * R_u * HEOS._T
165 * (ddelta_dxj__constT_V_xi(HEOS, j, xN_flag) * HEOS.dalphar_dDelta()
166 + HEOS._delta.pt() * d_dalpharddelta_dxj__constT_V_xi(HEOS, j, xN_flag));
167}
168
170 // Gernert Equation 3.134 (Catch test provided)
171 return HEOS.d2alphar_dDelta2() * ddelta_dxj__constT_V_xi(HEOS, j, xN_flag) + HEOS.d2alphar_dDelta_dTau() * dtau_dxj__constT_V_xi(HEOS, j, xN_flag)
172 + HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, j, xN_flag);
173}
174
176 //Gernert 3.119 (Catch test provided)
177 return HEOS.dalphar_dDelta() * ddelta_dxj__constT_V_xi(HEOS, j, xN_flag) + HEOS.dalphar_dTau() * dtau_dxj__constT_V_xi(HEOS, j, xN_flag)
178 + HEOS.residual_helmholtz->dalphar_dxi(HEOS, j, xN_flag);
179}
181 x_N_dependency_flag xN_flag) {
182 // Gernert 3.118
183 return d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, j, xN_flag)
184 + ddelta_dxj__constT_V_xi(HEOS, j, xN_flag) * d_ndalphardni_dDelta(HEOS, i, xN_flag)
185 + dtau_dxj__constT_V_xi(HEOS, j, xN_flag) * d_ndalphardni_dTau(HEOS, i, xN_flag);
186}
188 // Gernert 3.121 (Catch test provided)
189 return -HEOS._delta.pt() / HEOS._reducing.rhomolar * HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, j, xN_flag);
190}
192 // Gernert 3.122 (Catch test provided)
193 return 1 / HEOS._T * HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag);
194}
195
197 CoolPropDbl R_u = HEOS.gas_constant();
198 return HEOS._rhomolar * R_u * (1 + HEOS._delta.pt() * HEOS.dalphar_dDelta() - HEOS._delta.pt() * HEOS._tau.pt() * HEOS.d2alphar_dDelta_dTau());
199}
201 CoolPropDbl R_u = HEOS.gas_constant();
202 return R_u * HEOS._T * (1 + 2 * HEOS._delta.pt() * HEOS.dalphar_dDelta() + pow(HEOS._delta.pt(), 2) * HEOS.d2alphar_dDelta2());
203}
205 CoolPropDbl R_u = HEOS.gas_constant();
206 return -pow(HEOS._rhomolar, 2) * R_u * HEOS._T
207 * (1 + 2 * HEOS._delta.pt() * HEOS.dalphar_dDelta() + pow(HEOS._delta.pt(), 2) * HEOS.d2alphar_dDelta2());
208}
210 // Eqn 7.64 and 7.63
211 CoolPropDbl R_u = HEOS.gas_constant();
212 double ndrhorbar_dni__constnj = HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag);
213 double ndTr_dni__constnj = HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
214 double summer = 0;
215 std::size_t kmax = HEOS.mole_fractions.size();
216 if (xN_flag == XN_DEPENDENT) {
217 kmax--;
218 }
219 for (unsigned int k = 0; k < kmax; ++k) {
220 summer += HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, k, xN_flag);
221 }
222 double nd2alphar_dni_dDelta = HEOS._delta.pt() * HEOS.d2alphar_dDelta2() * (1 - 1 / HEOS._reducing.rhomolar * ndrhorbar_dni__constnj)
223 + HEOS._tau.pt() * HEOS.d2alphar_dDelta_dTau() / HEOS._reducing.T * ndTr_dni__constnj
224 + HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, i, xN_flag) - summer;
225 return HEOS._rhomolar * R_u * HEOS._T
226 * (1 + HEOS._delta.pt() * HEOS.dalphar_dDelta() * (2 - 1 / HEOS._reducing.rhomolar * ndrhorbar_dni__constnj)
227 + HEOS._delta.pt() * nd2alphar_dni_dDelta);
228}
229
231 double term1 = HEOS._delta.pt() * HEOS.dalphar_dDelta()
232 * (1 - 1 / HEOS._reducing.rhomolar * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
233 double term2 = HEOS._tau.pt() * HEOS.dalphar_dTau() * (1 / HEOS._reducing.T) * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
234
235 double s = 0;
236 std::size_t kmax = HEOS.mole_fractions.size();
237 if (xN_flag == XN_DEPENDENT) {
238 kmax--;
239 }
240 for (unsigned int k = 0; k < kmax; k++) {
241 s += HEOS.mole_fractions[k] * HEOS.residual_helmholtz->dalphar_dxi(HEOS, k, xN_flag);
242 }
243 double term3 = HEOS.residual_helmholtz->dalphar_dxi(HEOS, i, xN_flag);
244 return term1 + term2 + term3 - s;
245}
247 x_N_dependency_flag xN_flag) {
248 CoolPropDbl R_u = HEOS.gas_constant();
249 return nd2nalphardnidnj__constT_V(HEOS, j, i, xN_flag) + 1
250 - partial_molar_volume(HEOS, j, xN_flag) / (R_u * HEOS._T) * ndpdni__constT_V_nj(HEOS, i, xN_flag);
251}
253 return HEOS._delta.pt() - HEOS._delta.pt() / HEOS._reducing.rhomolar * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag);
254}
256 return 1 - 1 / HEOS._reducing.rhomolar * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag);
257}
259 x_N_dependency_flag xN_flag) {
260 double rhor = HEOS._reducing.rhomolar;
261 return -HEOS.delta() / rhor
262 * (HEOS.Reducing->d_ndrhorbardni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag)
263 - 1 / rhor * HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, j, xN_flag)
264 * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
265}
267 x_N_dependency_flag xN_flag) {
268 return d_nddeltadni_dxj__constdelta_tau(HEOS, i, j, xN_flag) / HEOS.delta();
269}
270
272 return HEOS._tau.pt() / HEOS._reducing.T * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
273}
274
276 return 1 / HEOS._reducing.T * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
277}
279 x_N_dependency_flag xN_flag) {
280 double Tr = HEOS._reducing.T;
281 return HEOS.tau() / Tr
282 * (HEOS.Reducing->d_ndTrdni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag)
283 - 1 / Tr * HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag)
284 * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag));
285}
287 x_N_dependency_flag xN_flag) {
288 return d_ndtaudni_dxj__constdelta_tau(HEOS, i, j, xN_flag) / HEOS.tau();
289}
291 x_N_dependency_flag xN_flag) {
292 double line1 = HEOS._delta.pt() * HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, j, xN_flag)
293 * (1 - 1 / HEOS._reducing.rhomolar * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
294 double line3 = HEOS._tau.pt() * HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, j, xN_flag) * (1 / HEOS._reducing.T)
295 * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
296 double line2 = -HEOS._delta.pt() * HEOS.dalphar_dDelta() * (1 / HEOS._reducing.rhomolar)
297 * (HEOS.Reducing->d_ndrhorbardni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag)
298 - 1 / HEOS._reducing.rhomolar * HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, j, xN_flag)
299 * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
300 double line4 = HEOS._tau.pt() * HEOS.dalphar_dTau() * (1 / HEOS._reducing.T)
301 * (HEOS.Reducing->d_ndTrdni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag)
302 - 1 / HEOS._reducing.T * HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag)
303 * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag));
304
305 double s = 0;
306 std::size_t kmax = HEOS.mole_fractions.size();
307 if (xN_flag == XN_DEPENDENT) {
308 kmax--;
309 }
310 for (unsigned int k = 0; k < kmax; k++) {
311 s += HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d2alphardxidxj(HEOS, j, k, xN_flag);
312 }
313 double line5 = HEOS.residual_helmholtz->d2alphardxidxj(HEOS, i, j, xN_flag) - HEOS.residual_helmholtz->dalphar_dxi(HEOS, j, xN_flag) - s;
314 return line1 + line2 + line3 + line4 + line5;
315}
316
318 x_N_dependency_flag xN_flag) {
319 return ndalphar_dni__constT_V_nj(HEOS, j, xN_flag) // First term from 7.46
320 + nd_ndalphardni_dnj__constT_V(HEOS, i, j, xN_flag); // 7.47
321}
322
324 x_N_dependency_flag xN_flag) {
325 double line1 = d_ndalphardni_dDelta(HEOS, i, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
326 double line2 = d_ndalphardni_dTau(HEOS, i, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
327 double line3 = d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, j, xN_flag);
328 std::size_t kmax = HEOS.mole_fractions.size();
329 if (xN_flag == XN_DEPENDENT) {
330 kmax--;
331 }
332 for (unsigned int k = 0; k < kmax; k++) {
333 line3 -= HEOS.mole_fractions[k] * d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, k, xN_flag);
334 }
335 return line1 + line2 + line3;
336}
338 x_N_dependency_flag xN_flag) {
339 double line1 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
340 double line2 = d2_ndalphardni_dTau2(HEOS, i, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag)
341 + d_ndalphardni_dTau(HEOS, i, xN_flag) * d_ndtaudni_dTau(HEOS, j, xN_flag);
342 double summer = 0;
343 std::size_t kmax = HEOS.mole_fractions.size();
344 if (xN_flag == XN_DEPENDENT) {
345 kmax--;
346 }
347 for (unsigned int k = 0; k < kmax; k++) {
348 summer += HEOS.mole_fractions[k] * d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag);
349 }
350 double line3 = d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, j, xN_flag) - summer;
351 return line1 + line2 + line3;
352}
354 x_N_dependency_flag xN_flag) {
355 double line1 = d3_ndalphardni_dDelta_dTau2(HEOS, i, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
356 double line2 = 2 * d2_ndalphardni_dTau2(HEOS, i, xN_flag) * d_ndtaudni_dTau(HEOS, j, xN_flag);
357 double line3 = d3_ndalphardni_dTau3(HEOS, i, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
358 double summer = 0;
359 std::size_t kmax = HEOS.mole_fractions.size();
360 if (xN_flag == XN_DEPENDENT) {
361 kmax--;
362 }
363 for (unsigned int k = 0; k < kmax; k++) {
364 summer += HEOS.mole_fractions[k] * d3_ndalphardni_dxj_dTau2__constdelta_xi(HEOS, i, k, xN_flag);
365 }
366 double line4 = d3_ndalphardni_dxj_dTau2__constdelta_xi(HEOS, i, j, xN_flag) - summer;
367 return line1 + line2 + line3 + line4;
368}
370 x_N_dependency_flag xN_flag) {
371 double line1 = d3_ndalphardni_dDelta3(HEOS, i, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
372 double line2 = 2 * d2_ndalphardni_dDelta2(HEOS, i, xN_flag) * d_nddeltadni_dDelta(HEOS, j, xN_flag);
373 double line3 = d3_ndalphardni_dDelta2_dTau(HEOS, i, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
374 double summer = 0;
375 std::size_t kmax = HEOS.mole_fractions.size();
376 if (xN_flag == XN_DEPENDENT) {
377 kmax--;
378 }
379 for (unsigned int k = 0; k < kmax; k++) {
380 summer += HEOS.mole_fractions[k] * d3_ndalphardni_dxj_dDelta2__consttau_xi(HEOS, i, k, xN_flag);
381 }
382 double line4 = d3_ndalphardni_dxj_dDelta2__consttau_xi(HEOS, i, j, xN_flag) - summer;
383 return line1 + line2 + line3 + line4;
384}
386 x_N_dependency_flag xN_flag) {
387 double line1 = d3_ndalphardni_dDelta2_dTau(HEOS, i, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
388 double line2 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag) * d_nddeltadni_dDelta(HEOS, j, xN_flag);
389 double line3 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag) * d_ndtaudni_dTau(HEOS, j, xN_flag);
390 double line4 = d3_ndalphardni_dDelta_dTau2(HEOS, i, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
391 double summer = 0;
392 std::size_t kmax = HEOS.mole_fractions.size();
393 if (xN_flag == XN_DEPENDENT) {
394 kmax--;
395 }
396 for (unsigned int k = 0; k < kmax; k++) {
397 summer += HEOS.mole_fractions[k] * d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, k, xN_flag);
398 }
399 double line5 = d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, j, xN_flag) - summer;
400 return line1 + line2 + line3 + line4 + line5;
401}
403 x_N_dependency_flag xN_flag) {
404 double line1 = d2_ndalphardni_dDelta2(HEOS, i, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag)
405 + d_ndalphardni_dDelta(HEOS, i, xN_flag) * d_nddeltadni_dDelta(HEOS, j, xN_flag);
406 double line2 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
407 double summer = 0;
408 std::size_t kmax = HEOS.mole_fractions.size();
409 if (xN_flag == XN_DEPENDENT) {
410 kmax--;
411 }
412 for (unsigned int k = 0; k < kmax; k++) {
413 summer += HEOS.mole_fractions[k] * d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, k, xN_flag);
414 }
415 double line3 = d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, j, xN_flag) - summer;
416 return line1 + line2 + line3;
417}
418
420 std::size_t k, x_N_dependency_flag xN_flag) {
421 double line1 = d_ndalphardni_dDelta(HEOS, i, xN_flag) * d_nddeltadni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
422 + d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, k, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
423
424 double line2 = d_ndalphardni_dTau(HEOS, i, xN_flag) * d_ndtaudni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
425 + d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
426
427 double line3 = d2_ndalphardni_dxj_dxk__constdelta_tau_xi(HEOS, i, j, k, xN_flag) - d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, k, xN_flag);
428
429 std::size_t mmax = HEOS.mole_fractions.size();
430 if (xN_flag == XN_DEPENDENT) {
431 mmax--;
432 }
433 for (unsigned int m = 0; m < mmax; m++) {
434 line3 -= HEOS.mole_fractions[m] * d2_ndalphardni_dxj_dxk__constdelta_tau_xi(HEOS, i, m, k, xN_flag);
435 }
436 return line1 + line2 + line3;
437}
439 std::size_t k, x_N_dependency_flag xN_flag) {
440 double line1 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag) * d_nddeltadni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
441 + d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, k, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
442
443 double line2 = d_ndalphardni_dTau(HEOS, i, xN_flag) * d2_ndtaudni_dxj_dTau__constdelta(HEOS, j, k, xN_flag)
444 + d2_ndalphardni_dTau2(HEOS, i, xN_flag) * d_ndtaudni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
445 + d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag) * d_ndtaudni_dTau(HEOS, j, xN_flag)
446 + d3_ndalphardni_dxj_dTau2__constdelta_xi(HEOS, i, k, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
447
448 double line3 = d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(HEOS, i, j, k, xN_flag) - d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag);
449
450 std::size_t mmax = HEOS.mole_fractions.size();
451 if (xN_flag == XN_DEPENDENT) {
452 mmax--;
453 }
454 for (unsigned int m = 0; m < mmax; m++) {
455 line3 -= HEOS.mole_fractions[m] * d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(HEOS, i, m, k, xN_flag);
456 }
457 return line1 + line2 + line3;
458}
460 std::size_t k, x_N_dependency_flag xN_flag) {
461 double line1 = d_ndalphardni_dDelta(HEOS, i, xN_flag) * d2_nddeltadni_dxj_dDelta__consttau(HEOS, j, k, xN_flag)
462 + d2_ndalphardni_dDelta2(HEOS, i, xN_flag) * d_nddeltadni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
463 + d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, k, xN_flag) * d_nddeltadni_dDelta(HEOS, j, xN_flag)
464 + d3_ndalphardni_dxj_dDelta2__consttau_xi(HEOS, i, k, xN_flag) * nddeltadni__constT_V_nj(HEOS, j, xN_flag);
465
466 double line2 = d2_ndalphardni_dDelta_dTau(HEOS, i, xN_flag) * d_ndtaudni_dxj__constdelta_tau(HEOS, j, k, xN_flag)
467 + d3_ndalphardni_dxj_dDelta_dTau__constxi(HEOS, i, k, xN_flag) * ndtaudni__constT_V_nj(HEOS, j, xN_flag);
468
469 double line3 = d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(HEOS, i, j, k, xN_flag) - d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, k, xN_flag);
470 std::size_t mmax = HEOS.mole_fractions.size();
471 if (xN_flag == XN_DEPENDENT) {
472 mmax--;
473 }
474 for (unsigned int m = 0; m < mmax; m++) {
475 line3 -= HEOS.mole_fractions[m] * d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(HEOS, i, m, k, xN_flag);
476 }
477 return line1 + line2 + line3;
478}
480 // The first line
481 double term1 = (HEOS._delta.pt() * HEOS.d2alphar_dDelta2() + HEOS.dalphar_dDelta())
482 * (1 - 1 / HEOS._reducing.rhomolar * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
483
484 // The second line
485 double term2 =
486 HEOS._tau.pt() * HEOS.d2alphar_dDelta_dTau() * (1 / HEOS._reducing.T) * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
487
488 // The third line
489 double term3 = HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, i, xN_flag);
490 std::size_t kmax = HEOS.mole_fractions.size();
491 if (xN_flag == XN_DEPENDENT) {
492 kmax--;
493 }
494 for (unsigned int k = 0; k < kmax; k++) {
495 term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, k, xN_flag);
496 }
497 return term1 + term2 + term3;
498}
500 double term1 = (2 * HEOS.d2alphar_dDelta2() + HEOS.delta() * HEOS.d3alphar_dDelta3()) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
501 double term2 = HEOS.tau() * HEOS.d3alphar_dDelta2_dTau() * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
502 double term3 = HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, i, xN_flag);
503 std::size_t kmax = HEOS.mole_fractions.size();
504 if (xN_flag == XN_DEPENDENT) {
505 kmax--;
506 }
507 for (unsigned int k = 0; k < kmax; k++) {
508 term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, k, xN_flag);
509 }
510 return term1 + term2 + term3;
511}
513 double term1 = (3 * HEOS.d3alphar_dDelta3() + HEOS.delta() * HEOS.d4alphar_dDelta4()) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
514 double term2 = HEOS.tau() * HEOS.d4alphar_dDelta3_dTau() * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
515 double term3 = HEOS.residual_helmholtz->d4alphar_dxi_dDelta3(HEOS, i, xN_flag);
516 std::size_t kmax = HEOS.mole_fractions.size();
517 if (xN_flag == XN_DEPENDENT) {
518 kmax--;
519 }
520 for (unsigned int k = 0; k < kmax; k++) {
521 term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dDelta3(HEOS, k, xN_flag);
522 }
523 return term1 + term2 + term3;
524}
526 double term1 =
527 (HEOS.d2alphar_dDelta_dTau() + HEOS.delta() * HEOS.d3alphar_dDelta2_dTau()) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
528 double term2 = (HEOS.tau() * HEOS.d3alphar_dDelta_dTau2() + HEOS.d2alphar_dDelta_dTau()) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
529 double term3 = HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, i, xN_flag);
530 std::size_t kmax = HEOS.mole_fractions.size();
531 if (xN_flag == XN_DEPENDENT) {
532 kmax--;
533 }
534 for (unsigned int k = 0; k < kmax; k++) {
535 term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, k, xN_flag);
536 }
537 return term1 + term2 + term3;
538}
540 double term1 =
541 (2 * HEOS.d3alphar_dDelta2_dTau() + HEOS.delta() * HEOS.d4alphar_dDelta3_dTau()) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
542 double term2 =
543 (HEOS.tau() * HEOS.d4alphar_dDelta2_dTau2() + HEOS.d3alphar_dDelta2_dTau()) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
544 double term3 = HEOS.residual_helmholtz->d4alphar_dxi_dDelta2_dTau(HEOS, i, xN_flag);
545 std::size_t kmax = HEOS.mole_fractions.size();
546 if (xN_flag == XN_DEPENDENT) {
547 kmax--;
548 }
549 for (unsigned int k = 0; k < kmax; k++) {
550 term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dDelta2_dTau(HEOS, k, xN_flag);
551 }
552 return term1 + term2 + term3;
553}
555 double term1 =
556 (HEOS.d3alphar_dDelta_dTau2() + HEOS.delta() * HEOS.d4alphar_dDelta2_dTau2()) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
557 double term2 =
558 (HEOS.tau() * HEOS.d4alphar_dDelta_dTau3() + 2 * HEOS.d3alphar_dDelta_dTau2()) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
559 double term3 = HEOS.residual_helmholtz->d4alphar_dxi_dDelta_dTau2(HEOS, i, xN_flag);
560 std::size_t kmax = HEOS.mole_fractions.size();
561 if (xN_flag == XN_DEPENDENT) {
562 kmax--;
563 }
564 for (unsigned int k = 0; k < kmax; k++) {
565 term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dDelta_dTau2(HEOS, k, xN_flag);
566 }
567 return term1 + term2 + term3;
568}
570 double term1 = HEOS.delta() * HEOS.d3alphar_dDelta_dTau2() * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
571 double term2 = (2 * HEOS.d2alphar_dTau2() + HEOS.tau() * HEOS.d3alphar_dTau3()) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
572 double term3 = HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, i, xN_flag);
573 std::size_t kmax = HEOS.mole_fractions.size();
574 if (xN_flag == XN_DEPENDENT) {
575 kmax--;
576 }
577 for (unsigned int k = 0; k < kmax; k++) {
578 term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, k, xN_flag);
579 }
580 return term1 + term2 + term3;
581}
583 double term1 = HEOS.delta() * HEOS.d4alphar_dDelta_dTau3() * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
584 double term2 = (3 * HEOS.d3alphar_dTau3() + HEOS.tau() * HEOS.d4alphar_dTau4()) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
585 double term3 = HEOS.residual_helmholtz->d4alphar_dxi_dTau3(HEOS, i, xN_flag);
586 std::size_t kmax = HEOS.mole_fractions.size();
587 if (xN_flag == XN_DEPENDENT) {
588 kmax--;
589 }
590 for (unsigned int k = 0; k < kmax; k++) {
591 term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dTau3(HEOS, k, xN_flag);
592 }
593 return term1 + term2 + term3;
594}
595
597 x_N_dependency_flag xN_flag) {
598 double term1 =
599 (HEOS.dalphar_dDelta() + HEOS.delta() * HEOS.d2alphar_dDelta2()) * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
600 double term2 = (HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, j, xN_flag)
601 + HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, j, xN_flag))
602 * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
603 double term3 = HEOS.tau() * HEOS.d2alphar_dDelta_dTau() * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
604 double term4 =
605 HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
606 double term5 = HEOS.residual_helmholtz->d3alphar_dxi_dxj_dDelta(HEOS, i, j, xN_flag);
607 std::size_t kmax = HEOS.mole_fractions.size();
608 if (xN_flag == XN_DEPENDENT) {
609 kmax--;
610 }
611 for (unsigned int k = 0; k < kmax; k++) {
612 term5 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d3alphar_dxi_dxj_dDelta(HEOS, k, j, xN_flag)
613 + Kronecker_delta(k, j) * HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, k, xN_flag);
614 }
615 return term1 + term2 + term3 + term4 + term5;
616}
618 x_N_dependency_flag xN_flag) {
619 double term1 = HEOS.delta() * HEOS.d2alphar_dDelta_dTau() * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
620 double term2 =
621 HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
622 double term3 = (HEOS.tau() * HEOS.d2alphar_dTau2() + HEOS.dalphar_dTau()) * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
623 double term4 =
624 (HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, j, xN_flag) + HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, j, xN_flag))
625 * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
626 double term5 = HEOS.residual_helmholtz->d3alphar_dxi_dxj_dTau(HEOS, i, j, xN_flag);
627 std::size_t kmax = HEOS.mole_fractions.size();
628 if (xN_flag == XN_DEPENDENT) {
629 kmax--;
630 }
631 for (unsigned int k = 0; k < kmax; k++) {
632 term5 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d3alphar_dxi_dxj_dTau(HEOS, k, j, xN_flag)
633 + Kronecker_delta(k, j) * HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, k, xN_flag);
634 }
635 return term1 + term2 + term3 + term4 + term5;
636}
638 x_N_dependency_flag xN_flag) {
639 double term1 = HEOS.delta() * HEOS.d3alphar_dDelta_dTau2() * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
640 double term2 =
641 HEOS.delta() * HEOS.residual_helmholtz->d4alphar_dxi_dDelta_dTau2(HEOS, j, xN_flag) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
642 double term3 = (HEOS.tau() * HEOS.d3alphar_dTau3() + 2 * HEOS.d2alphar_dTau2()) * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
643 double term4 =
644 (HEOS.tau() * HEOS.residual_helmholtz->d4alphar_dxi_dTau3(HEOS, j, xN_flag) + 2 * HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, j, xN_flag))
645 * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
646 double term5 = HEOS.residual_helmholtz->d4alphar_dxi_dxj_dTau2(HEOS, i, j, xN_flag);
647 std::size_t kmax = HEOS.mole_fractions.size();
648 if (xN_flag == XN_DEPENDENT) {
649 kmax--;
650 }
651 for (unsigned int k = 0; k < kmax; k++) {
652 term5 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dTau2(HEOS, k, j, xN_flag)
653 + Kronecker_delta(k, j) * HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, k, xN_flag);
654 }
655 return term1 + term2 + term3 + term4 + term5;
656}
658 x_N_dependency_flag xN_flag) {
659 double term1 =
660 (2 * HEOS.d2alphar_dDelta2() + HEOS.delta() * HEOS.d3alphar_dDelta3()) * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
661 double term2 = (2 * HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, j, xN_flag)
662 + HEOS.delta() * HEOS.residual_helmholtz->d4alphar_dxi_dDelta3(HEOS, j, xN_flag))
663 * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
664 double term3 = HEOS.tau() * HEOS.d3alphar_dDelta2_dTau() * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
665 double term4 =
666 HEOS.tau() * HEOS.residual_helmholtz->d4alphar_dxi_dDelta2_dTau(HEOS, j, xN_flag) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
667 double term5 = HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta2(HEOS, i, j, xN_flag);
668 std::size_t kmax = HEOS.mole_fractions.size();
669 if (xN_flag == XN_DEPENDENT) {
670 kmax--;
671 }
672 for (unsigned int k = 0; k < kmax; k++) {
673 term5 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta2(HEOS, k, j, xN_flag)
674 + Kronecker_delta(k, j) * HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, k, xN_flag);
675 }
676 return term1 + term2 + term3 + term4 + term5;
677}
679 x_N_dependency_flag xN_flag) {
680 double term1 =
681 (HEOS.d2alphar_dDelta_dTau() + HEOS.delta() * HEOS.d3alphar_dDelta2_dTau()) * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
682 double term2 = (HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag)
683 + HEOS.delta() * HEOS.residual_helmholtz->d4alphar_dxi_dDelta2_dTau(HEOS, j, xN_flag))
684 * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
685 double term3 =
686 (HEOS.tau() * HEOS.d3alphar_dDelta_dTau2() + HEOS.d2alphar_dDelta_dTau()) * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
687 double term4 = (HEOS.tau() * HEOS.residual_helmholtz->d4alphar_dxi_dDelta_dTau2(HEOS, j, xN_flag)
688 + HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag))
689 * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
690 double term5 = HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta_dTau(HEOS, i, j, xN_flag);
691 std::size_t kmax = HEOS.mole_fractions.size();
692 if (xN_flag == XN_DEPENDENT) {
693 kmax--;
694 }
695 for (unsigned int k = 0; k < kmax; k++) {
696 term5 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta_dTau(HEOS, k, j, xN_flag)
697 + Kronecker_delta(k, j) * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, k, xN_flag);
698 }
699 return term1 + term2 + term3 + term4 + term5;
700}
701
703 // The first line
704 double term1 = HEOS._delta.pt() * HEOS.d2alphar_dDelta_dTau()
705 * (1 - 1 / HEOS._reducing.rhomolar * HEOS.Reducing->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
706
707 // The second line
708 double term2 = (HEOS._tau.pt() * HEOS.d2alphar_dTau2() + HEOS.dalphar_dTau()) * (1 / HEOS._reducing.T)
709 * HEOS.Reducing->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
710
711 // The third line
712 double term3 = HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, i, xN_flag);
713 std::size_t kmax = HEOS.mole_fractions.size();
714 if (xN_flag == XN_DEPENDENT) {
715 kmax--;
716 }
717 for (unsigned int k = 0; k < kmax; k++) {
718 term3 -= HEOS.mole_fractions[k] * HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, k, xN_flag);
719 }
720 return term1 + term2 + term3;
721}
722
724 std::size_t k, x_N_dependency_flag xN_flag) {
725 double term1 =
726 HEOS.delta()
727 * (HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, j, xN_flag) * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, k, xN_flag)
728 + HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, k, xN_flag) * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag));
729 double term2 =
730 HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dxj_dDelta(HEOS, j, k, xN_flag) * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
731 double term3 = HEOS.delta() * HEOS.dalphar_dDelta() * HEOS.Reducing->d2_PSI_rho_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag);
732
733 double term4 =
734 HEOS.tau()
735 * (HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, j, xN_flag) * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, k, xN_flag)
736 + HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, k, xN_flag) * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag));
737 double term5 =
738 HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
739 double term6 = HEOS.tau() * HEOS.dalphar_dTau() * HEOS.Reducing->d2_PSI_T_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag);
740
741 double term7 =
742 HEOS.residual_helmholtz->d3alphardxidxjdxk(HEOS, i, j, k, xN_flag) - 2 * HEOS.residual_helmholtz->d2alphardxidxj(HEOS, j, k, xN_flag);
743 std::size_t mmax = HEOS.mole_fractions.size();
744 if (xN_flag == XN_DEPENDENT) {
745 mmax--;
746 }
747 for (unsigned int m = 0; m < mmax; m++) {
748 term7 -= HEOS.mole_fractions[m] * HEOS.residual_helmholtz->d3alphardxidxjdxk(HEOS, j, k, m, xN_flag);
749 }
750
751 return term1 + term2 + term3 + term4 + term5 + term6 + term7;
752}
753
755 std::size_t k, x_N_dependency_flag xN_flag) {
756 double term1a = HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag)
757 * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, k, xN_flag);
758 double term1b = HEOS.delta() * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta_dTau(HEOS, j, k, xN_flag)
759 * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
760 double term1c = HEOS.delta() * HEOS.d2alphar_dDelta_dTau() * HEOS.Reducing->d2_PSI_rho_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag);
761 double term1d = HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, k, xN_flag)
762 * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
763 double term1 = term1a + term1b + term1c + term1d;
764
765 double term2a =
766 (HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, j, xN_flag) + HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, j, xN_flag))
767 * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, k, xN_flag);
768 double term2b = (HEOS.tau() * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dTau2(HEOS, j, k, xN_flag)
769 + HEOS.residual_helmholtz->d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag))
770 * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
771 double term2c =
772 (HEOS.tau() * HEOS.d2alphar_dTau2() + HEOS.dalphar_dTau()) * HEOS.Reducing->d2_PSI_T_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag);
773 double term2d =
774 (HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dTau2(HEOS, k, xN_flag) + HEOS.residual_helmholtz->d2alphar_dxi_dTau(HEOS, k, xN_flag))
775 * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
776 double term2 = term2a + term2b + term2c + term2d;
777
778 double term3 = HEOS.residual_helmholtz->d4alphar_dxi_dxj_dxk_dTau(HEOS, i, j, k, xN_flag)
779 - 2 * HEOS.residual_helmholtz->d3alphar_dxi_dxj_dTau(HEOS, j, k, xN_flag);
780 std::size_t mmax = HEOS.mole_fractions.size();
781 if (xN_flag == XN_DEPENDENT) {
782 mmax--;
783 }
784 for (unsigned int m = 0; m < mmax; m++) {
785 term3 -= HEOS.mole_fractions[m] * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dxk_dTau(HEOS, j, k, m, xN_flag);
786 }
787
788 return term1 + term2 + term3;
789}
790
792 std::size_t k, x_N_dependency_flag xN_flag) {
793 double term1a = (HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, j, xN_flag)
794 + HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, j, xN_flag))
795 * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, k, xN_flag);
796 double term1b = (HEOS.delta() * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta2(HEOS, j, k, xN_flag)
797 + HEOS.residual_helmholtz->d3alphar_dxi_dxj_dDelta(HEOS, j, k, xN_flag))
798 * HEOS.Reducing->PSI_rho(HEOS.mole_fractions, i, xN_flag);
799 double term1c =
800 (HEOS.delta() * HEOS.d2alphar_dDelta2() + HEOS.dalphar_dDelta()) * HEOS.Reducing->d2_PSI_rho_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag);
801 double term1d = (HEOS.delta() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta2(HEOS, k, xN_flag)
802 + HEOS.residual_helmholtz->d2alphar_dxi_dDelta(HEOS, k, xN_flag))
803 * HEOS.Reducing->d_PSI_rho_dxj(HEOS.mole_fractions, i, j, xN_flag);
804 double term1 = term1a + term1b + term1c + term1d;
805
806 double term2a = HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, j, xN_flag)
807 * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, k, xN_flag);
808 double term2b =
809 HEOS.tau() * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dDelta_dTau(HEOS, j, k, xN_flag) * HEOS.Reducing->PSI_T(HEOS.mole_fractions, i, xN_flag);
810 double term2c = HEOS.tau() * HEOS.d2alphar_dDelta_dTau() * HEOS.Reducing->d2_PSI_T_dxj_dxk(HEOS.mole_fractions, i, j, k, xN_flag);
811 double term2d = HEOS.tau() * HEOS.residual_helmholtz->d3alphar_dxi_dDelta_dTau(HEOS, k, xN_flag)
812 * HEOS.Reducing->d_PSI_T_dxj(HEOS.mole_fractions, i, j, xN_flag);
813 double term2 = term2a + term2b + term2c + term2d;
814
815 double term3 = HEOS.residual_helmholtz->d4alphar_dxi_dxj_dxk_dDelta(HEOS, i, j, k, xN_flag)
816 - 2 * HEOS.residual_helmholtz->d3alphar_dxi_dxj_dDelta(HEOS, j, k, xN_flag);
817 std::size_t mmax = HEOS.mole_fractions.size();
818 if (xN_flag == XN_DEPENDENT) {
819 mmax--;
820 }
821 for (unsigned int m = 0; m < mmax; m++) {
822 term3 -= HEOS.mole_fractions[m] * HEOS.residual_helmholtz->d4alphar_dxi_dxj_dxk_dDelta(HEOS, j, k, m, xN_flag);
823 }
824 return term1 + term2 + term3;
825}
826
828 // Reducing values are constant for all components under consideration
829 double Tr = HEOS.T_reducing();
830 double rhor = HEOS.rhomolar_reducing();
831
832 // Values for the i-th component
833 double Tci = HEOS.get_fluid_constant(i, iT_critical);
834 double rhoci = HEOS.get_fluid_constant(i, irhomolar_critical);
835 double tau_oi = HEOS.tau() * Tci / Tr;
836 double delta_oi = HEOS.delta() * rhor / rhoci;
837 double Rratioi = 1; //HEOS.gas_constant()/HEOS.components[i].EOS().R_u;
838
839 double logxi = (std::abs(HEOS.mole_fractions[i]) > DBL_EPSILON) ? (log(HEOS.mole_fractions[i])) : 0;
840 double term = Rratioi * HEOS.components[i].EOS().alpha0.base(tau_oi, delta_oi) + logxi + 1;
841
842 std::size_t kmax = HEOS.mole_fractions.size();
843 if (xN_flag == XN_DEPENDENT) {
844 kmax--;
845 }
846 for (std::size_t k = 0; k < kmax; ++k) {
847 double xk = HEOS.mole_fractions[k];
848 double Tck = HEOS.get_fluid_constant(k, iT_critical);
849 double rhock = HEOS.get_fluid_constant(k, irhomolar_critical);
850 double tau_ok = HEOS.tau() * Tck / Tr;
851 double delta_ok = HEOS.delta() * rhor / rhock;
852 double dtauok_dxi = -tau_ok / Tr * HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, i, xN_flag); // (Gernert, supp, B.19)
853 double ddeltaok_dxi = delta_ok / rhor * HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, i, xN_flag); // (Gernert, supp. B.20)
854
855 double Rratiok = 1; //HEOS.gas_constant()/HEOS.components[k].EOS().R_u;
856 HelmholtzDerivatives alpha0kterms = HEOS.components[k].EOS().alpha0.all(tau_ok, delta_ok);
857 double dalpha0_ok_dxi = alpha0kterms.dalphar_dtau * dtauok_dxi + alpha0kterms.dalphar_ddelta * ddeltaok_dxi;
858 term += xk * (Rratiok * dalpha0_ok_dxi);
859 }
860 return term;
861}
862
864 // Reducing values are constant for all components under consideration
865 double Tr = HEOS.T_reducing();
866 double rhor = HEOS.rhomolar_reducing();
867
868 // Values for the i-th component
869 double Tci = HEOS.get_fluid_constant(i, iT_critical);
870 double rhoci = HEOS.get_fluid_constant(i, irhomolar_critical);
871 double tau_oi = HEOS.tau() * Tci / Tr;
872 double delta_oi = HEOS.delta() * rhor / rhoci;
873 double Rratioi = 1; //HEOS.gas_constant()/HEOS.components[i].EOS().R_u;
874
875 double term = rhor / rhoci * Rratioi * HEOS.components[i].EOS().alpha0.dDelta(tau_oi, delta_oi);
876
877 std::size_t kmax = HEOS.mole_fractions.size();
878 if (xN_flag == XN_DEPENDENT) {
879 kmax--;
880 }
881 for (std::size_t k = 0; k < kmax; ++k) {
882 double xk = HEOS.mole_fractions[k];
883 double Tck = HEOS.get_fluid_constant(k, iT_critical);
884 double rhock = HEOS.get_fluid_constant(k, irhomolar_critical);
885 double tau_ok = HEOS.tau() * Tck / Tr;
886 double delta_ok = HEOS.delta() * rhor / rhock;
887 double dtauok_dxi = -tau_ok / Tr * HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, i, xN_flag); // (Gernert, supp, B.19)
888 double drhor_dxi = HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, i, xN_flag);
889 double ddeltaok_dxi = delta_ok / rhor * drhor_dxi; // (Gernert, supp. B.20)
890
891 //double Rratiok = 1;//HEOS.gas_constant()/HEOS.components[k].EOS().R_u;
892 HelmholtzDerivatives alpha0kterms = HEOS.components[k].EOS().alpha0.all(tau_ok, delta_ok);
893 double dalpha0ok_ddeltaok = alpha0kterms.dalphar_ddelta;
894
895 double d_dalpha0ok_ddeltaok_dxi = alpha0kterms.d2alphar_ddelta_dtau * dtauok_dxi + alpha0kterms.d2alphar_ddelta2 * ddeltaok_dxi;
896 term += xk / rhock * (rhor * d_dalpha0ok_ddeltaok_dxi + drhor_dxi * dalpha0ok_ddeltaok);
897 }
898 return term;
899}
900
902 // Reducing values are constant for all components under consideration
903 double Tr = HEOS.T_reducing();
904 double rhor = HEOS.rhomolar_reducing();
905
906 // Values for the i-th component
907 double Tci = HEOS.get_fluid_constant(i, iT_critical);
908 double rhoci = HEOS.get_fluid_constant(i, irhomolar_critical);
909 double tau_oi = HEOS.tau() * Tci / Tr;
910 double delta_oi = HEOS.delta() * rhor / rhoci;
911 double Rratioi = 1; //HEOS.gas_constant()/HEOS.components[i].EOS().R_u;
912
913 double term = Tci / Tr * Rratioi * HEOS.components[i].EOS().alpha0.dTau(tau_oi, delta_oi);
914
915 std::size_t kmax = HEOS.mole_fractions.size();
916 if (xN_flag == XN_DEPENDENT) {
917 kmax--;
918 }
919 for (std::size_t k = 0; k < kmax; ++k) {
920 double xk = HEOS.mole_fractions[k];
921 double Tck = HEOS.get_fluid_constant(k, iT_critical);
922 double rhock = HEOS.get_fluid_constant(k, irhomolar_critical);
923 double tau_ok = HEOS.tau() * Tck / Tr;
924 double delta_ok = HEOS.delta() * rhor / rhock;
925 double dTr_dxi = HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, i, xN_flag);
926 double dtauok_dxi = -tau_ok / Tr * dTr_dxi; // (Gernert, supp, B.19)
927 double drhor_dxi = HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, i, xN_flag);
928 double ddeltaok_dxi = delta_ok / rhor * drhor_dxi; // (Gernert, supp. B.20)
929
930 //double Rratiok = 1;//HEOS.gas_constant()/HEOS.components[k].EOS().R_u;
931 HelmholtzDerivatives alpha0kterms = HEOS.components[k].EOS().alpha0.all(tau_ok, delta_ok);
932 double dalpha0ok_dtauok = alpha0kterms.dalphar_dtau;
933 double d_dalpha0ok_dTauok_dxi = alpha0kterms.d2alphar_dtau2 * dtauok_dxi + alpha0kterms.d2alphar_ddelta_dtau * ddeltaok_dxi;
934 term += xk * Tck * (1 / Tr * d_dalpha0ok_dTauok_dxi + -1 / POW2(Tr) * dTr_dxi * dalpha0ok_dtauok);
935 }
936 return term;
937}
938
940 // Reducing values are constant for all components under consideration
941 double Tr = HEOS.T_reducing();
942 double rhor = HEOS.rhomolar_reducing();
943
944 // Values for the i-th component
945 double Tci = HEOS.get_fluid_constant(i, iT_critical);
946 double rhoci = HEOS.get_fluid_constant(i, irhomolar_critical);
947 double tau_oi = HEOS.tau() * Tci / Tr;
948 double delta_oi = HEOS.delta() * rhor / rhoci;
949 double dTr_dxi = HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, i, xN_flag);
950 double drhor_dxi = HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, i, xN_flag);
951
952 // Values for the j-th component
953 double Tcj = HEOS.get_fluid_constant(j, iT_critical);
954 double rhocj = HEOS.get_fluid_constant(j, irhomolar_critical);
955 double tau_oj = HEOS.tau() * Tcj / Tr;
956 double delta_oj = HEOS.delta() * rhor / rhocj;
957 double dTr_dxj = HEOS.Reducing->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag);
958 double drhor_dxj = HEOS.Reducing->drhormolardxi__constxj(HEOS.mole_fractions, j, xN_flag);
959
960 // Cross-terms with both i & j
961 double dtauoi_dxj = -tau_oi / Tr * dTr_dxj; // (Gernert, supp, B.19)
962 double ddeltaoi_dxj = delta_oi / rhor * drhor_dxj; // (Gernert, supp. B.20)
963 double dtauoj_dxi = -tau_oj / Tr * dTr_dxi; // (Gernert, supp, B.19)
964 double ddeltaoj_dxi = delta_oj / rhor * drhor_dxi; // (Gernert, supp. B.20)
965 double d2Tr_dxidxj = HEOS.Reducing->d2Trdxidxj(HEOS.mole_fractions, i, j, xN_flag);
966 double d2rhor_dxidxj = HEOS.Reducing->d2rhormolardxidxj(HEOS.mole_fractions, i, j, xN_flag);
967
968 //double Rratioi = 1;//HEOS.gas_constant()/HEOS.components[i].EOS().R_u;
969 HelmholtzDerivatives alpha0iterms = HEOS.components[i].EOS().alpha0.all(tau_oi, delta_oi),
970 alpha0jterms = HEOS.components[j].EOS().alpha0.all(tau_oj, delta_oj);
971
972 double d_dalpha0oi_dxj = alpha0iterms.dalphar_dtau * dtauoi_dxj + alpha0iterms.dalphar_ddelta * ddeltaoi_dxj;
973 double d_dalpha0oj_dxi = alpha0jterms.dalphar_dtau * dtauoj_dxi + alpha0jterms.dalphar_ddelta * ddeltaoj_dxi;
974
975 double xi = HEOS.mole_fractions[i], xj = HEOS.mole_fractions[j];
976 double Kronecker_delta_over_xi = (xj > DBL_EPSILON && xi > DBL_EPSILON) ? Kronecker_delta(i, j) / xi : 0;
977 double term = d_dalpha0oi_dxj + d_dalpha0oj_dxi + Kronecker_delta_over_xi;
978
979 std::size_t kmax = HEOS.mole_fractions.size();
980 if (xN_flag == XN_DEPENDENT) {
981 kmax--;
982 }
983 for (std::size_t k = 0; k < kmax; ++k) {
984 // Values for the k-th component
985 double xk = HEOS.mole_fractions[k];
986 double Tck = HEOS.get_fluid_constant(k, iT_critical);
987 double rhock = HEOS.get_fluid_constant(k, irhomolar_critical);
988 double tau_ok = HEOS.tau() * Tck / Tr;
989 double delta_ok = HEOS.delta() * rhor / rhock;
990
991 double dtauok_dxj = -tau_ok / Tr * dTr_dxj; // (Gernert, supp, B.19)
992 double ddeltaok_dxj = delta_ok / rhor * drhor_dxj; // (Gernert, supp. B.20)
993 double dtauok_dxi = -tau_ok / Tr * dTr_dxi; // (Gernert, supp, B.19)
994 double ddeltaok_dxi = delta_ok / rhor * drhor_dxi; // (Gernert, supp. B.20)
995
996 HelmholtzDerivatives alpha0kterms = HEOS.components[k].EOS().alpha0.all(tau_ok, delta_ok);
997 double dalpha0ok_dtauok = alpha0kterms.dalphar_dtau;
998 double d2tauok_dxidxj = -Tck * HEOS.tau() * (POW2(Tr) * d2Tr_dxidxj - dTr_dxi * (2 * Tr * dTr_dxj)) / POW4(Tr);
999 double d_dalpha0ok_dtauok_dxj = alpha0kterms.d2alphar_dtau2 * dtauok_dxj + alpha0kterms.d2alphar_ddelta_dtau * ddeltaok_dxj;
1000
1001 double dalpha0ok_ddeltaok = alpha0kterms.dalphar_ddelta;
1002 double d2deltaok_dxidxj = HEOS.delta() / rhock * d2rhor_dxidxj;
1003 double d_dalpha0ok_ddeltaok_dxj = alpha0kterms.d2alphar_ddelta_dtau * dtauok_dxj + alpha0kterms.d2alphar_ddelta2 * ddeltaok_dxj;
1004
1005 term += xk
1006 * (dalpha0ok_dtauok * d2tauok_dxidxj + d_dalpha0ok_dtauok_dxj * dtauok_dxi + dalpha0ok_ddeltaok * d2deltaok_dxidxj
1007 + d_dalpha0ok_ddeltaok_dxj * ddeltaok_dxi);
1008 }
1009 return term;
1010}
1011
1014 return HEOS.rhomolar_reducing() * HEOS.gas_constant() * HEOS.T() * (HEOS.delta() * dalpha_dDelta(HEOS) + alpha(HEOS));
1015}
1016
1019 return HEOS.rhomolar_reducing() * HEOS.delta() * HEOS.gas_constant() * HEOS.T() / HEOS.tau() * (HEOS.tau() * dalpha_dTau(HEOS) - alpha(HEOS));
1020}
1021
1024 return HEOS.rhomolar_reducing() * HEOS.delta() * HEOS.gas_constant() * HEOS.T() / HEOS.tau() * (HEOS.tau() * dalphar_dTau(HEOS) - alphar(HEOS));
1025}
1026
1028 return HEOS.delta() * HEOS.gas_constant() / HEOS.tau()
1029 * (alpha(HEOS, xN_flag) * d_rhorTr_dxi(HEOS, i, xN_flag) + HEOS.rhomolar_reducing() * HEOS.T_reducing() * dalpha_dxi(HEOS, i, xN_flag));
1030}
1031
1033 return HEOS.delta() * HEOS.gas_constant() / HEOS.tau()
1034 * (alphar(HEOS, xN_flag) * d_rhorTr_dxi(HEOS, i, xN_flag) + HEOS.rhomolar_reducing() * HEOS.T_reducing() * dalphar_dxi(HEOS, i, xN_flag));
1035}
1036
1038 GERG2008ReducingFunction* GERG = static_cast<GERG2008ReducingFunction*>(HEOS.Reducing.get());
1039 return HEOS.rhomolar_reducing() * GERG->dTrdxi__constxj(HEOS.mole_fractions, i, xN_flag)
1040 + HEOS.T_reducing() * GERG->drhormolardxi__constxj(HEOS.mole_fractions, i, xN_flag);
1041}
1042
1044 GERG2008ReducingFunction* GERG = static_cast<GERG2008ReducingFunction*>(HEOS.Reducing.get());
1045 return (HEOS.rhomolar_reducing() * GERG->d2Trdxidxj(HEOS.mole_fractions, i, j, xN_flag)
1046 + GERG->drhormolardxi__constxj(HEOS.mole_fractions, j, xN_flag) * GERG->dTrdxi__constxj(HEOS.mole_fractions, i, xN_flag)
1047 + HEOS.T_reducing() * GERG->d2rhormolardxidxj(HEOS.mole_fractions, i, j, xN_flag)
1048 + GERG->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag) * GERG->drhormolardxi__constxj(HEOS.mole_fractions, i, xN_flag));
1049}
1050
1052 return HEOS.rhomolar_reducing() * HEOS.gas_constant() * HEOS.T() * (HEOS.delta() * d2alpha_dDelta2(HEOS) + 2 * dalpha_dDelta(HEOS));
1053}
1054
1056 return HEOS.rhomolar_reducing() * HEOS.gas_constant() * HEOS.T() / HEOS.tau()
1057 * (HEOS.tau() * dalpha_dTau(HEOS) - alpha(HEOS) - HEOS.delta() * dalpha_dDelta(HEOS)
1058 + HEOS.tau() * HEOS.delta() * d2alpha_dDelta_dTau(HEOS));
1059}
1061 return HEOS.rhomolar_reducing() * HEOS.gas_constant() * HEOS.T() / HEOS.tau()
1062 * (HEOS.tau() * dalphar_dTau(HEOS) - alphar(HEOS) - HEOS.delta() * dalphar_dDelta(HEOS)
1063 + HEOS.tau() * HEOS.delta() * d2alphar_dDelta_dTau(HEOS));
1064}
1066 double tau = HEOS.tau();
1067 return HEOS.rhomolar_reducing() * HEOS.delta() * HEOS.gas_constant() * HEOS.T() / POW2(tau)
1068 * (POW2(tau) * d2alpha_dTau2(HEOS) - 2 * tau * dalpha_dTau(HEOS) + 2 * alpha(HEOS));
1069}
1071 double tau = HEOS.tau();
1072 return HEOS.rhomolar_reducing() * HEOS.delta() * HEOS.gas_constant() * HEOS.T() / POW2(tau)
1073 * (POW2(tau) * d2alphar_dTau2(HEOS) - 2 * tau * dalphar_dTau(HEOS) + 2 * alphar(HEOS));
1074}
1076 return HEOS.gas_constant() / HEOS.tau()
1077 * (d_rhorTr_dxi(HEOS, i, xN_flag) * (HEOS.delta() * dalpha_dDelta(HEOS) + alpha(HEOS))
1078 + HEOS.rhomolar_reducing() * HEOS.T_reducing() * (HEOS.delta() * d2alpha_dxi_dDelta(HEOS, i, xN_flag) + dalpha_dxi(HEOS, i, xN_flag)));
1079}
1081 return HEOS.delta() * HEOS.gas_constant() / POW2(HEOS.tau())
1082 * (d_rhorTr_dxi(HEOS, i, xN_flag) * (HEOS.tau() * dalpha_dTau(HEOS) - alpha(HEOS))
1083 + HEOS.rhomolar_reducing() * HEOS.T_reducing() * (HEOS.tau() * d2alpha_dxi_dTau(HEOS, i, xN_flag) - dalpha_dxi(HEOS, i, xN_flag)));
1084}
1086 return HEOS.delta() * HEOS.gas_constant() / POW2(HEOS.tau())
1087 * (d_rhorTr_dxi(HEOS, i, xN_flag) * (HEOS.tau() * dalphar_dTau(HEOS) - alphar(HEOS))
1088 + HEOS.rhomolar_reducing() * HEOS.T_reducing() * (HEOS.tau() * d2alphar_dxi_dTau(HEOS, i, xN_flag) - dalphar_dxi(HEOS, i, xN_flag)));
1089}
1091 return HEOS.delta() * HEOS.gas_constant() / HEOS.tau()
1092 * (alpha(HEOS) * d2_rhorTr_dxidxj(HEOS, i, j, xN_flag) + dalpha_dxi(HEOS, i, xN_flag) * d_rhorTr_dxi(HEOS, j, xN_flag)
1093 + dalpha_dxi(HEOS, j, xN_flag) * d_rhorTr_dxi(HEOS, i, xN_flag)
1094 + HEOS.rhomolar_reducing() * HEOS.T_reducing() * d2alphadxidxj(HEOS, i, j, xN_flag));
1095}
1097 return HEOS.delta() * HEOS.gas_constant() / HEOS.tau()
1098 * (alphar(HEOS) * d2_rhorTr_dxidxj(HEOS, i, j, xN_flag) + dalphar_dxi(HEOS, i, xN_flag) * d_rhorTr_dxi(HEOS, j, xN_flag)
1099 + dalphar_dxi(HEOS, j, xN_flag) * d_rhorTr_dxi(HEOS, i, xN_flag)
1100 + HEOS.rhomolar_reducing() * HEOS.T_reducing() * d2alphardxidxj(HEOS, i, j, xN_flag));
1101}
1102
1103} /* namespace CoolProp */
1104
1105#ifdef ENABLE_CATCH
1106# include <catch2/catch_all.hpp>
1107
1110
1111using namespace CoolProp;
1112
1113double mix_deriv_err_func(double numeric, double analytic) {
1114 if (std::abs(analytic) < DBL_EPSILON) {
1115 return std::abs(numeric - analytic);
1116 } else {
1117 return std::abs(numeric / analytic - 1);
1118 }
1119}
1120
1122
1123enum derivative
1124{
1125 NO_DERIV = 0,
1126 TAU,
1127 DELTA,
1128 XI,
1129 XJ,
1130 XK,
1131 T_CONSTP,
1132 P_CONSTT,
1133 T_CONSTRHO,
1134 RHO_CONSTT,
1135 CONST_TAU_DELTA,
1136 CONST_TRHO
1137};
1138
1139typedef double (*zero_mole_fraction_pointer)(CoolProp::HelmholtzEOSMixtureBackend& HEOS, CoolProp::x_N_dependency_flag xN_flag);
1140typedef double (*one_mole_fraction_pointer)(CoolProp::HelmholtzEOSMixtureBackend& HEOS, std::size_t i, CoolProp::x_N_dependency_flag xN_flag);
1141typedef double (*two_mole_fraction_pointer)(CoolProp::HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j,
1143typedef double (*three_mole_fraction_pointer)(CoolProp::HelmholtzEOSMixtureBackend& HEOS, std::size_t i, std::size_t j, std::size_t k,
1145
1146template <class backend>
1147class DerivativeFixture
1148{
1149 public:
1150 shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS, HEOS_plus_tau, HEOS_minus_tau, HEOS_plus_delta, HEOS_minus_delta, HEOS_plus_T__constp,
1151 HEOS_minus_T__constp, HEOS_plus_p__constT, HEOS_minus_p__constT, HEOS_plus_T__constrho, HEOS_minus_T__constrho, HEOS_plus_rho__constT,
1152 HEOS_minus_rho__constT;
1153 std::vector<shared_ptr<CoolProp::HelmholtzEOSMixtureBackend>> HEOS_plus_z, HEOS_minus_z, HEOS_plus_z__constTrho, HEOS_minus_z__constTrho,
1154 HEOS_plus_n, HEOS_minus_n, HEOS_plus_2z, HEOS_minus_2z, HEOS_plus_2z__constTrho, HEOS_minus_2z__constTrho;
1156 double dtau, ddelta, dz, dn, tol, dT, drho, dp;
1157 DerivativeFixture() : xN(XN_INDEPENDENT) {
1158 dtau = 1e-6;
1159 ddelta = 1e-6;
1160 dz = 1e-6;
1161 dn = 1e-6;
1162 dT = 1e-3;
1163 drho = 1e-3;
1164 dp = 1;
1165 tol = 5e-6;
1166 std::vector<std::string> names;
1167 names.push_back("n-Pentane");
1168 names.push_back("Ethane");
1169 names.push_back("n-Propane");
1170 names.push_back("n-Butane");
1171 std::vector<CoolPropDbl> mole_fractions;
1172 mole_fractions.push_back(0.1);
1173 mole_fractions.push_back(0.12);
1174 mole_fractions.push_back(0.18);
1175 mole_fractions.push_back(0.6);
1176 HEOS.reset(new backend(names));
1177 HEOS->set_mole_fractions(mole_fractions);
1178 HEOS->specify_phase(CoolProp::iphase_gas);
1179 HEOS->update_DmolarT_direct(300, 300);
1180
1181 init_state(HEOS_plus_tau);
1182 init_state(HEOS_minus_tau);
1183 init_state(HEOS_plus_delta);
1184 init_state(HEOS_minus_delta);
1185 init_state(HEOS_plus_T__constp);
1186 init_state(HEOS_minus_T__constp);
1187 init_state(HEOS_plus_p__constT);
1188 init_state(HEOS_minus_p__constT);
1189 init_state(HEOS_plus_T__constrho);
1190 init_state(HEOS_minus_T__constrho);
1191 init_state(HEOS_plus_rho__constT);
1192 init_state(HEOS_minus_rho__constT);
1193 // Constant composition, varying state
1194 HEOS_plus_tau->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS->rhomolar_reducing(), HEOS->T_reducing() / (HEOS->tau() + dtau));
1195 HEOS_minus_tau->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS->rhomolar_reducing(), HEOS->T_reducing() / (HEOS->tau() - dtau));
1196 HEOS_plus_delta->update(CoolProp::DmolarT_INPUTS, (HEOS->delta() + ddelta) * HEOS->rhomolar_reducing(), HEOS->T_reducing() / HEOS->tau());
1197 HEOS_minus_delta->update(CoolProp::DmolarT_INPUTS, (HEOS->delta() - ddelta) * HEOS->rhomolar_reducing(), HEOS->T_reducing() / HEOS->tau());
1198 HEOS_plus_T__constp->update(CoolProp::PT_INPUTS, HEOS->p(), HEOS->T() + dT);
1199 HEOS_minus_T__constp->update(CoolProp::PT_INPUTS, HEOS->p(), HEOS->T() - dT);
1200 HEOS_plus_p__constT->update(CoolProp::PT_INPUTS, HEOS->p() + dp, HEOS->T());
1201 HEOS_minus_p__constT->update(CoolProp::PT_INPUTS, HEOS->p() - dp, HEOS->T());
1202 HEOS_plus_T__constrho->update(CoolProp::DmolarT_INPUTS, HEOS->rhomolar(), HEOS->T() + dT);
1203 HEOS_minus_T__constrho->update(CoolProp::DmolarT_INPUTS, HEOS->rhomolar(), HEOS->T() - dT);
1204 HEOS_plus_rho__constT->update(CoolProp::DmolarT_INPUTS, HEOS->rhomolar() + drho, HEOS->T());
1205 HEOS_minus_rho__constT->update(CoolProp::DmolarT_INPUTS, HEOS->rhomolar() - drho, HEOS->T());
1206
1207 // Varying mole fractions
1208 HEOS_plus_z.resize(4);
1209 HEOS_minus_z.resize(4);
1210 HEOS_plus_z__constTrho.resize(4);
1211 HEOS_minus_z__constTrho.resize(4);
1212 HEOS_plus_2z.resize(4);
1213 HEOS_minus_2z.resize(4);
1214 HEOS_plus_2z__constTrho.resize(4);
1215 HEOS_minus_2z__constTrho.resize(4);
1216 for (int i = 0; i < HEOS_plus_z.size(); ++i) {
1217 init_state(HEOS_plus_z[i]);
1218 init_state(HEOS_plus_2z[i]);
1219 init_state(HEOS_plus_z__constTrho[i]);
1220 std::vector<double> zz = HEOS->get_mole_fractions(), zz2 = HEOS->get_mole_fractions();
1221 zz[i] += dz;
1222 zz2[i] += 2 * dz;
1223 if (xN == CoolProp::XN_DEPENDENT) {
1224 zz[zz.size() - 1] -= dz;
1225 zz2[zz2.size() - 1] -= 2 * dz;
1226 }
1227 HEOS_plus_z[i]->set_mole_fractions(zz);
1228 HEOS_plus_z[i]->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS_plus_z[i]->rhomolar_reducing(),
1229 HEOS_plus_z[i]->T_reducing() / HEOS->tau());
1230 HEOS_plus_2z[i]->set_mole_fractions(zz2);
1231 HEOS_plus_2z[i]->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS_plus_2z[i]->rhomolar_reducing(),
1232 HEOS_plus_2z[i]->T_reducing() / HEOS->tau());
1233 HEOS_plus_z__constTrho[i]->set_mole_fractions(zz);
1234 HEOS_plus_z__constTrho[i]->update(CoolProp::DmolarT_INPUTS, HEOS->rhomolar(), HEOS->T());
1235 }
1236 for (int i = 0; i < HEOS_minus_z.size(); ++i) {
1237 init_state(HEOS_minus_z[i]);
1238 init_state(HEOS_minus_2z[i]);
1239 init_state(HEOS_minus_z__constTrho[i]);
1240 std::vector<double> zz = HEOS->get_mole_fractions(), zz2 = HEOS->get_mole_fractions();
1241 zz[i] -= dz;
1242 zz2[i] -= 2 * dz;
1243 if (xN == CoolProp::XN_DEPENDENT) {
1244 zz[zz.size() - 1] += dz;
1245 zz2[zz2.size() - 1] += 2 * dz;
1246 }
1247 HEOS_minus_z[i]->set_mole_fractions(zz);
1248 HEOS_minus_z[i]->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS_minus_z[i]->rhomolar_reducing(),
1249 HEOS_minus_z[i]->T_reducing() / HEOS->tau());
1250 HEOS_minus_2z[i]->set_mole_fractions(zz2);
1251 HEOS_minus_2z[i]->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS_minus_2z[i]->rhomolar_reducing(),
1252 HEOS_minus_2z[i]->T_reducing() / HEOS->tau());
1253 HEOS_minus_z__constTrho[i]->set_mole_fractions(zz);
1254 HEOS_minus_z__constTrho[i]->update(CoolProp::DmolarT_INPUTS, HEOS->rhomolar(), HEOS->T());
1255 }
1256
1257 // Varying mole numbers
1258 HEOS_plus_n.resize(4);
1259 HEOS_minus_n.resize(4);
1260 for (int i = 0; i < HEOS_plus_n.size(); ++i) {
1261 init_state(HEOS_plus_n[i]);
1262 std::vector<double> zz = HEOS->get_mole_fractions();
1263 zz[i] += dn;
1264 for (int j = 0; j < HEOS_minus_n.size(); ++j) {
1265 zz[i] /= 1 + dn;
1266 }
1267 HEOS_plus_n[i]->set_mole_fractions(zz);
1268 HEOS_plus_n[i]->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS_plus_n[i]->rhomolar_reducing(),
1269 HEOS_plus_n[i]->T_reducing() / HEOS->tau());
1270 }
1271 for (int i = 0; i < HEOS_plus_z.size(); ++i) {
1272 init_state(HEOS_minus_n[i]);
1273 std::vector<double> zz = HEOS->get_mole_fractions();
1274 zz[i] -= dn;
1275 for (int j = 0; j < HEOS_minus_n.size(); ++j) {
1276 zz[i] /= 1 - dn;
1277 }
1278 HEOS_minus_n[i]->set_mole_fractions(zz);
1279 HEOS_minus_n[i]->update(CoolProp::DmolarT_INPUTS, HEOS->delta() * HEOS_minus_n[i]->rhomolar_reducing(),
1280 HEOS_minus_n[i]->T_reducing() / HEOS->tau());
1281 }
1282 };
1283 void init_state(shared_ptr<CoolProp::HelmholtzEOSMixtureBackend>& other) {
1284 other.reset(HEOS->get_copy());
1285 other->set_mole_fractions(HEOS->get_mole_fractions());
1286 other->specify_phase(CoolProp::iphase_gas); // Something homogeneous
1287 }
1288 void zero(const std::string& name, zero_mole_fraction_pointer f, zero_mole_fraction_pointer g = NULL, derivative wrt = NO_DERIV) {
1289 double analytic = f(*HEOS, xN);
1290 double numeric = 0;
1291 if (wrt == TAU) {
1292 numeric = (g(*HEOS_plus_tau, xN) - g(*HEOS_minus_tau, xN)) / (2 * dtau);
1293 } else if (wrt == DELTA) {
1294 numeric = (g(*HEOS_plus_delta, xN) - g(*HEOS_minus_delta, xN)) / (2 * ddelta);
1295 }
1296 CAPTURE(name);
1297 CAPTURE(analytic);
1298 CAPTURE(numeric);
1299 CAPTURE(xN);
1300 double error = mix_deriv_err_func(numeric, analytic);
1301 CAPTURE(error);
1302 CHECK(error < tol);
1303 }
1304 void one(const std::string& name, one_mole_fraction_pointer f, one_mole_fraction_pointer g = NULL, derivative wrt = NO_DERIV) {
1305 for (int i = 0; i < 4; ++i) {
1306 double analytic = f(*HEOS, i, xN);
1307 double numeric = 0;
1308 if (wrt == TAU) {
1309 numeric = (g(*HEOS_plus_tau, i, xN) - g(*HEOS_minus_tau, i, xN)) / (2 * dtau);
1310 } else if (wrt == DELTA) {
1311 numeric = (g(*HEOS_plus_delta, i, xN) - g(*HEOS_minus_delta, i, xN)) / (2 * ddelta);
1312 } else if (wrt == T_CONSTP) {
1313 numeric = (g(*HEOS_plus_T__constp, i, xN) - g(*HEOS_minus_T__constp, i, xN)) / (2 * dT);
1314 } else if (wrt == P_CONSTT) {
1315 numeric = (g(*HEOS_plus_p__constT, i, xN) - g(*HEOS_minus_p__constT, i, xN)) / (2 * dp);
1316 } else if (wrt == T_CONSTRHO) {
1317 double g1 = g(*HEOS_plus_T__constrho, i, xN), g2 = g(*HEOS_minus_T__constrho, i, xN);
1318 numeric = (g1 - g2) / (2 * dT);
1319 } else if (wrt == RHO_CONSTT) {
1320 numeric = (g(*HEOS_plus_rho__constT, i, xN) - g(*HEOS_minus_rho__constT, i, xN)) / (2 * drho);
1321 }
1322 CAPTURE(name);
1323 CAPTURE(i);
1324 CAPTURE(analytic);
1325 CAPTURE(numeric);
1326 CAPTURE(xN);
1327 double error = mix_deriv_err_func(numeric, analytic);
1328 CAPTURE(error);
1329 CHECK(error < tol);
1330 }
1331 };
1332 void one_comp(const std::string& name, one_mole_fraction_pointer f, zero_mole_fraction_pointer g, derivative wrt = CONST_TAU_DELTA) {
1333 for (int i = 0; i < 4; ++i) {
1334 double analytic;
1335 CHECK_NOTHROW(analytic = f(*HEOS, i, xN));
1336 double numeric = -10000;
1337 if (wrt == CONST_TAU_DELTA) {
1338 if (HEOS->get_mole_fractions()[i] > dz) {
1339 CHECK_NOTHROW(numeric = (g(*HEOS_plus_z[i], xN) - g(*HEOS_minus_z[i], xN)) / (2 * dz));
1340 } else {
1341 CHECK_NOTHROW(numeric = (-3 * g(*HEOS, xN) + 4 * g(*HEOS_plus_z[i], xN) - g(*HEOS_plus_2z[i], xN)) / (2 * dz));
1342 }
1343 } else if (wrt == CONST_TRHO) {
1344 CHECK_NOTHROW(numeric = (g(*HEOS_plus_z__constTrho[i], xN) - g(*HEOS_minus_z__constTrho[i], xN)) / (2 * dz));
1345 }
1346
1347 CAPTURE(name);
1348 CAPTURE(i);
1349 CAPTURE(analytic);
1350 CAPTURE(numeric);
1351 CAPTURE(xN);
1352 double error = mix_deriv_err_func(numeric, analytic);
1353 CAPTURE(error);
1354 CHECK(error < tol);
1355 }
1356 }
1357 void two(const std::string& name, two_mole_fraction_pointer f, two_mole_fraction_pointer g = NULL, derivative wrt = NO_DERIV) {
1358 for (int i = 0; i < 4; ++i) {
1359 for (int j = 0; j < 4; ++j) {
1360 double analytic = f(*HEOS, i, j, xN);
1361 bool is_other = false;
1362 double numeric = 0;
1363 if (wrt == TAU) {
1364 numeric = (g(*HEOS_plus_tau, i, j, xN) - g(*HEOS_minus_tau, i, j, xN)) / (2 * dtau);
1365 is_other = true;
1366 } else if (wrt == DELTA) {
1367 numeric = (g(*HEOS_plus_delta, i, j, xN) - g(*HEOS_minus_delta, i, j, xN)) / (2 * ddelta);
1368 is_other = true;
1369 }
1370 CAPTURE(name);
1371 CAPTURE(i);
1372 CAPTURE(j);
1373 CAPTURE(analytic);
1374 CAPTURE(numeric);
1375 CAPTURE(xN);
1376 double error = mix_deriv_err_func(numeric, analytic);
1377 CAPTURE(error);
1378 CHECK(error < tol);
1379 }
1380 }
1381 }
1382 void two_comp(const std::string& name, two_mole_fraction_pointer f, one_mole_fraction_pointer g, derivative wrt = NO_DERIV) {
1383 for (int i = 0; i < 4; ++i) {
1384 for (int j = 0; j < 4; ++j) {
1385 double analytic = f(*HEOS, i, j, xN);
1386 double numeric = 500;
1387 if (HEOS->get_mole_fractions()[j] > 2 * dz) {
1388 // Second order centered difference in composition
1389 CHECK_NOTHROW(numeric = (g(*HEOS_plus_z[j], i, xN) - g(*HEOS_minus_z[j], i, xN)) / (2 * dz));
1390 } else {
1391 // Forward difference in composition
1392 CHECK_NOTHROW(numeric = (-3 * g(*HEOS, i, xN) + 4 * g(*HEOS_plus_z[j], i, xN) - g(*HEOS_plus_2z[j], i, xN)) / (2 * dz));
1393 }
1394 CAPTURE(name);
1395 CAPTURE(i);
1396 CAPTURE(j);
1397 CAPTURE(analytic);
1398 CAPTURE(numeric);
1399 CAPTURE(xN);
1400 double error = mix_deriv_err_func(numeric, analytic);
1401 CAPTURE(error);
1402 CHECK(error < tol);
1403 }
1404 }
1405 }
1406 void three(const std::string& name, three_mole_fraction_pointer f, three_mole_fraction_pointer g = NULL, derivative wrt = NO_DERIV) {
1407 for (int i = 0; i < 4; ++i) {
1408 for (int j = 0; j < 4; ++j) {
1409 for (int k = 0; k < 4; ++k) {
1410 double analytic = f(*HEOS, i, j, k, xN);
1411 double numeric = 0;
1412 if (wrt == TAU) {
1413 numeric = (g(*HEOS_plus_tau, i, j, k, xN) - g(*HEOS_minus_tau, i, j, k, xN)) / (2 * dtau);
1414 } else if (wrt == DELTA) {
1415 numeric = (g(*HEOS_plus_delta, i, j, k, xN) - g(*HEOS_minus_delta, i, j, k, xN)) / (2 * ddelta);
1416 }
1417 CAPTURE(name);
1418 CAPTURE(i);
1419 CAPTURE(j);
1420 CAPTURE(analytic);
1421 CAPTURE(numeric);
1422 CAPTURE(xN);
1423 double error = mix_deriv_err_func(numeric, analytic);
1424 CAPTURE(error);
1425 CHECK(error < tol);
1426 }
1427 }
1428 }
1429 }
1430 void three_comp(const std::string& name, three_mole_fraction_pointer f, two_mole_fraction_pointer g, derivative wrt = NO_DERIV) {
1431 for (int i = 0; i < 4; ++i) {
1432 for (int j = 0; j < 4; ++j) {
1433 for (int k = 0; k < 4; ++k) {
1434 double analytic = f(*HEOS, i, j, k, xN);
1435 double numeric;
1436 if (HEOS->get_mole_fractions()[i] > 2 * dz) {
1437 CHECK_NOTHROW(numeric = (g(*HEOS_plus_z[k], i, j, xN) - g(*HEOS_minus_z[k], i, j, xN)) / (2 * dz));
1438 } else {
1439 CHECK_NOTHROW(numeric =
1440 (-3 * g(*HEOS, i, j, xN) + 4 * g(*HEOS_plus_z[k], i, j, xN) - g(*HEOS_plus_2z[k], i, j, xN)) / (2 * dz));
1441 }
1442 CAPTURE(name);
1443 CAPTURE(i);
1444 CAPTURE(j);
1445 CAPTURE(k);
1446 CAPTURE(analytic);
1447 CAPTURE(numeric);
1448 CAPTURE(xN);
1449 double error = mix_deriv_err_func(numeric, analytic);
1450 CAPTURE(error);
1451 CHECK(error < tol);
1452 }
1453 }
1454 }
1455 }
1456 Eigen::MatrixXd get_matrix(CoolProp::HelmholtzEOSMixtureBackend& HEOS, const std::string name) {
1457 Eigen::MatrixXd Lstar = MixtureDerivatives::Lstar(HEOS, xN);
1458 if (name == "Lstar") {
1459 return Lstar;
1460 } else if (name == "Mstar") {
1461 return MixtureDerivatives::Mstar(HEOS, xN, Lstar);
1462 } else {
1463 throw ValueError("Must be one of Lstar or Mstar");
1464 }
1465 }
1466 void matrix_derivative(const std::string name, const std::string& wrt) {
1467 CAPTURE(name);
1468 CAPTURE(wrt);
1469 double err = 10000;
1470 Eigen::MatrixXd analytic, numeric, Lstar, dLstar_dTau, dLstar_dDelta;
1471 if (name == "Mstar") {
1472 Lstar = MixtureDerivatives::Lstar(*HEOS, xN);
1473 dLstar_dDelta = MixtureDerivatives::dLstar_dX(*HEOS, xN, CoolProp::iDelta);
1474 dLstar_dTau = MixtureDerivatives::dLstar_dX(*HEOS, xN, CoolProp::iTau);
1475 }
1476 if (wrt == "tau") {
1477 if (name == "Lstar") {
1478 analytic = MixtureDerivatives::dLstar_dX(*HEOS, xN, CoolProp::iTau);
1479 } else if (name == "Mstar") {
1480 analytic = MixtureDerivatives::dMstar_dX(*HEOS, xN, CoolProp::iTau, Lstar, dLstar_dTau);
1481 } else {
1482 throw ValueError("argument not understood");
1483 }
1484 numeric = (get_matrix(*HEOS_plus_tau, name) - get_matrix(*HEOS_minus_tau, name)) / (2 * dtau);
1485 } else if (wrt == "delta") {
1486 if (name == "Lstar") {
1487 analytic = MixtureDerivatives::dLstar_dX(*HEOS, xN, CoolProp::iDelta);
1488 } else if (name == "Mstar") {
1489 analytic = MixtureDerivatives::dMstar_dX(*HEOS, xN, CoolProp::iDelta, Lstar, dLstar_dDelta);
1490 } else {
1491 throw ValueError("argument not understood");
1492 }
1493 numeric = (get_matrix(*HEOS_plus_delta, name) - get_matrix(*HEOS_minus_delta, name)) / (2 * ddelta);
1494 } else {
1495 throw ValueError("argument not understood");
1496 }
1497 CAPTURE(analytic);
1498 CAPTURE(numeric);
1499 Eigen::MatrixXd rel_error = ((analytic - numeric).cwiseQuotient(analytic));
1500 CAPTURE(rel_error);
1501 err = rel_error.squaredNorm();
1502 CAPTURE(err);
1503 CHECK(err < 1e-8);
1504 }
1505
1506 void run_checks() {
1507
1508 two_comp("d_PSI_rho_dxj", MD::d_PSI_rho_dxj, MD::PSI_rho);
1509 two_comp("d_PSI_T_dxj", MD::d_PSI_T_dxj, MD::PSI_T);
1510
1511 one("d_ndalphardni_dDelta", MD::d_ndalphardni_dDelta, MD::ndalphar_dni__constT_V_nj, DELTA);
1512 one("d2_ndalphardni_dDelta2", MD::d2_ndalphardni_dDelta2, MD::d_ndalphardni_dDelta, DELTA);
1513 one("d3_ndalphardni_dDelta3", MD::d3_ndalphardni_dDelta3, MD::d2_ndalphardni_dDelta2, DELTA);
1514 one("d_ndalphardni_dTau", MD::d_ndalphardni_dTau, MD::ndalphar_dni__constT_V_nj, TAU);
1515 one("d2_ndalphardni_dTau2", MD::d2_ndalphardni_dTau2, MD::d_ndalphardni_dTau, TAU);
1516 one("d3_ndalphardni_dTau3", MD::d3_ndalphardni_dTau3, MD::d2_ndalphardni_dTau2, TAU);
1517 one("d2_ndalphardni_dDelta_dTau", MD::d2_ndalphardni_dDelta_dTau, MD::d_ndalphardni_dDelta, TAU);
1518 one("d3_ndalphardni_dDelta2_dTau", MD::d3_ndalphardni_dDelta2_dTau, MD::d2_ndalphardni_dDelta2, TAU);
1519 one("d3_ndalphardni_dDelta_dTau2", MD::d3_ndalphardni_dDelta_dTau2, MD::d2_ndalphardni_dDelta_dTau, TAU);
1520
1521 zero("dalphar_dDelta", MD::dalphar_dDelta, MD::alphar, DELTA);
1522 zero("d2alphar_dDelta2", MD::d2alphar_dDelta2, MD::dalphar_dDelta, DELTA);
1523 zero("dalphar_dTau", MD::dalphar_dTau, MD::alphar, TAU);
1524 zero("d2alphar_dTau2", MD::d2alphar_dTau2, MD::dalphar_dTau, TAU);
1525
1526 zero("dalpha0_dDelta", MD::dalpha0_dDelta, MD::alpha0, DELTA);
1527 zero("d2alpha0_dDelta2", MD::d2alpha0_dDelta2, MD::dalpha0_dDelta, DELTA);
1528 zero("dalpha0_dTau", MD::dalpha0_dTau, MD::alpha0, TAU);
1529 zero("d2alpha0_dTau2", MD::d2alpha0_dTau2, MD::dalpha0_dTau, TAU);
1530
1531 one_comp("dalpha0_dxi", MD::dalpha0_dxi, MD::alpha0);
1532 one("d2alpha0_dxi_dDelta", MD::d2alpha0_dxi_dDelta, MD::dalpha0_dxi, DELTA);
1533 one("d2alpha0_dxi_dTau", MD::d2alpha0_dxi_dTau, MD::dalpha0_dxi, TAU);
1534 two_comp("d2alpha0dxidxj", MD::d2alpha0dxidxj, MD::dalpha0_dxi);
1535
1536 one_comp("dalpha_dxi", MD::dalpha_dxi, MD::alpha);
1537 one("d2alpha_dxi_dDelta", MD::d2alpha_dxi_dDelta, MD::dalpha_dxi, DELTA);
1538 one("d2alpha_dxi_dTau", MD::d2alpha_dxi_dTau, MD::dalpha_dxi, TAU);
1539 two_comp("d2alphadxidxj", MD::d2alphadxidxj, MD::dalpha_dxi);
1540
1541 zero("dpsi_dDelta", MD::dpsi_dDelta, MD::psi, DELTA);
1542 zero("dpsi_dTau", MD::dpsi_dTau, MD::psi, TAU);
1543 zero("d2psi_dDelta2", MD::d2psi_dDelta2, MD::dpsi_dDelta, DELTA);
1544 zero("d2psi_dDelta_dTau", MD::d2psi_dDelta_dTau, MD::dpsi_dDelta, TAU);
1545 zero("d2psi_dTau2", MD::d2psi_dTau2, MD::dpsi_dTau, TAU);
1546 one_comp("dpsi_dxi", MD::dpsi_dxi, MD::psi);
1547 one_comp("d2psi_dxi_dDelta", MD::d2psi_dxi_dDelta, MD::dpsi_dDelta);
1548 one_comp("d2psi_dxi_dTau", MD::d2psi_dxi_dTau, MD::dpsi_dTau);
1549 two_comp("d2psi_dxi_dxj", MD::d2psi_dxi_dxj, MD::dpsi_dxi);
1550
1551 //two_comp("d_ndalphardni_dxj__constT_V_xi", MD::d_ndalphardni_dxj__constT_V_xi, MD::ndalphar_dni__constT_V_nj);
1552
1553 one_comp("dalphar_dxi", MD::dalphar_dxi, MD::alphar);
1554 two_comp("d2alphardxidxj", MD::d2alphardxidxj, MD::dalphar_dxi);
1555 three_comp("d3alphardxidxjdxk", MD::d3alphardxidxjdxk, MD::d2alphardxidxj);
1556 one("d2alphar_dxi_dTau", MD::d2alphar_dxi_dTau, MD::dalphar_dxi, TAU);
1557 one("d2alphar_dxi_dDelta", MD::d2alphar_dxi_dDelta, MD::dalphar_dxi, DELTA);
1558 one("d3alphar_dxi_dDelta2", MD::d3alphar_dxi_dDelta2, MD::d2alphar_dxi_dDelta, DELTA);
1559 one("d3alphar_dxi_dTau2", MD::d3alphar_dxi_dTau2, MD::d2alphar_dxi_dTau, TAU);
1560 one("d4alphar_dxi_dTau3", MD::d4alphar_dxi_dTau3, MD::d3alphar_dxi_dTau2, TAU);
1561 one("d3alphar_dxi_dDelta_dTau", MD::d3alphar_dxi_dDelta_dTau, MD::d2alphar_dxi_dDelta, TAU);
1562 one("d4alphar_dxi_dDelta_dTau2", MD::d4alphar_dxi_dDelta_dTau2, MD::d3alphar_dxi_dDelta_dTau, TAU);
1563 one("d4alphar_dxi_dDelta2_dTau", MD::d4alphar_dxi_dDelta2_dTau, MD::d3alphar_dxi_dDelta2, TAU);
1564 two("d3alphar_dxi_dxj_dDelta", MD::d3alphar_dxi_dxj_dDelta, MD::d2alphardxidxj, DELTA);
1565 two("d4alphar_dxi_dxj_dDelta2", MD::d4alphar_dxi_dxj_dDelta2, MD::d3alphar_dxi_dxj_dDelta, DELTA);
1566 two("d4alphar_dxi_dxj_dDelta_dTau", MD::d4alphar_dxi_dxj_dDelta_dTau, MD::d3alphar_dxi_dxj_dDelta, TAU);
1567 two("d3alphar_dxi_dxj_dTau", MD::d3alphar_dxi_dxj_dTau, MD::d2alphardxidxj, TAU);
1568 two("d4alphar_dxi_dxj_dTau2", MD::d4alphar_dxi_dxj_dTau2, MD::d3alphar_dxi_dxj_dTau, TAU);
1569 one_comp("d_dalpharddelta_dxj__constT_V_xi", MD::d_dalpharddelta_dxj__constT_V_xi, MD::dalphar_dDelta, CONST_TRHO);
1570
1571 two_comp("d_ndalphardni_dxj__constdelta_tau_xi", MD::d_ndalphardni_dxj__constdelta_tau_xi, MD::ndalphar_dni__constT_V_nj);
1572 two("d2_ndalphardni_dxj_dDelta__consttau_xi", MD::d2_ndalphardni_dxj_dDelta__consttau_xi, MD::d_ndalphardni_dxj__constdelta_tau_xi, DELTA);
1573 two("d3_ndalphardni_dxj_dDelta2__consttau_xi", MD::d3_ndalphardni_dxj_dDelta2__consttau_xi, MD::d2_ndalphardni_dxj_dDelta__consttau_xi,
1574 DELTA);
1575 two("d2_ndalphardni_dxj_dTau__constdelta_xi", MD::d2_ndalphardni_dxj_dTau__constdelta_xi, MD::d_ndalphardni_dxj__constdelta_tau_xi, TAU);
1576 two("d3_ndalphardni_dxj_dTau2__constdelta_xi", MD::d3_ndalphardni_dxj_dTau2__constdelta_xi, MD::d2_ndalphardni_dxj_dTau__constdelta_xi, TAU);
1577 two("d3_ndalphardni_dxj_dDelta_dTau__constxi", MD::d3_ndalphardni_dxj_dDelta_dTau__constxi, MD::d2_ndalphardni_dxj_dDelta__consttau_xi, TAU);
1578 three_comp("d2_ndalphardni_dxj_dxk__constdelta_tau_xi", MD::d2_ndalphardni_dxj_dxk__constdelta_tau_xi,
1579 MD::d_ndalphardni_dxj__constdelta_tau_xi);
1580 three("d3_ndalphardni_dxj_dxk_dTau__constdelta_xi", MD::d3_ndalphardni_dxj_dxk_dTau__constdelta_xi,
1581 MD::d2_ndalphardni_dxj_dxk__constdelta_tau_xi, TAU);
1582 three("d3_ndalphardni_dxj_dxk_dDelta__consttau_xi", MD::d3_ndalphardni_dxj_dxk_dDelta__consttau_xi,
1583 MD::d2_ndalphardni_dxj_dxk__constdelta_tau_xi, DELTA);
1584
1585 // two("nd_ndalphardni_dnj__constT_V", MD::nd_ndalphardni_dnj__constT_V);
1586 two("d_nd_ndalphardni_dnj_dTau__constdelta_x", MD::d_nd_ndalphardni_dnj_dTau__constdelta_x, MD::nd_ndalphardni_dnj__constT_V, TAU);
1587 two("d2_nd_ndalphardni_dnj_dTau2__constdelta_x", MD::d2_nd_ndalphardni_dnj_dTau2__constdelta_x, MD::d_nd_ndalphardni_dnj_dTau__constdelta_x,
1588 TAU);
1589 two("d_nd_ndalphardni_dnj_dDelta__consttau_x", MD::d_nd_ndalphardni_dnj_dDelta__consttau_x, MD::nd_ndalphardni_dnj__constT_V, DELTA);
1590 two("d2_nd_ndalphardni_dnj_dDelta_dTau__constx", MD::d2_nd_ndalphardni_dnj_dDelta_dTau__constx, MD::d_nd_ndalphardni_dnj_dDelta__consttau_x,
1591 TAU);
1592 two("d2_nd_ndalphardni_dnj_dDelta2__consttau_x", MD::d2_nd_ndalphardni_dnj_dDelta2__consttau_x, MD::d_nd_ndalphardni_dnj_dDelta__consttau_x,
1593 DELTA);
1594 three_comp("d_nd_ndalphardni_dnj_dxk__consttau_delta", MD::d_nd_ndalphardni_dnj_dxk__consttau_delta, MD::nd_ndalphardni_dnj__constT_V);
1595 three("d2_nd_ndalphardni_dnj_dxk_dTau__constdelta", MD::d2_nd_ndalphardni_dnj_dxk_dTau__constdelta,
1596 MD::d_nd_ndalphardni_dnj_dxk__consttau_delta, TAU);
1597 three("d2_nd_ndalphardni_dnj_dxk_dDelta__consttau", MD::d2_nd_ndalphardni_dnj_dxk_dDelta__consttau,
1598 MD::d_nd_ndalphardni_dnj_dxk__consttau_delta, DELTA);
1599
1600 // Xn-dep only two_comp("dln_fugacity_dxj__constT_rho_xi", MD::dln_fugacity_dxj__constT_rho_xi, MD::ln_fugacity);
1601 three("d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau", MD::d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau,
1602 MD::d_ndln_fugacity_i_dnj_ddxk__consttau_delta, DELTA);
1603 three("d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta", MD::d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta,
1604 MD::d_ndln_fugacity_i_dnj_ddxk__consttau_delta, TAU);
1605 two("d2_ndln_fugacity_i_dnj_ddelta_dtau__constx", MD::d2_ndln_fugacity_i_dnj_ddelta_dtau__constx,
1606 MD::d_ndln_fugacity_i_dnj_ddelta__consttau_x, TAU);
1607 two("d_ndln_fugacity_i_dnj_ddelta__consttau_x", MD::d_ndln_fugacity_i_dnj_ddelta__consttau_x, MD::ndln_fugacity_i_dnj__constT_V_xi, DELTA);
1608 two("d_ndln_fugacity_i_dnj_dtau__constdelta_x", MD::d_ndln_fugacity_i_dnj_dtau__constdelta_x, MD::ndln_fugacity_i_dnj__constT_V_xi, TAU);
1609 three_comp("d_ndln_fugacity_i_dnj_ddxk__consttau_delta", MD::d_ndln_fugacity_i_dnj_ddxk__consttau_delta, MD::ndln_fugacity_i_dnj__constT_V_xi,
1610 TAU);
1611 one("dln_fugacity_i_dT__constrho_n", MD::dln_fugacity_i_dT__constrho_n, MD::ln_fugacity, T_CONSTRHO);
1612 one("dln_fugacity_i_drho__constT_n", MD::dln_fugacity_i_drho__constT_n, MD::ln_fugacity, RHO_CONSTT);
1613 one("dln_fugacity_i_dT__constp_n", MD::dln_fugacity_i_dT__constp_n, MD::ln_fugacity, T_CONSTP);
1614 one("dln_fugacity_i_dp__constT_n", MD::dln_fugacity_i_dp__constT_n, MD::ln_fugacity, P_CONSTT);
1615 one("dln_fugacity_coefficient_dT__constp_n", MD::dln_fugacity_coefficient_dT__constp_n, MD::ln_fugacity_coefficient, T_CONSTP);
1616 one("dln_fugacity_coefficient_dp__constT_n", MD::dln_fugacity_coefficient_dp__constT_n, MD::ln_fugacity_coefficient, P_CONSTT);
1617
1618 three_comp("d2_PSI_T_dxj_dxk", MD::d2_PSI_T_dxj_dxk, MD::d_PSI_T_dxj);
1619 three_comp("d2_PSI_rho_dxj_dxk", MD::d2_PSI_rho_dxj_dxk, MD::d_PSI_rho_dxj);
1620
1621 three("d_n2Aijk_dTau", MD::d_n2Aijk_dTau, MD::n2Aijk, TAU);
1622 three("d_n2Aijk_dDelta", MD::d_n2Aijk_dDelta, MD::n2Aijk, DELTA);
1623 two("d_nAij_dTau", MD::d_nAij_dTau, MD::nAij, TAU);
1624 two("d_nAij_dDelta", MD::d_nAij_dDelta, MD::nAij, DELTA);
1625
1626 two_comp("d_nddeltadni_dxj__constdelta_tau", MD::d_nddeltadni_dxj__constdelta_tau, MD::nddeltadni__constT_V_nj);
1627 two_comp("d_ndtaudni_dxj__constdelta_tau", MD::d_ndtaudni_dxj__constdelta_tau, MD::ndtaudni__constT_V_nj);
1628 two("d2_ndtaudni_dxj_dTau__constdelta", MD::d2_ndtaudni_dxj_dTau__constdelta, MD::d_ndtaudni_dxj__constdelta_tau, TAU);
1629
1630 one_comp("dTrdxi__constxj", MD::dTrdxi__constxj, MD::Tr);
1631 one_comp("d2Tr_dxidbetaT", MD::d2Tr_dxidbetaT, MD::dTr_dbetaT);
1632 one_comp("d2Tr_dxidgammaT", MD::d2Tr_dxidgammaT, MD::dTr_dgammaT);
1633 // (??)two_comp("d2Trdxi2__constxj", MD::d2Trdxi2__constxj, MD::dTrdxi__constxj);
1634 two_comp("d2Trdxidxj", MD::d2Trdxidxj, MD::dTrdxi__constxj);
1635 three_comp("d3Trdxidxjdxk", MD::d3Trdxidxjdxk, MD::d2Trdxidxj);
1636 one_comp("dtaudxj__constT_V_xi", MD::dtau_dxj__constT_V_xi, MD::tau, CONST_TRHO);
1637 two_comp("d_ndTrdni_dxj__constxi", MD::d_ndTrdni_dxj__constxi, MD::ndTrdni__constnj);
1638 three_comp("d2_ndTrdni_dxj_dxk__constxi", MD::d2_ndTrdni_dxj_dxk__constxi, MD::d_ndTrdni_dxj__constxi);
1639
1640 one_comp("drhormolardxi__constxj", MD::drhormolardxi__constxj, MD::rhormolar);
1641 one_comp("d2rhormolar_dxidbetaV", MD::d2rhormolar_dxidbetaV, MD::drhormolar_dbetaV);
1642 one_comp("d2rhormolar_dxidgammaV", MD::d2rhormolar_dxidgammaV, MD::drhormolar_dgammaV);
1643 // (??) two_comp("d2rhormolardxi2__constxj", MD::d2rhormolardxi2__constxj, MD::drhormolardxi__constxj);
1644 two_comp("d2rhormolardxidxj", MD::d2rhormolardxidxj, MD::drhormolardxi__constxj);
1645 three_comp("d3rhormolardxidxjdxk", MD::d3rhormolardxidxjdxk, MD::d2rhormolardxidxj);
1646 one_comp("ddelta_dxj__constT_V_xi", MD::ddelta_dxj__constT_V_xi, MD::delta, CONST_TRHO);
1647 two_comp("d_ndrhorbardni_dxj__constxi", MD::d_ndrhorbardni_dxj__constxi, MD::ndrhorbardni__constnj);
1648 three_comp("d2_ndrhorbardni_dxj_dxk__constxi", MD::d2_ndrhorbardni_dxj_dxk__constxi, MD::d_ndrhorbardni_dxj__constxi);
1649
1650 one_comp("dpdxj__constT_V_xi", MD::dpdxj__constT_V_xi, MD::p, CONST_TRHO);
1651
1652 matrix_derivative("Lstar", "tau");
1653 matrix_derivative("Lstar", "delta");
1654 matrix_derivative("Mstar", "tau");
1655 matrix_derivative("Mstar", "delta");
1656 }
1657};
1658
1659TEST_CASE_METHOD(DerivativeFixture<HelmholtzEOSMixtureBackend>, "Check derivatives for HEOS", "[mixture_derivs2]") {
1660 run_checks();
1661};
1662TEST_CASE_METHOD(DerivativeFixture<PengRobinsonBackend>, "Check derivatives for Peng-Robinson", "[mixture_derivs2]") {
1663 tol = 1e-4; // Relax the tolerance a bit
1664 run_checks();
1665};
1666TEST_CASE_METHOD(DerivativeFixture<SRKBackend>, "Check derivatives for SRK", "[mixture_derivs2]") {
1667 tol = 1e-4; // Relax the tolerance a bit
1668 run_checks();
1669};
1670 // Make sure you set the VTPR UNIFAC path with something like set_config_string(VTPR_UNIFAC_PATH, "/Users/ian/Code/CUBAC/dev/unifaq/");
1671 //TEST_CASE_METHOD(DerivativeFixture<VTPRBackend>, "Check derivatives for VTPR", "[mixture_derivs2]")
1672 //{
1673 // tol = 1e-4; // Relax the tolerance a bit
1674 // run_checks();
1675 //};
1676
1677#endif