CoolProp 8.0.0
An open-source fluid property and humid air property database
TransportRoutines.cpp
Go to the documentation of this file.
1
2#include "TransportRoutines.h"
3
4#include <cmath>
5#include "CoolProp/CoolPropFluid.h"
6
7namespace CoolProp {
8
9CoolPropDbl TransportRoutines::viscosity_dilute_kinetic_theory(HelmholtzEOSMixtureBackend& HEOS) {
10 if (HEOS.is_pure_or_pseudopure) {
11 CoolPropDbl Tstar = HEOS.T() / HEOS.components[0].transport.epsilon_over_k;
12 CoolPropDbl sigma_nm = HEOS.components[0].transport.sigma_eta * 1e9; // 1e9 to convert from m to nm
13 CoolPropDbl molar_mass_kgkmol = HEOS.molar_mass() * 1000; // 1000 to convert from kg/mol to kg/kmol
14
15 // The nondimensional empirical collision integral from Neufeld
16 // Neufeld, P. D.; Janzen, A. R.; Aziz, R. A. Empirical Equations to Calculate 16 of the Transport Collision Integrals (l,s)*
17 // for the Lennard-Jones (12-6) Potential. J. Chem. Phys. 1972, 57, 1100-1102
18 CoolPropDbl OMEGA22 =
19 1.16145 * pow(Tstar, static_cast<CoolPropDbl>(-0.14874)) + 0.52487 * exp(-0.77320 * Tstar) + 2.16178 * exp(-2.43787 * Tstar);
20
21 // The dilute gas component -
22 return 26.692e-9 * sqrt(molar_mass_kgkmol * HEOS.T()) / (pow(sigma_nm, 2) * OMEGA22); // Pa-s
23 } else {
24 throw NotImplementedError("TransportRoutines::viscosity_dilute_kinetic_theory is only for pure and pseudo-pure");
25 }
26}
27
28CoolPropDbl TransportRoutines::viscosity_dilute_collision_integral(HelmholtzEOSMixtureBackend& HEOS) {
29 if (HEOS.is_pure_or_pseudopure) {
30 // Retrieve values from the state class
31 CoolProp::ViscosityDiluteGasCollisionIntegralData& data = HEOS.components[0].transport.viscosity_dilute.collision_integral;
32 const std::vector<CoolPropDbl>&a = data.a, &t = data.t;
33 const CoolPropDbl C = data.C, molar_mass = data.molar_mass;
34
35 CoolPropDbl S = NAN;
36 // Unit conversions and variable definitions
37 const CoolPropDbl Tstar = HEOS.T() / HEOS.components[0].transport.epsilon_over_k;
38 const CoolPropDbl sigma_nm = HEOS.components[0].transport.sigma_eta * 1e9; // 1e9 to convert from m to nm
39 const CoolPropDbl molar_mass_kgkmol = molar_mass * 1000; // 1000 to convert from kg/mol to kg/kmol
40
44 CoolPropDbl summer = 0, lnTstar = log(Tstar);
45 for (std::size_t i = 0; i < a.size(); ++i) {
46 summer += a[i] * pow(lnTstar, t[i]);
47 }
48 S = exp(summer);
49
50 // The dilute gas component
51 return C * sqrt(molar_mass_kgkmol * HEOS.T()) / (pow(sigma_nm, 2) * S); // Pa-s
52 } else {
53 throw NotImplementedError("TransportRoutines::viscosity_dilute_collision_integral is only for pure and pseudo-pure");
54 }
55}
56
57CoolPropDbl TransportRoutines::viscosity_dilute_powers_of_T(HelmholtzEOSMixtureBackend& HEOS) {
58 if (HEOS.is_pure_or_pseudopure) {
59 // Retrieve values from the state class
60 CoolProp::ViscosityDiluteGasPowersOfT& data = HEOS.components[0].transport.viscosity_dilute.powers_of_T;
61 const std::vector<CoolPropDbl>&a = data.a, &t = data.t;
62
63 CoolPropDbl summer = 0, T = HEOS.T();
64 for (std::size_t i = 0; i < a.size(); ++i) {
65 summer += a[i] * pow(T, t[i]);
66 }
67 return summer;
68 } else {
69 throw NotImplementedError("TransportRoutines::viscosity_dilute_powers_of_T is only for pure and pseudo-pure");
70 }
71}
72CoolPropDbl TransportRoutines::viscosity_dilute_powers_of_Tr(HelmholtzEOSMixtureBackend& HEOS) {
73 if (HEOS.is_pure_or_pseudopure) {
74 // Retrieve values from the state class
75 CoolProp::ViscosityDiluteGasPowersOfTr& data = HEOS.components[0].transport.viscosity_dilute.powers_of_Tr;
76 const std::vector<CoolPropDbl>&a = data.a, &t = data.t;
77 CoolPropDbl summer = 0, Tr = HEOS.T() / data.T_reducing;
78 for (std::size_t i = 0; i < a.size(); ++i) {
79 summer += a[i] * pow(Tr, t[i]);
80 }
81 return summer;
82 } else {
83 throw NotImplementedError("TransportRoutines::viscosity_dilute_powers_of_Tr is only for pure and pseudo-pure");
84 }
85}
86
87CoolPropDbl TransportRoutines::viscosity_dilute_collision_integral_powers_of_T(HelmholtzEOSMixtureBackend& HEOS) {
88 if (HEOS.is_pure_or_pseudopure) {
89 // Retrieve values from the state class
90 CoolProp::ViscosityDiluteCollisionIntegralPowersOfTstarData& data =
91 HEOS.components[0].transport.viscosity_dilute.collision_integral_powers_of_Tstar;
92 const std::vector<CoolPropDbl>&a = data.a, &t = data.t;
93
94 CoolPropDbl summer = 0, Tstar = HEOS.T() / data.T_reducing;
95 for (std::size_t i = 0; i < a.size(); ++i) {
96 summer += a[i] * pow(Tstar, t[i]);
97 }
98 return data.C * sqrt(HEOS.T()) / summer;
99 } else {
100 throw NotImplementedError("TransportRoutines::viscosity_dilute_collision_integral_powers_of_T is only for pure and pseudo-pure");
101 }
102}
103CoolPropDbl TransportRoutines::viscosity_higher_order_modified_Batschinski_Hildebrand(HelmholtzEOSMixtureBackend& HEOS) {
104 if (HEOS.is_pure_or_pseudopure) {
105 CoolProp::ViscosityModifiedBatschinskiHildebrandData& HO =
106 HEOS.components[0].transport.viscosity_higher_order.modified_Batschinski_Hildebrand;
107
108 CoolPropDbl delta = HEOS.rhomolar() / HO.rhomolar_reduce, tau = HO.T_reduce / HEOS.T();
109
110 // The first term that is formed of powers of tau (Tc/T) and delta (rho/rhoc)
111 CoolPropDbl S = 0;
112 for (unsigned int i = 0; i < HO.a.size(); ++i) {
113 S += HO.a[i] * pow(delta, HO.d1[i]) * pow(tau, HO.t1[i]) * exp(HO.gamma[i] * pow(delta, HO.l[i]));
114 }
115
116 // For the terms that multiplies the bracketed term with delta and delta0
117 CoolPropDbl F = 0;
118 for (unsigned int i = 0; i < HO.f.size(); ++i) {
119 F += HO.f[i] * pow(delta, HO.d2[i]) * pow(tau, HO.t2[i]);
120 }
121
122 // for delta_0
123 CoolPropDbl summer_numer = 0;
124 for (unsigned int i = 0; i < HO.g.size(); ++i) {
125 summer_numer += HO.g[i] * pow(tau, HO.h[i]);
126 }
127 CoolPropDbl summer_denom = 0;
128 for (unsigned int i = 0; i < HO.p.size(); ++i) {
129 summer_denom += HO.p[i] * pow(tau, HO.q[i]);
130 }
131 CoolPropDbl delta0 = summer_numer / summer_denom;
132
133 // The higher-order-term component
134 return S + F * (1 / (delta0 - delta) - 1 / delta0); // Pa-s
135 } else {
136 throw NotImplementedError("TransportRoutines::viscosity_higher_order_modified_Batschinski_Hildebrand is only for pure and pseudo-pure");
137 }
138}
139
140CoolPropDbl TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend(HelmholtzEOSMixtureBackend& HEOS) {
141 if (HEOS.is_pure_or_pseudopure) {
142 // Retrieve values from the state class
143 CoolProp::ViscosityRainWaterFriendData& data = HEOS.components[0].transport.viscosity_initial.rainwater_friend;
144 const std::vector<CoolPropDbl>&b = data.b, &t = data.t;
145
146 CoolPropDbl B_eta = NAN, B_eta_star = NAN;
147 CoolPropDbl Tstar = HEOS.T() / HEOS.components[0].transport.epsilon_over_k; // [no units]
148 CoolPropDbl sigma = HEOS.components[0].transport.sigma_eta; // [m]
149
150 CoolPropDbl summer = 0;
151 for (unsigned int i = 0; i < b.size(); ++i) {
152 summer += b[i] * pow(Tstar, t[i]);
153 }
154 B_eta_star = summer; // [no units]
155 B_eta = 6.02214129e23 * pow(sigma, 3) * B_eta_star; // [m^3/mol]
156 return B_eta; // [m^3/mol]
157 } else {
158 throw NotImplementedError("TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend is only for pure and pseudo-pure");
159 }
160}
161
162CoolPropDbl TransportRoutines::viscosity_initial_density_dependence_empirical(HelmholtzEOSMixtureBackend& HEOS) {
163 // Inspired by the form from Tariq, JPCRD, 2014
164 if (HEOS.is_pure_or_pseudopure) {
165 // Retrieve values from the state class
166 CoolProp::ViscosityInitialDensityEmpiricalData& data = HEOS.components[0].transport.viscosity_initial.empirical;
167 const std::vector<CoolPropDbl>&n = data.n, &d = data.d, &t = data.t;
168
169 CoolPropDbl tau = data.T_reducing / HEOS.T(); // [no units]
170 CoolPropDbl delta = HEOS.rhomolar() / data.rhomolar_reducing; // [no units]
171
172 CoolPropDbl summer = 0;
173 for (unsigned int i = 0; i < n.size(); ++i) {
174 summer += n[i] * pow(delta, d[i]) * pow(tau, t[i]);
175 }
176 return summer; // [Pa-s]
177 } else {
178 throw NotImplementedError("TransportRoutines::viscosity_initial_density_dependence_empirical is only for pure and pseudo-pure");
179 }
180}
181
182static void visc_Helper(double Tbar, double rhobar, double* mubar_0, double* mubar_1) {
183 std::vector<std::vector<CoolPropDbl>> H(6, std::vector<CoolPropDbl>(7, 0));
184 double sum = NAN;
185 int i = 0, j = 0;
186
187 // Dilute-gas component
188 *mubar_0 = 100.0 * sqrt(Tbar) / (1.67752 + 2.20462 / Tbar + 0.6366564 / powInt(Tbar, 2) - 0.241605 / powInt(Tbar, 3));
189
190 //Fill in zeros in H
191 for (i = 0; i <= 5; i++) {
192 for (j = 0; j <= 6; j++) {
193 H[i][j] = 0;
194 }
195 }
196
197 //Set non-zero parameters of H
198 H[0][0] = 5.20094e-1;
199 H[1][0] = 8.50895e-2;
200 H[2][0] = -1.08374;
201 H[3][0] = -2.89555e-1;
202
203 H[0][1] = 2.22531e-1;
204 H[1][1] = 9.99115e-1;
205 H[2][1] = 1.88797;
206 H[3][1] = 1.26613;
207 H[5][1] = 1.20573e-1;
208
209 H[0][2] = -2.81378e-1;
210 H[1][2] = -9.06851e-1;
211 H[2][2] = -7.72479e-1;
212 H[3][2] = -4.89837e-1;
213 H[4][2] = -2.57040e-1;
214
215 H[0][3] = 1.61913e-1;
216 H[1][3] = 2.57399e-1;
217
218 H[0][4] = -3.25372e-2;
219 H[3][4] = 6.98452e-2;
220
221 H[4][5] = 8.72102e-3;
222
223 H[3][6] = -4.35673e-3;
224 H[5][6] = -5.93264e-4;
225
226 // Finite density component
227 sum = 0;
228 for (i = 0; i <= 5; i++) {
229 for (j = 0; j <= 6; j++) {
230 sum += powInt(1 / Tbar - 1, i) * (H[i][j] * powInt(rhobar - 1, j));
231 }
232 }
233 *mubar_1 = exp(rhobar * sum);
234}
235CoolPropDbl TransportRoutines::viscosity_heavywater_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
236 double Tbar = HEOS.T() / 643.847, rhobar = HEOS.rhomass() / 358;
237 double A[] = {1.000000, 0.940695, 0.578377, -0.202044};
238 int I[] = {0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 0, 1, 2, 5, 0, 1, 2, 3, 0, 1, 3, 5, 0, 1, 5, 3};
239 int J[] = {0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6};
240 double Bij[] = {0.4864192, -0.2448372, -0.8702035, 0.8716056, -1.051126, 0.3458395, 0.3509007, 1.315436, 1.297752,
241 1.353448, -0.2847572, -1.037026, -1.287846, -0.02148229, 0.07013759, 0.4660127, 0.2292075, -0.4857462,
242 0.01641220, -0.02884911, 0.1607171, -0.009603846, -0.01163815, -0.008239587, 0.004559914, -0.003886659};
243 double mu0 = sqrt(Tbar) / (A[0] + A[1] / Tbar + A[2] / POW2(Tbar) + A[3] / POW3(Tbar));
244 double summer = 0;
245 for (int i = 0; i < 26; ++i) {
246 summer += Bij[i] * pow(1 / Tbar - 1, I[i]) * pow(rhobar - 1, J[i]);
247 }
248 double mu1 = exp(rhobar * summer);
249 double mubar = mu0 * mu1;
250 return 55.2651e-6 * mubar;
251}
252CoolPropDbl TransportRoutines::viscosity_water_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
253 double x_mu = 0.068, qc = 1 / 1.9, qd = 1 / 1.1, nu = 0.630, gamma = 1.239, zeta_0 = 0.13, LAMBDA_0 = 0.06, Tbar_R = 1.5, pstar = NAN,
254 Tstar = NAN, rhostar = NAN;
255 double delta = NAN, tau = NAN, mubar_0 = NAN, mubar_1 = NAN, mubar_2 = NAN, drhodp = NAN, drhodp_R = NAN, DeltaChibar = NAN, zeta = NAN, w = NAN,
256 L = NAN, Y = NAN, psi_D = NAN, Tbar = NAN, rhobar = NAN;
257 double drhobar_dpbar = NAN, drhobar_dpbar_R = NAN, R_Water = NAN;
258
259 pstar = 22.064e6; // [Pa]
260 Tstar = 647.096; // [K]
261 rhostar = 322; // [kg/m^3]
262 Tbar = HEOS.T() / Tstar;
263 rhobar = HEOS.rhomass() / rhostar;
264 R_Water = HEOS.gas_constant() / HEOS.molar_mass(); // [J/kg/K]
265
266 // Dilute and finite gas portions
267 visc_Helper(Tbar, rhobar, &mubar_0, &mubar_1);
268
269 // **********************************************************************
270 // ************************ Critical Enhancement ************************
271 // **********************************************************************
272 delta = rhobar;
273 // "Normal" calculation
274 drhodp = 1 / (R_Water * HEOS.T() * (1 + 2 * delta * HEOS.dalphar_dDelta() + delta * delta * HEOS.d2alphar_dDelta2()));
275 drhobar_dpbar = pstar / rhostar * drhodp;
276 // "Reducing" calculation
277 tau = 1 / Tbar_R;
278 drhodp_R = 1
279 / (R_Water * Tbar_R * Tstar
280 * (1 + 2 * rhobar * HEOS.calc_alphar_deriv_nocache(0, 1, HEOS.mole_fractions, tau, delta)
281 + delta * delta * HEOS.calc_alphar_deriv_nocache(0, 2, HEOS.mole_fractions, tau, delta)));
282 drhobar_dpbar_R = pstar / rhostar * drhodp_R;
283
284 DeltaChibar = rhobar * (drhobar_dpbar - drhobar_dpbar_R * Tbar_R / Tbar);
285 if (DeltaChibar < 0) DeltaChibar = 0;
286 zeta = zeta_0 * pow(DeltaChibar / LAMBDA_0, nu / gamma);
287 if (zeta < 0.3817016416) {
288 Y = 1.0 / 5.0 * qc * zeta * powInt(qd * zeta, 5) * (1 - qc * zeta + powInt(qc * zeta, 2) - 765.0 / 504.0 * powInt(qd * zeta, 2));
289 } else {
290 psi_D = acos(pow(1 + powInt(qd * zeta, 2), -1.0 / 2.0));
291 w = sqrt(std::abs((qc * zeta - 1) / (qc * zeta + 1))) * tan(psi_D / 2.0);
292 if (qc * zeta > 1) {
293 L = log((1 + w) / (1 - w));
294 } else {
295 L = 2 * atan(std::abs(w));
296 }
297 Y = 1.0 / 12.0 * sin(3 * psi_D) - 1 / (4 * qc * zeta) * sin(2 * psi_D)
298 + 1.0 / powInt(qc * zeta, 2) * (1 - 5.0 / 4.0 * powInt(qc * zeta, 2)) * sin(psi_D)
299 - 1.0 / powInt(qc * zeta, 3) * ((1 - 3.0 / 2.0 * powInt(qc * zeta, 2)) * psi_D - pow(std::abs(powInt(qc * zeta, 2) - 1), 3.0 / 2.0) * L);
300 }
301 mubar_2 = exp(x_mu * Y);
302
303 return (mubar_0 * mubar_1 * mubar_2) / 1e6;
304}
305CoolPropDbl TransportRoutines::viscosity_toluene_higher_order_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
306 CoolPropDbl Tr = HEOS.T() / 591.75, rhor = HEOS.keyed_output(CoolProp::iDmass) / 291.987;
307 CoolPropDbl c[] = {19.919216, -2.6557905, -135.904211, -7.9962719, -11.014795, -10.113817};
308 return 1e-6 * pow(static_cast<double>(rhor), 2.0 / 3.0) * sqrt(Tr)
309 * ((c[0] * rhor + c[1] * pow(rhor, 4)) / Tr + c[2] * rhor * rhor * rhor / (rhor * rhor + c[3] + c[4] * Tr) + c[5] * rhor);
310}
311
312CoolPropDbl TransportRoutines::viscosity_hydrogen_higher_order_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
313 CoolPropDbl Tr = HEOS.T() / 33.145, rhor = HEOS.keyed_output(CoolProp::iDmass) * 0.011;
314 CoolPropDbl c[] = {0, 6.43449673e-6, 4.56334068e-2, 2.32797868e-1, 9.58326120e-1, 1.27941189e-1, 3.63576595e-1};
315 return c[1] * pow(rhor, 2) * exp(c[2] * Tr + c[3] / Tr + c[4] * pow(rhor, 2) / (c[5] + Tr) + c[6] * pow(rhor, 6));
316}
317CoolPropDbl TransportRoutines::viscosity_benzene_higher_order_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
318 CoolPropDbl Tr = HEOS.T() / 562.02, rhor = HEOS.rhomass() / 304.792;
319 CoolPropDbl c[] = {-9.98945, 86.06260, 2.74872, 1.11130, -1.0, -134.1330, -352.473, 6.60989, 88.4174};
320 return 1e-6 * pow(rhor, static_cast<CoolPropDbl>(2.0 / 3.0)) * sqrt(Tr)
321 * (c[0] * pow(rhor, 2) + c[1] * rhor / (c[2] + c[3] * Tr + c[4] * rhor)
322 + (c[5] * rhor + c[6] * pow(rhor, 2)) / (c[7] + c[8] * pow(rhor, 2)));
323}
324CoolPropDbl TransportRoutines::viscosity_hexane_higher_order_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
325
326 CoolPropDbl Tr = HEOS.T() / 507.82, rhor = HEOS.keyed_output(CoolProp::iDmass) / 233.182;
327
328 // Output is in Pa-s
329 double c[] = {2.53402335 / 1e6, -9.724061002 / 1e6, 0.469437316, 158.5571631, 72.42916856 / 1e6,
330 10.60751253, 8.628373915, -6.61346441, -2.212724566};
331 return pow(rhor, static_cast<CoolPropDbl>(2.0 / 3.0)) * sqrt(Tr)
332 * (c[0] / Tr + c[1] / (c[2] + Tr + c[3] * rhor * rhor)
333 + c[4] * (1 + rhor) / (c[5] + c[6] * Tr + c[7] * rhor + rhor * rhor + c[8] * rhor * Tr));
334}
335
336CoolPropDbl TransportRoutines::viscosity_heptane_higher_order_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
338 CoolPropDbl Tr = HEOS.T() / 540.13, rhor = HEOS.rhomass() / 232;
339
340 // Output is in Pa-s
341 double c[] = {0, 22.15000 / 1e6, -15.00870 / 1e6, 3.71791 / 1e6, 77.72818 / 1e6, 9.73449, 9.51900, -6.34076, -2.51909};
342 return pow(rhor, static_cast<CoolPropDbl>(2.0L / 3.0L)) * sqrt(Tr)
343 * (c[1] * rhor + c[2] * pow(rhor, 2) + c[3] * pow(rhor, 3)
344 + c[4] * rhor / (c[5] + c[6] * Tr + c[7] * rhor + rhor * rhor + c[8] * rhor * Tr));
345}
346
347CoolPropDbl TransportRoutines::viscosity_CO2_higher_order_hardcoded_LaeseckeJPCRD2017(HelmholtzEOSMixtureBackend& HEOS) {
348 double c1 = 0.360603235428487, c2 = 0.121550806591497, gamma = 8.06282737481277;
349 double Tt = HEOS.Ttriple(), rho_tL = 1178.53;
350 double Tr = HEOS.T() / Tt, rhor = HEOS.rhomass() / rho_tL;
351 // Eq. (9) from Laesecke, JPCRD, 2017
352 double eta_tL = pow(rho_tL, 2.0 / 3.0) * sqrt(HEOS.gas_constant() * Tt) / (pow(HEOS.molar_mass(), 1.0 / 6.0) * 84446887.43579945);
353 // Eq. (8) from Laesecke, JPCRD, 2017
354 double residual = eta_tL * (c1 * Tr * pow(rhor, 3) + (pow(rhor, 2) + pow(rhor, gamma)) / (Tr - c2));
355 return residual;
356}
357
358CoolPropDbl TransportRoutines::viscosity_higher_order_friction_theory(HelmholtzEOSMixtureBackend& HEOS) {
359 if (HEOS.is_pure_or_pseudopure) {
360 CoolProp::ViscosityFrictionTheoryData& F = HEOS.components[0].transport.viscosity_higher_order.friction_theory;
361
362 CoolPropDbl tau = F.T_reduce / HEOS.T(), kii = 0, krrr = 0, kaaa = 0, krr = NAN, kdrdr = NAN;
363
364 double psi1 = exp(tau) - F.c1;
365 double psi2 = exp(pow(tau, 2)) - F.c2;
366
367 double ki = (F.Ai[0] + F.Ai[1] * psi1 + F.Ai[2] * psi2) * tau;
368
369 double ka = (F.Aa[0] + F.Aa[1] * psi1 + F.Aa[2] * psi2) * pow(tau, F.Na);
370 double kr = (F.Ar[0] + F.Ar[1] * psi1 + F.Ar[2] * psi2) * pow(tau, F.Nr);
371 double kaa = (F.Aaa[0] + F.Aaa[1] * psi1 + F.Aaa[2] * psi2) * pow(tau, F.Naa);
372 if (F.Arr.empty()) {
373 krr = 0;
374 kdrdr = (F.Adrdr[0] + F.Adrdr[1] * psi1 + F.Adrdr[2] * psi2) * pow(tau, F.Nrr);
375 } else {
376 krr = (F.Arr[0] + F.Arr[1] * psi1 + F.Arr[2] * psi2) * pow(tau, F.Nrr);
377 kdrdr = 0;
378 }
379 if (!F.Aii.empty()) {
380 kii = (F.Aii[0] + F.Aii[1] * psi1 + F.Aii[2] * psi2) * pow(tau, F.Nii);
381 }
382 if (!F.Arrr.empty() && !F.Aaaa.empty()) {
383 krrr = (F.Arrr[0] + F.Arrr[1] * psi1 + F.Arrr[2] * psi2) * pow(tau, F.Nrrr);
384 kaaa = (F.Aaaa[0] + F.Aaaa[1] * psi1 + F.Aaaa[2] * psi2) * pow(tau, F.Naaa);
385 }
386
387 double p = HEOS.p() / 1e5; // [bar]; 1e5 for conversion from Pa -> bar
388 double pr =
389 HEOS.T() * HEOS.first_partial_deriv(CoolProp::iP, CoolProp::iT, CoolProp::iDmolar) / 1e5; // [bar/K]; 1e5 for conversion from Pa -> bar
390 double pa = p - pr; //[bar]
391 double pid = HEOS.rhomolar() * HEOS.gas_constant() * HEOS.T() / 1e5; // [bar]; 1e5 for conversion from Pa -> bar
392 double deltapr = pr - pid;
393
394 double eta_f = ka * pa + kr * deltapr + ki * pid + kaa * pa * pa + kdrdr * deltapr * deltapr + krr * pr * pr + kii * pid * pid
395 + krrr * pr * pr * pr + kaaa * pa * pa * pa;
396
397 return eta_f; //[Pa-s]
398 } else {
399 throw NotImplementedError("TransportRoutines::viscosity_higher_order_friction_theory is only for pure and pseudo-pure");
400 }
401}
402
403CoolPropDbl TransportRoutines::viscosity_helium_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
404 double eta_0 = NAN, eta_0_slash = NAN, eta_E_slash = NAN, B = NAN, C = NAN, D = NAN, ln_eta = NAN, x = NAN;
405 //
406 // Arp, V.D., McCarty, R.D., and Friend, D.G.,
407 // "Thermophysical Properties of Helium-4 from 0.8 to 1500 K with Pressures to 2000 MPa",
408 // NIST Technical Note 1334 (revised), 1998.
409 //
410 // Using Arp NIST report
411 // Report is not clear on viscosity, referring to REFPROP source code for clarity
412
413 // Correlation wants density in g/cm^3; kg/m^3 --> g/cm^3, divide by 1000
414 CoolPropDbl rho = HEOS.keyed_output(CoolProp::iDmass) / 1000.0, T = HEOS.T();
415
416 if (T <= 300) {
417 x = log(T);
418 } else {
419 x = log(300.0);
420 }
421 // Evaluate the terms B,C,D
422 B = -47.5295259 / x + 87.6799309 - 42.0741589 * x + 8.33128289 * x * x - 0.589252385 * x * x * x;
423 C = 547.309267 / x - 904.870586 + 431.404928 * x - 81.4504854 * x * x + 5.37008433 * x * x * x;
424 D = -1684.39324 / x + 3331.08630 - 1632.19172 * x + 308.804413 * x * x - 20.2936367 * x * x * x;
425 eta_0_slash = -0.135311743 / x + 1.00347841 + 1.20654649 * x - 0.149564551 * x * x + 0.012520841 * x * x * x;
426 eta_E_slash = rho * B + rho * rho * C + rho * rho * rho * D;
427
428 if (T <= 100) {
429 ln_eta = eta_0_slash + eta_E_slash;
430 // Correlation yields viscosity in micro g/(cm-s); to get Pa-s, divide by 10 to get micro Pa-s, then another 1e6 to get Pa-s
431 return exp(ln_eta) / 10.0 / 1e6;
432 } else {
433 ln_eta = eta_0_slash + eta_E_slash;
434 eta_0 = 196 * pow(T, static_cast<CoolPropDbl>(0.71938)) * exp(12.451 / T - 295.67 / T / T - 4.1249);
435 // Correlation yields viscosity in micro g/(cm-s); to get Pa-s, divide by 10 to get micro Pa-s, then another 1e6 to get Pa-s
436 return (exp(ln_eta) + eta_0 - exp(eta_0_slash)) / 10.0 / 1e6;
437 }
438}
439
440CoolPropDbl TransportRoutines::viscosity_methanol_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
441 CoolPropDbl B_eta = NAN, C_eta = NAN, epsilon_over_k = 577.87, /* [K]*/
442 sigma0 = 0.3408e-9, /* [m] */
443 delta = 0.4575, /* NOT the reduced density, that is rhor here*/
444 N_A = 6.02214129e23, M = 32.04216, /* kg/kmol */
445 T = HEOS.T();
446 CoolPropDbl rhomolar = HEOS.rhomolar();
447
448 CoolPropDbl B_eta_star = NAN, C_eta_star = NAN;
449 CoolPropDbl Tstar = T / epsilon_over_k; // [no units]
450 CoolPropDbl rhor = HEOS.rhomass() / 273;
451 CoolPropDbl Tr = T / 512.6;
452
453 // Rainwater-Friend initial density terms
454 { // Scoped here so that we can re-use the b variable
455 CoolPropDbl b[9] = {-19.572881, 219.73999, -1015.3226, 2471.01251, -3375.1717, 2491.6597, -787.26086, 14.085455, -0.34664158};
456 CoolPropDbl t[9] = {0, -0.25, -0.5, -0.75, -1.0, -1.25, -1.5, -2.5, -5.5};
457 CoolPropDbl summer = 0;
458 for (unsigned int i = 0; i < 9; ++i) {
459 summer += b[i] * pow(Tstar, t[i]);
460 }
461 B_eta_star = summer; // [no units]
462 B_eta = N_A * pow(sigma0, 3) * B_eta_star; // [m^3/mol]
463
464 CoolPropDbl c[2] = {1.86222085e-3, 9.990338};
465 C_eta_star = c[0] * pow(Tstar, 3) * exp(c[1] * pow(Tstar, static_cast<CoolPropDbl>(-0.5))); // [no units]
466 C_eta = pow(N_A * pow(sigma0, 3), 2) * C_eta_star; // [m^6/mol^2]
467 }
468
469 CoolPropDbl eta_g = 1 + B_eta * rhomolar + C_eta * rhomolar * rhomolar;
470 CoolPropDbl a[13] = {1.16145, -0.14874, 0.52487, -0.77320, 2.16178, -2.43787, 0.95976e-3,
471 0.10225, -0.97346, 0.10657, -0.34528, -0.44557, -2.58055};
472 CoolPropDbl d[7] = {-1.181909, 0.5031030, -0.6268461, 0.5169312, -0.2351349, 5.3980235e-2, -4.9069617e-3};
473 CoolPropDbl e[10] = {0, 4.018368, -4.239180, 2.245110, -0.5750698, 2.3021026e-2, 2.5696775e-2, -6.8372749e-3, 7.2707189e-4, -2.9255711e-5};
474
475 CoolPropDbl OMEGA_22_star_LJ = a[0] * pow(Tstar, a[1]) + a[2] * exp(a[3] * Tstar) + a[4] * exp(a[5] * Tstar);
476 CoolPropDbl OMEGA_22_star_delta = a[7] * pow(Tstar, a[8]) + a[9] * exp(a[10] * Tstar) + a[11] * exp(a[12] * Tstar);
477 CoolPropDbl OMEGA_22_star_SM = OMEGA_22_star_LJ * (1 + pow(delta, 2) / (1 + a[6] * pow(delta, 6)) * OMEGA_22_star_delta);
478 CoolPropDbl eta_0 = 2.66957e-26 * sqrt(M * T) / (pow(sigma0, 2) * OMEGA_22_star_SM);
479
480 CoolPropDbl summerd = 0;
481 for (int i = 0; i < 7; ++i) {
482 summerd += d[i] / pow(Tr, i);
483 }
484 for (int j = 1; j < 10; ++j) {
485 summerd += e[j] * pow(rhor, j);
486 }
487 CoolPropDbl sigmac = 0.7193422e-9; // [m]
488 CoolPropDbl sigma_HS = summerd * sigmac; // [m]
489 CoolPropDbl b = 2 * M_PI * N_A * pow(sigma_HS, 3) / 3; // [m^3/mol]
490 CoolPropDbl zeta = b * rhomolar / 4; // [-]
491 CoolPropDbl g_sigma_HS = (1 - 0.5 * zeta) / pow(1 - zeta, 3); // [-]
492 CoolPropDbl eta_E = 1 / g_sigma_HS + 0.8 * b * rhomolar + 0.761 * g_sigma_HS * pow(b * rhomolar, 2); // [-]
493
494 CoolPropDbl f = 1 / (1 + exp(5 * (rhor - 1)));
495 return eta_0 * (f * eta_g + (1 - f) * eta_E);
496}
497
498CoolPropDbl TransportRoutines::viscosity_R23_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
499 double C1 = 1.3163, //
500 C2 = 0.1832, DeltaGstar = 771.23, rhoL = 32.174, rhocbar = 7.5114, Tc = 299.2793, DELTAeta_max = 3.967, Ru = 8.31451, molar_mass = 70.014;
501
502 double a[] = {0.4425728, -0.5138403, 0.1547566, -0.02821844, 0.001578286};
503 double e_k = 243.91, sigma = 0.4278;
504 double Tstar = HEOS.T() / e_k;
505 double logTstar = log(Tstar);
506 double Omega = exp(a[0] + a[1] * logTstar + a[2] * pow(logTstar, 2) + a[3] * pow(logTstar, 3) + a[4] * pow(logTstar, 4));
507 double eta_DG = 1.25 * 0.021357 * sqrt(molar_mass * HEOS.T()) / (sigma * sigma * Omega); // uPa-s
508
509 double rhobar = HEOS.rhomolar() / 1000; // [mol/L]
510 double eta_L = C2 * (rhoL * rhoL) / (rhoL - rhobar) * sqrt(HEOS.T()) * exp(rhobar / (rhoL - rhobar) * DeltaGstar / (Ru * HEOS.T()));
511
512 double chi = rhobar - rhocbar;
513 double tau = HEOS.T() - Tc;
514
515 double DELTAeta_c = 4 * DELTAeta_max / ((exp(chi) + exp(-chi)) * (exp(tau) + exp(-tau)));
516
517 return (pow((rhoL - rhobar) / rhoL, C1) * eta_DG + pow(rhobar / rhoL, C1) * eta_L + DELTAeta_c) / 1e6;
518}
519
520CoolPropDbl TransportRoutines::viscosity_o_xylene_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
521 // From CAO, JPCRD, 2016
522 double D[] = {-2.05581e-3, 2.38762, 0, 10.4497, 15.9587};
523 double n[] = {10.3, 3.3, 25, 0.7, 0.4};
524 double E[] = {2.65651e-3, 0, 1.77616e-12, -18.2446, 0};
525 double k[] = {0.8, 0, 4.4};
526 double Tr = HEOS.T() / 630.259, rhor = HEOS.rhomolar() / 1000.0 / 2.6845;
527
528 double A0 = -1.4933, B0 = 473.2, C0 = -57033, T = HEOS.T();
529 double ln_Seta = A0 + B0 / T + C0 / (T * T);
530 double eta0 = 0.22225 * sqrt(T) / exp(ln_Seta); // [uPa-s]
531
532 double A1 = 13.2814, B1 = -10862.4, C1 = 1664060, rho_molL = HEOS.rhomolar() / 1000.0;
533 double eta1 = (A1 + B1 / T + C1 / (T * T)) * rho_molL; // [uPa-s]
534
535 double f = (D[0] + E[0] * pow(Tr, -k[0])) * pow(rhor, n[0]) + D[1] * pow(rhor, n[1]) + E[2] * pow(rhor, n[2]) / pow(Tr, k[2])
536 + (D[3] * rhor + E[3] * Tr) * pow(rhor, n[3]) + D[4] * pow(rhor, n[4]);
537 double DELTAeta = pow(rhor, 2.0 / 3.0) * sqrt(Tr) * f; // [uPa-s]
538
539 return (eta0 + eta1 + DELTAeta) / 1e6;
540}
541CoolPropDbl TransportRoutines::viscosity_m_xylene_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
542 // From CAO, JPCRD, 2016
543 double D[] = {-0.268950, -0.0290018, 0, 14.7728, 17.1128};
544 double n[] = {6.8, 3.3, 22.0, 0.6, 0.4};
545 double E[] = {0.320971, 0, 1.72866e-10, -18.9852, 0};
546 double k[] = {0.3, 0, 3.2};
547 double Tr = HEOS.T() / 616.89, rhor = HEOS.rhomolar() / 1000.0 / 2.665;
548
549 double A0 = -1.4933, B0 = 473.2, C0 = -57033, T = HEOS.T();
550 double ln_Seta = A0 + B0 / T + C0 / (T * T);
551 double eta0 = 0.22115 * sqrt(T) / exp(ln_Seta); // [uPa-s]
552
553 double A1 = 13.2814, B1 = -10862.4, C1 = 1664060, rho_molL = HEOS.rhomolar() / 1000.0;
554 double eta1 = (A1 + B1 / T + C1 / (T * T)) * rho_molL; // [uPa-s]
555
556 double f = (D[0] + E[0] * pow(Tr, -k[0])) * pow(rhor, n[0]) + D[1] * pow(rhor, n[1]) + E[2] * pow(rhor, n[2]) / pow(Tr, k[2])
557 + (D[3] * rhor + E[3] * Tr) * pow(rhor, n[3]) + D[4] * pow(rhor, n[4]);
558 double DELTAeta = pow(rhor, 2.0 / 3.0) * sqrt(Tr) * f; // [uPa-s]
559
560 return (eta0 + eta1 + DELTAeta) / 1e6; // [Pa-s]
561}
562CoolPropDbl TransportRoutines::viscosity_p_xylene_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
563 // From Balogun, JPCRD, 2016
564 double Tr = HEOS.T() / 616.168, rhor = HEOS.rhomolar() / 1000.0 / 2.69392;
565
566 double A0 = -1.4933, B0 = 473.2, C0 = -57033, T = HEOS.T();
567 double ln_Seta = A0 + B0 / T + C0 / (T * T);
568 double eta0 = 0.22005 * sqrt(T) / exp(ln_Seta); // [uPa-s]
569
570 double A1 = 13.2814, B1 = -10862.4, C1 = 1664060, rho_molL = HEOS.rhomolar() / 1000.0;
571 double eta1 = (A1 + B1 / T + C1 / (T * T)) * rho_molL; // [uPa-s]
572
573 double sum1 = 122.919 * pow(rhor, 1.5) - 282.329 * pow(rhor, 2) + 279.348 * pow(rhor, 3) - 146.776 * pow(rhor, 4) + 28.361 * pow(rhor, 5)
574 - 0.004585 * pow(rhor, 11);
575 double sum2 = 15.337 * pow(rhor, 1.5) - 0.0004382 * pow(rhor, 11) + 0.00002307 * pow(rhor, 15);
576 double DELTAeta = pow(rhor, 2.0 / 3.0) * (sum1 + 1 / sqrt(Tr) * sum2);
577
578 return (eta0 + eta1 + DELTAeta) / 1e6; // [Pa-s]
579}
580
581CoolPropDbl TransportRoutines::viscosity_dilute_ethane(HelmholtzEOSMixtureBackend& HEOS) {
582 double C[] = {
583 0, -3.0328138281, 16.918880086, -37.189364917, 41.288861858, -24.615921140, 8.9488430959, -1.8739245042, 0.20966101390, -9.6570437074e-3};
584 double OMEGA_2_2 = 0, e_k = 245, Tstar = NAN;
585
586 Tstar = HEOS.T() / e_k;
587 for (int i = 1; i <= 9; i++) {
588 OMEGA_2_2 += C[i] * pow(Tstar, (i - 1) / 3.0 - 1);
589 }
590
591 return 12.0085 * sqrt(Tstar) * OMEGA_2_2 / 1e6; //[Pa-s]
592}
593CoolPropDbl TransportRoutines::viscosity_dilute_cyclohexane(HelmholtzEOSMixtureBackend& HEOS) {
594 // From Tariq, JPCRD, 2014
595 CoolPropDbl T = HEOS.T();
596 CoolPropDbl S_eta = exp(-1.5093 + 364.87 / T - 39537 / pow(T, 2)); //[nm^2]
597 return 0.19592 * sqrt(T) / S_eta / 1e6; //[Pa-s]
598}
599
600CoolPropDbl TransportRoutines::viscosity_dilute_CO2_LaeseckeJPCRD2017(HelmholtzEOSMixtureBackend& HEOS) {
601 // From Laesecke, JPRCD, 2016
602 double eta0 = NAN, den = NAN;
603 double T = HEOS.T();
604
605 double a[] = {1749.354893188350, -369.069300007128, 5423856.34887691, -2.21283852168356, -269503.247933569, 73145.021531826, 5.34368649509278};
606
607 // Eq. (4) from Laesecke, JPRCD, 2016
608 den = a[0] + a[1] * pow(T, 1.0 / 6.0) + a[2] * exp(a[3] * pow(T, 1.0 / 3.0)) + (a[4] + a[5] * pow(T, 1.0 / 3.0)) / exp(pow(T, 1.0 / 3.0))
609 + a[6] * sqrt(T);
610 eta0 = 0.0010055 * sqrt(T) / den; // [Pa-s]
611 return eta0;
612}
613
614CoolPropDbl TransportRoutines::viscosity_ethane_higher_order_hardcoded(HelmholtzEOSMixtureBackend& HEOS) {
615 double r[] = {0, 1, 1, 2, 2, 2, 3, 3, 4, 4, 1, 1};
616 double s[] = {0, 0, 1, 0, 1, 1.5, 0, 2, 0, 1, 0, 1};
617 double g[] = {0, 0.47177003, -0.23950311, 0.39808301, -0.27343335, 0.35192260,
618 -0.21101308, -0.00478579, 0.07378129, -0.030435255, -0.30435286, 0.001215675};
619
620 double sum1 = 0, sum2 = 0, tau = 305.33 / HEOS.T(), delta = HEOS.rhomolar() / 6870;
621
622 for (int i = 1; i <= 9; ++i) {
623 sum1 += g[i] * pow(delta, r[i]) * pow(tau, s[i]);
624 }
625 for (int i = 10; i <= 11; ++i) {
626 sum2 += g[i] * pow(delta, r[i]) * pow(tau, s[i]);
627 }
628 return 15.977 * sum1 / (1 + sum2) / 1e6;
629}
630CoolPropDbl TransportRoutines::viscosity_Chung(HelmholtzEOSMixtureBackend& HEOS) {
631 // Retrieve values from the state class
632 CoolProp::ViscosityChungData& data = HEOS.components[0].transport.viscosity_Chung;
633
634 double a0[] = {0, 6.32402, 0.12102e-2, 5.28346, 6.62263, 19.74540, -1.89992, 24.27450, 0.79716, -0.23816, 0.68629e-1};
635 double a1[] = {0, 50.41190, -0.11536e-2, 254.20900, 38.09570, 7.63034, -12.53670, 3.44945, 1.11764, 0.67695e-1, 0.34793};
636 double a2[] = {0, -51.68010, -0.62571e-2, -168.48100, -8.46414, -14.35440, 4.98529, -11.29130, 0.12348e-1, -0.81630, 0.59256};
637 double a3[] = {0, 1189.02000, 0.37283e-1, 3898.27000, 31.41780, 31.52670, -18.15070, 69.34660, -4.11661, 4.02528, -0.72663};
638 double A[11];
639
640 if (HEOS.is_pure_or_pseudopure) {
641 double Vc_cm3mol = 1 / (data.rhomolar_critical / 1e6); // [cm^3/mol]
642 double acentric = data.acentric; // [-]
643 double M_gmol = data.molar_mass * 1000.0; // [g/mol]
644 double Tc = data.T_critical; // [K]
645 double mu_D = data.dipole_moment_D; // [D]
646 double kappa = 0;
647
648 double mu_r = 131.3 * mu_D / sqrt(Vc_cm3mol * Tc); // [-]
649
650 for (int i = 1; i <= 10; ++i) {
651 A[i] = a0[i] + a1[i] * acentric + a2[i] * pow(mu_r, 4) + a3[i] * kappa;
652 }
653 double F_c = 1 - 0.2756 * acentric + 0.059035 * pow(mu_r, 4) + kappa; // [-]
654 double epsilon_over_k = Tc / 1.2593; // [K]
655
656 double rho_molcm3 = HEOS.rhomolar() / 1e6;
657 double T = HEOS.T();
658 double Tstar = T / epsilon_over_k;
659 double Omega_2_2 = 1.16145 * pow(Tstar, -0.14874) + 0.52487 * exp(-0.77320 * Tstar) + 2.16178 * exp(-2.43787 * Tstar)
660 - 6.435e-4 * pow(Tstar, 0.14874) * sin(18.0323 * pow(Tstar, -0.76830) - 7.27371); // [-]
661 double eta0_P = 4.0785e-5 * sqrt(M_gmol * T) / (pow(Vc_cm3mol, 2.0 / 3.0) * Omega_2_2) * F_c; // [P]
662
663 double Y = rho_molcm3 * Vc_cm3mol / 6.0;
664 double G_1 = (1.0 - 0.5 * Y) / pow(1 - Y, 3);
665 double G_2 = (A[1] * (1 - exp(-A[4] * Y)) / Y + A[2] * G_1 * exp(A[5] * Y) + A[3] * G_1) / (A[1] * A[4] + A[2] + A[3]);
666 double eta_k_P = eta0_P * (1 / G_2 + A[6] * Y); // [P]
667
668 double eta_p_P = (36.344e-6 * sqrt(M_gmol * Tc) / pow(Vc_cm3mol, 2.0 / 3.0)) * A[7] * pow(Y, 2) * G_2
669 * exp(A[8] + A[9] / Tstar + A[10] / pow(Tstar, 2)); // [P]
670
671 return (eta_k_P + eta_p_P) / 10.0; // [P] -> [Pa*s]
672 } else {
673 throw NotImplementedError("TransportRoutines::viscosity_Chung is only for pure and pseudo-pure");
674 }
675}
676
677CoolPropDbl TransportRoutines::conductivity_dilute_ratio_polynomials(HelmholtzEOSMixtureBackend& HEOS) {
678 if (HEOS.is_pure_or_pseudopure) {
679 // Retrieve values from the state class
680 CoolProp::ConductivityDiluteRatioPolynomialsData& data = HEOS.components[0].transport.conductivity_dilute.ratio_polynomials;
681
682 CoolPropDbl summer1 = 0, summer2 = 0, Tr = HEOS.T() / data.T_reducing;
683 for (std::size_t i = 0; i < data.A.size(); ++i) {
684 summer1 += data.A[i] * pow(Tr, data.n[i]);
685 }
686 for (std::size_t i = 0; i < data.B.size(); ++i) {
687 summer2 += data.B[i] * pow(Tr, data.m[i]);
688 }
689
690 return summer1 / summer2;
691 } else {
692 throw NotImplementedError("TransportRoutines::conductivity_dilute_ratio_polynomials is only for pure and pseudo-pure");
693 }
694};
695
696CoolPropDbl TransportRoutines::conductivity_residual_polynomial(HelmholtzEOSMixtureBackend& HEOS) {
697 if (HEOS.is_pure_or_pseudopure) {
698 // Retrieve values from the state class
699 CoolProp::ConductivityResidualPolynomialData& data = HEOS.components[0].transport.conductivity_residual.polynomials;
700
701 CoolPropDbl summer = 0, tau = data.T_reducing / HEOS.T(), delta = HEOS.keyed_output(CoolProp::iDmass) / data.rhomass_reducing;
702 for (std::size_t i = 0; i < data.B.size(); ++i) {
703 summer += data.B[i] * pow(tau, data.t[i]) * pow(delta, data.d[i]);
704 }
705 return summer;
706 } else {
707 throw NotImplementedError("TransportRoutines::conductivity_residual_polynomial is only for pure and pseudo-pure");
708 }
709};
710
711CoolPropDbl TransportRoutines::conductivity_residual_polynomial_and_exponential(HelmholtzEOSMixtureBackend& HEOS) {
712 if (HEOS.is_pure_or_pseudopure) {
713 // Retrieve values from the state class
714 CoolProp::ConductivityResidualPolynomialAndExponentialData& data =
715 HEOS.components[0].transport.conductivity_residual.polynomial_and_exponential;
716
717 CoolPropDbl summer = 0, tau = HEOS.tau(), delta = HEOS.delta();
718 for (std::size_t i = 0; i < data.A.size(); ++i) {
719 summer += data.A[i] * pow(tau, data.t[i]) * pow(delta, data.d[i]) * exp(-data.gamma[i] * pow(delta, data.l[i]));
720 }
721 return summer;
722 } else {
723 throw NotImplementedError("TransportRoutines::conductivity_residual_polynomial_and_exponential is only for pure and pseudo-pure");
724 }
725};
726
727CoolPropDbl TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers(HelmholtzEOSMixtureBackend& HEOS) {
728 if (HEOS.is_pure_or_pseudopure) {
729 // Olchowy and Sengers cross-over term
730
731 // Retrieve values from the state class
732 CoolProp::ConductivityCriticalSimplifiedOlchowySengersData& data = HEOS.components[0].transport.conductivity_critical.Olchowy_Sengers;
733
734 double k = data.k, R0 = data.R0, nu = data.nu, gamma = data.gamma, GAMMA = data.GAMMA, zeta0 = data.zeta0, qD = data.qD,
735 Tc = HEOS.get_reducing_state().T, // [K]
736 rhoc = HEOS.get_reducing_state().rhomolar, // [mol/m^3]
737 Pcrit = HEOS.get_reducing_state().p, // [Pa]
738 Tref = NAN, // [K]
739 cp = NAN, cv = NAN, delta = NAN, num = NAN, zeta = NAN, mu = NAN, pi = M_PI, OMEGA_tilde = NAN, OMEGA_tilde0 = NAN;
740
741 if (ValidNumber(data.T_ref))
742 Tref = data.T_ref;
743 else
744 Tref = 1.5 * Tc;
745
746 delta = HEOS.delta();
747
748 double dp_drho = HEOS.gas_constant() * HEOS.T() * (1 + 2 * delta * HEOS.dalphar_dDelta() + delta * delta * HEOS.d2alphar_dDelta2());
749 double X = Pcrit / pow(rhoc, 2) * HEOS.rhomolar() / dp_drho;
750
751 double tau_ref = Tc / Tref;
752 double dp_drho_ref = HEOS.gas_constant() * Tref
753 * (1 + 2 * delta * HEOS.calc_alphar_deriv_nocache(0, 1, HEOS.mole_fractions, tau_ref, delta)
754 + delta * delta * HEOS.calc_alphar_deriv_nocache(0, 2, HEOS.mole_fractions, tau_ref, delta));
755 double Xref = Pcrit / pow(rhoc, 2) * HEOS.rhomolar() / dp_drho_ref * Tref / HEOS.T();
756 num = X - Xref;
757
758 // No critical enhancement if numerator is negative, zero, or just a tiny bit positive due to roundoff
759 // See also Lemmon, IJT, 2004, page 27
760 if (num < DBL_EPSILON * 10)
761 return 0.0;
762 else
763 zeta = zeta0 * pow(num / GAMMA, nu / gamma); //[m]
764
765 cp = HEOS.cpmolar(); //[J/mol/K]
766 cv = HEOS.cvmolar(); //[J/mol/K]
767 mu = HEOS.viscosity(); //[Pa-s]
768
769 OMEGA_tilde = 2.0 / pi * ((cp - cv) / cp * atan(zeta * qD) + cv / cp * (zeta * qD)); //[-]
770 OMEGA_tilde0 = 2.0 / pi * (1.0 - exp(-1.0 / (1.0 / (qD * zeta) + 1.0 / 3.0 * (zeta * qD) * (zeta * qD) / delta / delta))); //[-]
771
772 double lambda = HEOS.rhomolar() * cp * R0 * k * HEOS.T() / (6 * pi * mu * zeta) * (OMEGA_tilde - OMEGA_tilde0); //[W/m/K]
773 return lambda; //[W/m/K]
774 } else {
775 throw NotImplementedError("TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers is only for pure and pseudo-pure");
776 }
777};
778
779CoolPropDbl TransportRoutines::conductivity_critical_hardcoded_R123(HelmholtzEOSMixtureBackend& HEOS) {
780 double a13 = 0.486742e-2, a14 = -100, a15 = -7.08535;
781 return a13 * exp(a14 * pow(HEOS.tau() - 1, 4) + a15 * pow(HEOS.delta() - 1, 2));
782};
783
784CoolPropDbl TransportRoutines::conductivity_critical_hardcoded_CO2_ScalabrinJPCRD2006(HelmholtzEOSMixtureBackend& HEOS) {
785 CoolPropDbl nc = 0.775547504e-3 * 4.81384, Tr = HEOS.T() / 304.1282, alpha = NAN, rhor = HEOS.keyed_output(iDmass) / 467.6;
786 static CoolPropDbl a[] = {0.0, 3.0, 6.70697, 0.94604, 0.30, 0.30, 0.39751, 0.33791, 0.77963, 0.79857, 0.90, 0.02, 0.20};
787
788 // Equation 6 from Scalabrin
789 alpha = 1 - a[10] * acosh(1 + a[11] * pow(pow(1 - Tr, 2), a[12]));
790
791 // Equation 5 from Scalabrin
792 CoolPropDbl numer = rhor * exp(-pow(rhor, a[1]) / a[1] - pow(a[2] * (Tr - 1), 2) - pow(a[3] * (rhor - 1), 2));
793 CoolPropDbl braced = (1 - 1 / Tr) + a[4] * pow(pow(rhor - 1, 2), 0.5 / a[5]);
794 CoolPropDbl denom = pow(pow(pow(braced, 2), a[6]) + pow(pow(a[7] * (rhor - alpha), 2), a[8]), a[9]);
795 return nc * numer / denom;
796}
797
798CoolPropDbl TransportRoutines::conductivity_dilute_hardcoded_CO2(HelmholtzEOSMixtureBackend& HEOS) {
799
800 double e_k = 251.196, Tstar = NAN;
801 double b[] = {0.4226159, 0.6280115, -0.5387661, 0.6735941, 0, 0, -0.4362677, 0.2255388};
802 double c[] = {0, 2.387869e-2, 4.350794, -10.33404, 7.981590, -1.940558};
803
804 //Vesovic Eq. 31 [no units]
805 double summer = 0;
806 for (int i = 1; i <= 5; i++)
807 summer += c[i] * pow(HEOS.T() / 100.0, 2 - i);
808 double cint_k = 1.0 + exp(-183.5 / HEOS.T()) * summer;
809
810 //Vesovic Eq. 12 [no units]
811 double r = sqrt(2.0 / 5.0 * cint_k);
812
813 // According to REFPROP, 1+r^2 = cp-2.5R. This is unclear to me but seems to suggest that cint/k is the difference
814 // between the ideal gas specific heat and a monatomic specific heat of 5/2*R. Using the form of cint/k from Vesovic
815 // does not yield exactly the correct values
816
817 Tstar = HEOS.T() / e_k;
818 //Vesovic Eq. 30 [no units]
819 summer = 0;
820 for (int i = 0; i <= 7; i++)
821 summer += b[i] / pow(Tstar, i);
822 double Gstar_lambda = summer;
823
824 //Vesovic Eq. 29 [W/m/K]
825 double lambda_0 = 0.475598e-3 * sqrt(HEOS.T()) * (1 + r * r) / Gstar_lambda;
826
827 return lambda_0;
828}
829
830CoolPropDbl TransportRoutines::conductivity_dilute_hardcoded_CO2_HuberJPCRD2016(HelmholtzEOSMixtureBackend& HEOS) {
831
832 double tau = HEOS.tau();
833 double l[] = {0.0151874307, 0.0280674040, 0.0228564190, -0.00741624210};
834 // Huber 2016 Eq. (3)
835 double lambda_0 = pow(tau, -0.5) / (l[0] + l[1] * tau + l[2] * pow(tau, 2) + l[3] * pow(tau, 3)); // [mW/m/K]
836
837 return lambda_0 / 1000;
838}
839
840CoolPropDbl TransportRoutines::conductivity_dilute_hardcoded_ethane(HelmholtzEOSMixtureBackend& HEOS) {
841
842 double e_k = 245.0;
843 double tau = 305.33 / HEOS.T(), Tstar = HEOS.T() / e_k;
844 double fint = 1.7104147 - 0.6936482 / Tstar;
845 double lambda_0 = 0.276505e-3 * (HEOS.calc_viscosity_dilute() * 1e6) * (3.75 - fint * (tau * tau * HEOS.d2alpha0_dTau2() + 1.5)); //[W/m/K]
846
847 return lambda_0;
848}
849
850CoolPropDbl TransportRoutines::conductivity_dilute_eta0_and_poly(HelmholtzEOSMixtureBackend& HEOS) {
851
852 if (HEOS.is_pure_or_pseudopure) {
853 CoolProp::ConductivityDiluteEta0AndPolyData& E = HEOS.components[0].transport.conductivity_dilute.eta0_and_poly;
854
855 double eta0_uPas = HEOS.calc_viscosity_dilute() * 1e6; // [uPa-s]
856 double summer = E.A[0] * eta0_uPas;
857 for (std::size_t i = 1; i < E.A.size(); ++i)
858 summer += E.A[i] * pow(static_cast<CoolPropDbl>(HEOS.tau()), E.t[i]);
859 return summer;
860 } else {
861 throw NotImplementedError("TransportRoutines::conductivity_dilute_eta0_and_poly is only for pure and pseudo-pure");
862 }
863}
864
865CoolPropDbl TransportRoutines::conductivity_hardcoded_heavywater(HelmholtzEOSMixtureBackend& HEOS) {
866 double Tbar = HEOS.T() / 643.847, rhobar = HEOS.rhomass() / 358;
867 double A[] = {1.00000, 37.3223, 22.5485, 13.0465, 0, -2.60735};
868 double lambda0 = A[0] + A[1] * Tbar + A[2] * POW2(Tbar) + A[3] * POW3(Tbar) + A[4] * POW4(Tbar) + A[5] * POW5(Tbar);
869 double Be = -2.506, B[] = {-167.310, 483.656, -191.039, 73.0358, -7.57467};
870 double DELTAlambda = B[0] * (1 - exp(Be * rhobar)) + B[1] * rhobar + B[2] * POW2(rhobar) + B[3] * POW3(rhobar) + B[4] * POW4(rhobar);
871 double f_1 = exp(0.144847 * Tbar + -5.64493 * POW2(Tbar));
872 double f_2 = exp(-2.80000 * POW2(rhobar - 1)) - 0.080738543 * exp(-17.9430 * POW2(rhobar - 0.125698));
873 double tau = Tbar / (std::abs(Tbar - 1.1) + 1.1);
874 double f_3 = 1 + exp(60 * (tau - 1) + 20);
875 double f_4 = 1 + exp(100 * (tau - 1) + 15);
876 double DELTAlambda_c = 35429.6 * f_1 * f_2 * (1 + POW2(f_2) * (5000.0e6 * POW4(f_1) / f_3 + 3.5 * f_2 / f_4));
877 double DELTAlambda_L = -741.112 * pow(f_1, 1.2) * (1 - exp(-pow(rhobar / 2.5, 10)));
878 double lambdabar = lambda0 + DELTAlambda + DELTAlambda_c + DELTAlambda_L;
879 return lambdabar * 0.742128e-3;
880}
881
882CoolPropDbl TransportRoutines::conductivity_hardcoded_water(HelmholtzEOSMixtureBackend& HEOS) {
883
884 double L[5][6] = {{1.60397357, -0.646013523, 0.111443906, 0.102997357, -0.0504123634, 0.00609859258},
885 {2.33771842, -2.78843778, 1.53616167, -0.463045512, 0.0832827019, -0.00719201245},
886 {2.19650529, -4.54580785, 3.55777244, -1.40944978, 0.275418278, -0.0205938816},
887 {-1.21051378, 1.60812989, -0.621178141, 0.0716373224, 0, 0},
888 {-2.7203370, 4.57586331, -3.18369245, 1.1168348, -0.19268305, 0.012913842}};
889
890 double lambdabar_0 = NAN, lambdabar_1 = NAN, lambdabar_2 = NAN, rhobar = NAN, Tbar = NAN, sum = NAN;
891 double Tstar = 647.096, rhostar = 322, pstar = 22064000, lambdastar = 1e-3, mustar = 1e-6;
892 double xi = NAN;
893 int i = 0, j = 0;
894 double R = 461.51805; //[J/kg/K]
895
896 Tbar = HEOS.T() / Tstar;
897 rhobar = HEOS.keyed_output(CoolProp::iDmass) / rhostar;
898
899 // Dilute gas contribution
900 lambdabar_0 =
901 sqrt(Tbar) / (2.443221e-3 + 1.323095e-2 / Tbar + 6.770357e-3 / pow(Tbar, 2) - 3.454586e-3 / pow(Tbar, 3) + 4.096266e-4 / pow(Tbar, 4));
902
903 sum = 0;
904 for (i = 0; i <= 4; i++) {
905 for (j = 0; j <= 5; j++) {
906 sum += L[i][j] * powInt(1.0 / Tbar - 1.0, i) * powInt(rhobar - 1, j);
907 }
908 }
909 // Finite density contribution
910 lambdabar_1 = exp(rhobar * sum);
911
912 double nu = 0.630, GAMMA = 177.8514, gamma = 1.239, xi_0 = 0.13, Lambda_0 = 0.06, Tr_bar = 1.5, qd_bar = 1 / 0.4, pi = 3.141592654,
913 delta = HEOS.delta();
914
915 double drhodp = 1 / (R * HEOS.T() * (1 + 2 * rhobar * HEOS.dalphar_dDelta() + rhobar * rhobar * HEOS.d2alphar_dDelta2()));
916 double drhobar_dpbar = pstar / rhostar * drhodp;
917 double drhodp_Trbar = 1
918 / (R * Tr_bar * Tstar
919 * (1 + 2 * rhobar * HEOS.calc_alphar_deriv_nocache(0, 1, HEOS.mole_fractions, 1 / Tr_bar, delta)
920 + delta * delta * HEOS.calc_alphar_deriv_nocache(0, 2, HEOS.mole_fractions, 1 / Tr_bar, delta)));
921 double drhobar_dpbar_Trbar = pstar / rhostar * drhodp_Trbar;
922 double cp = HEOS.cpmass(); // [J/kg/K]
923 double cv = HEOS.cvmass(); // [J/kg/K]
924 double cpbar = cp / R; //[-]
925 double mubar = HEOS.viscosity() / mustar;
926 double DELTAchibar_T = rhobar * (drhobar_dpbar - drhobar_dpbar_Trbar * Tr_bar / Tbar);
927 if (DELTAchibar_T < 0)
928 xi = 0;
929 else
930 xi = xi_0 * pow(DELTAchibar_T / Lambda_0, nu / gamma);
931 double y = qd_bar * xi;
932
933 double Z = NAN;
934 double kappa = cp / cv;
935 if (y < 1.2e-7)
936 Z = 0;
937 else
938 Z = 2 / (pi * y) * (((1 - 1 / kappa) * atan(y) + y / kappa) - (1 - exp(-1 / (1 / y + y * y / 3 / rhobar / rhobar))));
939
940 lambdabar_2 = GAMMA * rhobar * cpbar * Tbar / mubar * Z;
941
942 return (lambdabar_0 * lambdabar_1 + lambdabar_2) * lambdastar;
943}
944
945CoolPropDbl TransportRoutines::conductivity_hardcoded_R23(HelmholtzEOSMixtureBackend& HEOS) {
946
947 double B1 = -2.5370, // [mW/m/K]
948 B2 = 0.05366, // [mW/m/K^2]
949 C1 = 0.94215, // [-]
950 C2 = 0.14914, // [mW/m/K^2]
951 DeltaGstar = 2508.58, //[J/mol]
952 rhoL = 68.345, // [mol/dm^3] = [mol/L]
953 rhocbar = 7.5114, // [mol/dm^3]
954 DELTAlambda_max = 25, //[mW/m/K]
955 Ru = 8.31451, // [J/mol/K]
956 Tc = 299.2793, //[K]
957 T = HEOS.T(); //[K]
958
959 double lambda_DG = B1 + B2 * T;
960
961 double rhobar = HEOS.rhomolar() / 1000; // [mol/L]
962 double lambda_L = C2 * (rhoL * rhoL) / (rhoL - rhobar) * sqrt(T) * exp(rhobar / (rhoL - rhobar) * DeltaGstar / (Ru * T));
963
964 double chi = rhobar - rhocbar;
965 double tau = T - Tc;
966
967 double DELTAlambda_c = 4 * DELTAlambda_max / ((exp(chi) + exp(-chi)) * (exp(tau) + exp(-tau)));
968
969 return (pow((rhoL - rhobar) / rhoL, C1) * lambda_DG + pow(rhobar / rhoL, C1) * lambda_L + DELTAlambda_c) / 1e3;
970}
971
972CoolPropDbl TransportRoutines::conductivity_critical_hardcoded_ammonia(HelmholtzEOSMixtureBackend& HEOS) {
973
974 /*
975 From "Thermal Conductivity of Ammonia in a Large
976 Temperature and Pressure Range Including the Critical Region"
977 by R. Tufeu, D.Y. Ivanov, Y. Garrabos, B. Le Neindre,
978 Bereicht der Bunsengesellschaft Phys. Chem. 88 (1984) 422-427
979 */
980
981 double T = HEOS.T(), Tc = 405.4, rhoc = 235, rho = NAN;
982 double LAMBDA = 1.2, nu = 0.63, gamma = 1.24, DELTA = 0.50, t = NAN, zeta_0_plus = 1.34e-10, a_zeta = 1, GAMMA_0_plus = 0.423e-8;
983 double pi = 3.141592654, a_chi = NAN, k_B = 1.3806504e-23, X_T = NAN, DELTA_lambda = NAN, dPdT = NAN, eta_B = NAN, DELTA_lambda_id = NAN,
984 DELTA_lambda_i = NAN;
985
986 rho = HEOS.keyed_output(CoolProp::iDmass);
987 t = std::abs((T - Tc) / Tc);
988 a_chi = a_zeta / 0.7;
989 eta_B = (2.60 + 1.6 * t) * 1e-5;
990 dPdT = (2.18 - 0.12 / exp(17.8 * t)) * 1e5; // [Pa-K]
991 X_T = 0.61 * rhoc + 16.5 * log(t);
992 // Along the critical isochore (only a function of temperature) (Eq. 9)
993 DELTA_lambda_i = LAMBDA * (k_B * T * T) / (6 * pi * eta_B * (zeta_0_plus * pow(t, -nu) * (1 + a_zeta * pow(t, DELTA)))) * dPdT * dPdT
994 * GAMMA_0_plus * pow(t, -gamma) * (1 + a_chi * pow(t, DELTA));
995 DELTA_lambda_id = DELTA_lambda_i * exp(-36 * t * t);
996 if (rho < 0.6 * rhoc) {
997 DELTA_lambda = DELTA_lambda_id * (X_T * X_T) / (X_T * X_T + powInt(0.6 * rhoc - 0.96 * rhoc, 2)) * powInt(rho, 2) / powInt(0.6 * rhoc, 2);
998 } else {
999 DELTA_lambda = DELTA_lambda_id * (X_T * X_T) / (X_T * X_T + powInt(rho - 0.96 * rhoc, 2));
1000 }
1001
1002 return DELTA_lambda;
1003}
1004
1005CoolPropDbl TransportRoutines::conductivity_hardcoded_helium(HelmholtzEOSMixtureBackend& HEOS) {
1006 /*
1007 What an incredibly annoying formulation! Implied coefficients?? Not cool.
1008 */
1009 double rhoc = 68.0, lambda_e = NAN, lambda_c = NAN, T = HEOS.T(), rho = HEOS.rhomass();
1010 double summer = 3.739232544 / T - 2.620316969e1 / T / T + 5.982252246e1 / T / T / T - 4.926397634e1 / T / T / T / T;
1011 double lambda_0 = 2.7870034e-3 * pow(T, 7.034007057e-1) * exp(summer);
1012 double c[] = {1.862970530e-4, -7.275964435e-7, -1.427549651e-4, 3.290833592e-5, -5.213335363e-8, 4.492659933e-8,
1013 -5.924416513e-9, 7.087321137e-6, -6.013335678e-6, 8.067145814e-7, 3.995125013e-7};
1014 // Equation 17
1015 lambda_e = (c[0] + c[1] * T + c[2] * pow(T, 1 / 3.0) + c[3] * pow(T, 2.0 / 3.0)) * rho
1016 + (c[4] + c[5] * pow(T, 1.0 / 3.0) + c[6] * pow(T, 2.0 / 3.0)) * rho * rho * rho
1017 + (c[7] + c[8] * pow(T, 1.0 / 3.0) + c[9] * pow(T, 2.0 / 3.0) + c[10] / T) * rho * rho * log(rho / rhoc);
1018
1019 // Critical component
1020 lambda_c = 0.0;
1021
1022 if (3.5 < T && T < 12) {
1023 double x0 = 0.392, E1 = 2.8461, E2 = 0.27156, beta = 0.3554, gamma = 1.1743, delta = 4.304, rhoc_crit = 69.158, Tc = 5.18992, pc = 2.2746e5;
1024
1025 double DeltaT = std::abs(1 - T / Tc), DeltaRho = std::abs(1 - rho / rhoc_crit);
1026 double eta = HEOS.viscosity(); // [Pa-s]
1027 double K_T = HEOS.isothermal_compressibility(), K_Tprime = NAN, K_Tbar = NAN;
1028 double dpdT = HEOS.first_partial_deriv(CoolProp::iP, CoolProp::iT, CoolProp::iDmolar);
1029
1030 double W = pow(DeltaT / 0.2, 2) + pow(DeltaRho / 0.25, 2);
1031
1032 if (W > 1) {
1033 K_Tbar = K_T;
1034 } else {
1035 double x = pow(DeltaT / DeltaRho, 1 / beta);
1036 double h = E1 * (1 + x / x0) * pow(1 + E2 * pow(1 + x / x0, 2 / beta), (gamma - 1) / (2 * beta));
1037
1045 double dhdx =
1046 E1
1047 * (E2 * pow((x + x0) / x0, 2 / beta) * (gamma - 1) * pow(E2 * pow((x + x0) / x0, 2 / beta) + 1, (1.0 / 2.0) * (gamma - 1) / beta)
1048 + pow(beta, 2) * pow(E2 * pow((x + x0) / x0, 2 / beta) + 1, (1.0 / 2.0) * (2 * beta + gamma - 1) / beta))
1049 / (pow(beta, 2) * x0 * (E2 * pow((x + x0) / x0, 2 / beta) + 1));
1050 // Right-hand-side of Equation 9
1051 double RHS = pow(DeltaRho, delta - 1) * (delta * h - x / beta * dhdx);
1052 K_Tprime = 1 / (RHS * pow(rho / rhoc_crit, 2) * pc);
1053 K_Tbar = W * K_T + (1 - W) * K_Tprime;
1054 }
1055
1056 // 3.4685233d-17 and 3.726229668d0 are "magical" coefficients that are present in the REFPROP source to yield the right values. Not clear why these values are needed.
1057 // Also, the form of the critical term in REFPROP does not agree with Hands paper. EL and MH from NIST are not sure where these coefficients come from.
1058 lambda_c =
1059 3.4685233e-17 * 3.726229668 * sqrt(K_Tbar) * pow(T, 2) / rho / eta * pow(dpdT, 2) * exp(-18.66 * pow(DeltaT, 2) - 4.25 * pow(DeltaRho, 4));
1060 }
1061 return lambda_0 + lambda_e + lambda_c;
1062}
1063CoolPropDbl TransportRoutines::conductivity_hardcoded_methane(HelmholtzEOSMixtureBackend& HEOS) {
1064
1065 double delta = HEOS.rhomolar() / 10139.0, tau = 190.55 / HEOS.T();
1066 double lambda_dilute = NAN, lambda_residual = NAN, lambda_critical = NAN;
1067
1068 // Viscosity formulation from Friend, JPCRD, 1989
1069 // Dilute
1070 double C[] = {
1071 0, -3.0328138281, 16.918880086, -37.189364917, 41.288861858, -24.615921140, 8.9488430959, -1.8739245042, 0.20966101390, -9.6570437074e-3};
1072 double OMEGA22_summer = 0;
1073 double t = HEOS.T() / 174.0;
1074 for (int i = 1; i <= 9; ++i) {
1075 OMEGA22_summer += C[i] * pow(t, (i - 1.0) / 3.0 - 1.0);
1076 }
1077 double eta_dilute = 10.50 * sqrt(t) * OMEGA22_summer;
1078 double re[] = {0, 1, 1, 2, 2, 2, 3, 3, 4, 4, 1, 1};
1079 double se[] = {0, 0, 1, 0, 1, 1.5, 0, 2, 0, 1, 0, 1};
1080 double ge[] = {0, 0.41250137, -0.14390912, 0.10366993, 0.40287464, -0.24903524,
1081 -0.12953131, 0.06575776, 0.02566628, -0.03716526, -0.38798341, 0.03533815};
1082 double summer1 = 0, summer2 = 0;
1083 for (int i = 1; i <= 9; ++i) {
1084 summer1 += ge[i] * pow(delta, re[i]) * pow(tau, se[i]);
1085 }
1086 for (int i = 10; i <= 11; ++i) {
1087 summer2 += ge[i] * pow(delta, re[i]) * pow(tau, se[i]);
1088 }
1089 double eta_residual = 12.149 * summer1 / (1 + summer2);
1090 double eta = eta_residual + eta_dilute;
1091
1092 // Dilute
1093 double f_int = 1.458850 - 0.4377162 / t;
1094 lambda_dilute = 0.51828 * eta_dilute * (3.75 - f_int * (POW2(HEOS.tau()) * HEOS.d2alpha0_dTau2() + 1.5)); // [mW/m/K]
1095 // Residual
1096 double rl[] = {0, 1, 3, 4, 4, 5, 5, 2};
1097 double sl[] = {0, 0, 0, 0, 1, 0, 1, 0};
1098 double jl[] = {0, 2.4149207, 0.55166331, -0.52837734, 0.073809553, 0.24465507, -0.047613626, 1.5554612};
1099 double summer = 0;
1100 for (int i = 1; i <= 6; ++i) {
1101 summer += jl[i] * pow(delta, rl[i]) * pow(tau, sl[i]);
1102 }
1103 double delta_sigma_star = 1.0; // Looks like a typo in Friend - should be 1 instead of 11
1104 if (HEOS.T() < HEOS.T_critical() && HEOS.rhomolar() < HEOS.rhomolar_critical()) {
1105 delta_sigma_star = HEOS.saturation_ancillary(iDmolar, 1, iT, HEOS.T()) / HEOS.keyed_output(CoolProp::irhomolar_critical);
1106 }
1107 lambda_residual = 6.29638 * (summer + jl[7] * POW2(delta) / delta_sigma_star); // [mW/m/K]
1108 // Critical region
1109 double Tstar = 1 - 1 / tau;
1110 double rhostar = 1 - delta;
1111 double F_T = 2.646, F_rho = 2.678, F_A = -0.637;
1112 double F = exp(-F_T * sqrt(std::abs(Tstar)) - F_rho * POW2(rhostar) - F_A * rhostar);
1113 double CHI_T_star = NAN;
1114 if (std::abs(Tstar) < 0.03) {
1115 if (std::abs(rhostar) < 1e-16) {
1116 // Equation 26
1117 const double LAMBDA = 0.0801, gamma = 1.190;
1118 CHI_T_star = LAMBDA * pow(std::abs(Tstar), -gamma);
1119 } else if (std::abs(rhostar) < 0.03) {
1120 // Equation 23
1121 const double beta = 0.355, W = -1.401, S = -6.098, E = 0.287, a = 3.352, b = 0.732, R = 0.535, Q = 0.1133;
1122 double OMEGA = W * Tstar * pow(std::abs(rhostar), -1 / beta);
1123 double theta = 1;
1124 if (Tstar < -pow(std::abs(rhostar), -1 / beta) / S) {
1125 theta = 1 + E * pow(1 + S * Tstar * pow(std::abs(rhostar), -1 / beta), 2 * beta);
1126 }
1127 CHI_T_star = Q * pow(std::abs(rhostar), -a) * pow(theta, b) / (theta + OMEGA * (theta + R));
1128 } else {
1129 // Equation 19a
1130 CHI_T_star = 0.28631 * delta * tau / (1 + 2 * delta * HEOS.dalphar_dDelta() + POW2(delta) * HEOS.d2alphar_dDelta2());
1131 }
1132 } else {
1133 // Equation 19a
1134 CHI_T_star = 0.28631 * delta * tau / (1 + 2 * delta * HEOS.dalphar_dDelta() + POW2(delta) * HEOS.d2alphar_dDelta2());
1135 }
1136
1137 lambda_critical = 91.855 / (eta * POW2(tau)) * POW2(1 + delta * HEOS.dalphar_dDelta() - delta * tau * HEOS.d2alphar_dDelta_dTau())
1138 * pow(CHI_T_star, 0.4681) * F; //[mW/m/K]
1139 return (lambda_dilute + lambda_residual + lambda_critical) * 0.001;
1140}
1141
1142void TransportRoutines::conformal_state_solver(HelmholtzEOSMixtureBackend& HEOS, HelmholtzEOSMixtureBackend& HEOS_Reference, CoolPropDbl& T0,
1143 CoolPropDbl& rhomolar0) {
1144 int iter = 0;
1145 double resid = 9e30, resid_old = 9e30;
1146 CoolPropDbl alphar = HEOS.alphar();
1147 CoolPropDbl Z = HEOS.keyed_output(iZ);
1148
1149 Eigen::Vector2d r;
1150 Eigen::Matrix2d J;
1151 HEOS_Reference.specify_phase(iphase_gas); // Something homogeneous, not checked
1152 // Update the reference fluid with the conformal state
1153 HEOS_Reference.update_DmolarT_direct(rhomolar0, T0);
1154 do {
1155 CoolPropDbl dtau_dT = -HEOS_Reference.T_critical() / (T0 * T0);
1156 CoolPropDbl ddelta_drho = 1 / HEOS_Reference.rhomolar_critical();
1157 // Independent variables are T0 and rhomolar0, residuals are matching alphar and Z
1158 r(0) = HEOS_Reference.alphar() - alphar;
1159 r(1) = HEOS_Reference.keyed_output(iZ) - Z;
1160 J(0, 0) = HEOS_Reference.dalphar_dTau() * dtau_dT;
1161 J(0, 1) = HEOS_Reference.dalphar_dDelta() * ddelta_drho;
1162 // Z = 1+delta*dalphar_ddelta(tau,delta)
1163 // dZ_dT
1164 J(1, 0) = HEOS_Reference.delta() * HEOS_Reference.d2alphar_dDelta_dTau() * dtau_dT;
1165 // dZ_drho
1166 J(1, 1) = (HEOS_Reference.delta() * HEOS_Reference.d2alphar_dDelta2() + HEOS_Reference.dalphar_dDelta()) * ddelta_drho;
1167 // Step in v obtained from Jv = -r
1168 Eigen::Vector2d v = J.colPivHouseholderQr().solve(-r);
1169 bool good_solution = false;
1170 double T0_init = HEOS_Reference.T(), rhomolar0_init = HEOS_Reference.rhomolar();
1171 // Calculate the old residual after the last step
1172 resid_old = sqrt(POW2(r(0)) + POW2(r(1)));
1173 // Geometric Newton-step halving (~10 iters from 1.0 to ~1/1024).
1174 for (double frac = 1.0; frac > 0.001; frac /= 2) { // NOLINT(cert-flp30-c)
1175 try {
1176 // Calculate new values
1177 double T_new = T0_init + frac * v(0);
1178 double rhomolar_new = rhomolar0_init + frac * v(1);
1179 // Update state with step
1180 HEOS_Reference.update_DmolarT_direct(rhomolar_new, T_new);
1181 resid = sqrt(POW2(HEOS_Reference.alphar() - alphar) + POW2(HEOS_Reference.keyed_output(iZ) - Z));
1182 if (resid > resid_old) {
1183 continue;
1184 }
1185 good_solution = true;
1186 T0 = T_new;
1187 rhomolar0 = rhomolar_new;
1188 break;
1189 } catch (...) {
1190 continue;
1191 }
1192 }
1193 if (!good_solution) {
1194 throw ValueError(format("Not able to get a solution"));
1195 }
1196 iter++;
1197 if (iter > 50) {
1198 throw ValueError(format("conformal_state_solver took too many iterations; residual is %g; prior was %g", resid, resid_old));
1199 }
1200 } while (std::abs(resid) > 1e-9);
1201}
1202
1203CoolPropDbl TransportRoutines::viscosity_ECS(HelmholtzEOSMixtureBackend& HEOS, HelmholtzEOSMixtureBackend& HEOS_Reference) {
1204 // Collect some parameters
1205 CoolPropDbl M = HEOS.molar_mass(), M0 = HEOS_Reference.molar_mass(), Tc = HEOS.T_critical(), Tc0 = HEOS_Reference.T_critical(),
1206 rhocmolar = HEOS.rhomolar_critical(), rhocmolar0 = HEOS_Reference.rhomolar_critical();
1207
1208 // Get a reference to the ECS data
1209 CoolProp::ViscosityECSVariables& ECS = HEOS.components[0].transport.viscosity_ecs;
1210
1211 // The correction polynomial psi_eta
1212 double psi = 0;
1213 for (std::size_t i = 0; i < ECS.psi_a.size(); i++) {
1214 psi += ECS.psi_a[i] * pow(HEOS.rhomolar() / ECS.psi_rhomolar_reducing, ECS.psi_t[i]);
1215 }
1216
1217 // The dilute gas portion for the fluid of interest [Pa-s]
1218 CoolPropDbl eta_dilute = viscosity_dilute_kinetic_theory(HEOS);
1219
1220 // ************************************
1221 // Start with a guess for theta and phi
1222 // ************************************
1223 CoolPropDbl theta = 1;
1224 CoolPropDbl phi = 1;
1225
1226 // The equivalent substance reducing ratios
1227 CoolPropDbl f = Tc / Tc0 * theta;
1228 CoolPropDbl h = rhocmolar0 / rhocmolar * phi; // Must be the ratio of MOLAR densities!!
1229
1230 // To be solved for
1231 CoolPropDbl T0 = HEOS.T() / f;
1232 CoolPropDbl rhomolar0 = HEOS.rhomolar() * h;
1233
1234 // **************************
1235 // Solver for conformal state
1236 // **************************
1237
1238 //
1239 HEOS_Reference.specify_phase(iphase_gas); // something homogeneous
1240
1241 conformal_state_solver(HEOS, HEOS_Reference, T0, rhomolar0);
1242
1243 // Update the reference fluid with the updated conformal state
1244 HEOS_Reference.update_DmolarT_direct(rhomolar0 * psi, T0);
1245
1246 // Recalculate ESRR
1247 f = HEOS.T() / T0;
1248 h = rhomolar0 / HEOS.rhomolar(); // Must be the ratio of MOLAR densities!!
1249
1250 // **********************
1251 // Remaining calculations
1252 // **********************
1253
1254 // The reference fluid's contribution to the viscosity [Pa-s]
1255 CoolPropDbl eta_resid = HEOS_Reference.calc_viscosity_background();
1256
1257 // The F factor
1258 CoolPropDbl F_eta = sqrt(f) * pow(h, -static_cast<CoolPropDbl>(2.0L / 3.0L)) * sqrt(M / M0);
1259
1260 // The total viscosity considering the contributions of the fluid of interest and the reference fluid [Pa-s]
1261 CoolPropDbl eta = eta_dilute + eta_resid * F_eta;
1262
1263 return eta;
1264}
1265
1266CoolPropDbl TransportRoutines::viscosity_rhosr(HelmholtzEOSMixtureBackend& HEOS) {
1267
1268 // Get a reference to the data
1269 const CoolProp::ViscosityRhoSrVariables& data = HEOS.components[0].transport.viscosity_rhosr;
1270
1271 // The dilute gas portion for the fluid of interest [Pa-s]
1272 CoolPropDbl eta_dilute = viscosity_dilute_kinetic_theory(HEOS);
1273
1274 // Calculate x
1275 double x = HEOS.rhomolar() * HEOS.gas_constant() * (HEOS.tau() * HEOS.dalphar_dTau() - HEOS.alphar()) / data.rhosr_critical;
1276
1277 // Crossover variable
1278 double psi_liq = 1 / (1 + exp(-100.0 * (x - 2)));
1279
1280 // Evaluated using Horner's method
1281 const std::vector<double>&cL = data.c_liq, cV = data.c_vap;
1282 double f_liq = cL[0] + x * (cL[1] + x * (cL[2] + x * (cL[3])));
1283 double f_vap = cV[0] + x * (cV[1] + x * (cV[2] + x * (cV[3])));
1284
1285 // Evaluate the reference fluid
1286 double etastar_ref = exp(psi_liq * f_liq + (1 - psi_liq) * f_vap);
1287
1288 // Get the non-dimensionalized viscosity
1289 double etastar_fluid = 1 + data.C * (etastar_ref - 1);
1290
1291 return etastar_fluid * eta_dilute;
1292}
1293
1294CoolPropDbl TransportRoutines::conductivity_ECS(HelmholtzEOSMixtureBackend& HEOS, HelmholtzEOSMixtureBackend& HEOS_Reference) {
1295 // Collect some parameters
1296 CoolPropDbl M = HEOS.molar_mass(), M_kmol = M * 1000, M0 = HEOS_Reference.molar_mass(), Tc = HEOS.T_critical(), Tc0 = HEOS_Reference.T_critical(),
1297 rhocmolar = HEOS.rhomolar_critical(), rhocmolar0 = HEOS_Reference.rhomolar_critical(), R_u = HEOS.gas_constant(),
1298 R = HEOS.gas_constant() / HEOS.molar_mass(), //[J/kg/K]
1299 R_kJkgK = R_u / M_kmol;
1300
1301 // Get a reference to the ECS data
1302 CoolProp::ConductivityECSVariables& ECS = HEOS.components[0].transport.conductivity_ecs;
1303
1304 // The correction polynomial psi_eta in rho/rho_red
1305 double psi = 0;
1306 for (std::size_t i = 0; i < ECS.psi_a.size(); ++i) {
1307 psi += ECS.psi_a[i] * pow(HEOS.rhomolar() / ECS.psi_rhomolar_reducing, ECS.psi_t[i]);
1308 }
1309
1310 // The correction polynomial f_int in T/T_red
1311 double fint = 0;
1312 for (std::size_t i = 0; i < ECS.f_int_a.size(); ++i) {
1313 fint += ECS.f_int_a[i] * pow(HEOS.T() / ECS.f_int_T_reducing, ECS.f_int_t[i]);
1314 }
1315
1316 // The dilute gas density for the fluid of interest [uPa-s]
1317 CoolPropDbl eta_dilute = viscosity_dilute_kinetic_theory(HEOS) * 1e6;
1318
1319 // The mass specific ideal gas constant-pressure specific heat [J/kg/K]
1320 CoolPropDbl cp0 = HEOS.calc_cpmolar_idealgas() / HEOS.molar_mass();
1321
1322 // The internal contribution to the thermal conductivity [W/m/K]
1323 CoolPropDbl lambda_int = fint * eta_dilute * (cp0 - 2.5 * R) / 1e3;
1324
1325 // The dilute gas contribution to the thermal conductivity [W/m/K]
1326 CoolPropDbl lambda_dilute = 15.0e-3 / 4.0 * R_kJkgK * eta_dilute;
1327
1328 // ************************************
1329 // Start with a guess for theta and phi
1330 // ************************************
1331
1332 CoolPropDbl theta = 1;
1333 CoolPropDbl phi = 1;
1334
1335 // The equivalent substance reducing ratios
1336 CoolPropDbl f = Tc / Tc0 * theta;
1337 CoolPropDbl h = rhocmolar0 / rhocmolar * phi; // Must be the ratio of MOLAR densities!!
1338
1339 // Initial values for the conformal state
1340 CoolPropDbl T0 = HEOS.T() / f;
1341 CoolPropDbl rhomolar0 = HEOS.rhomolar() * h;
1342
1343 // **************************
1344 // Solver for conformal state
1345 // **************************
1346
1347 try {
1348 conformal_state_solver(HEOS, HEOS_Reference, T0, rhomolar0);
1349 } catch (std::exception& e) {
1350 throw ValueError(format("Conformal state solver failed; error was %s", e.what()));
1351 }
1352
1353 // Update the reference fluid with the conformal state
1354 HEOS_Reference.update(DmolarT_INPUTS, rhomolar0 * psi, T0);
1355
1356 // Recalculate ESRR
1357 f = HEOS.T() / T0;
1358 h = rhomolar0 / HEOS.rhomolar(); // Must be the ratio of MOLAR densities!!
1359
1360 // The reference fluid's contribution to the conductivity [W/m/K]
1361 CoolPropDbl lambda_resid = HEOS_Reference.calc_conductivity_background();
1362
1363 // The F factor
1364 CoolPropDbl F_lambda = sqrt(f) * pow(h, static_cast<CoolPropDbl>(-2.0 / 3.0)) * sqrt(M0 / M);
1365
1366 // The critical contribution from the fluid of interest [W/m/K]
1367 CoolPropDbl lambda_critical = conductivity_critical_simplified_Olchowy_Sengers(HEOS);
1368
1369 // The total thermal conductivity considering the contributions of the fluid of interest and the reference fluid [Pa-s]
1370 CoolPropDbl lambda = lambda_int + lambda_dilute + lambda_resid * F_lambda + lambda_critical;
1371
1372 return lambda;
1373}
1374
1375}; /* namespace CoolProp */
Generated by doxygen 1.9.4