17 1.16145 * pow(Tstar,
static_cast<CoolPropDbl>(-0.14874)) + 0.52487 * exp(-0.77320 * Tstar) + 2.16178 * exp(-2.43787 * Tstar);
20 return 26.692e-9 * sqrt(molar_mass_kgkmol * HEOS.
T()) / (pow(sigma_nm, 2) * OMEGA22);
22 throw NotImplementedError(
"TransportRoutines::viscosity_dilute_kinetic_theory is only for pure and pseudo-pure");
30 const std::vector<CoolPropDbl>&a = data.
a, &t = data.
t;
37 const CoolPropDbl molar_mass_kgkmol = molar_mass * 1000;
43 for (std::size_t i = 0; i < a.size(); ++i) {
44 summer += a[i] * pow(lnTstar, t[i]);
49 return C * sqrt(molar_mass_kgkmol * HEOS.
T()) / (pow(sigma_nm, 2) * S);
51 throw NotImplementedError(
"TransportRoutines::viscosity_dilute_collision_integral is only for pure and pseudo-pure");
59 const std::vector<CoolPropDbl>&a = data.
a, &t = data.
t;
62 for (std::size_t i = 0; i < a.size(); ++i) {
63 summer += a[i] * pow(T, t[i]);
67 throw NotImplementedError(
"TransportRoutines::viscosity_dilute_powers_of_T is only for pure and pseudo-pure");
74 const std::vector<CoolPropDbl>&a = data.
a, &t = data.
t;
76 for (std::size_t i = 0; i < a.size(); ++i) {
77 summer += a[i] * pow(Tr, t[i]);
81 throw NotImplementedError(
"TransportRoutines::viscosity_dilute_powers_of_Tr is only for pure and pseudo-pure");
89 HEOS.
components[0].transport.viscosity_dilute.collision_integral_powers_of_Tstar;
90 const std::vector<CoolPropDbl>&a = data.
a, &t = data.
t;
93 for (std::size_t i = 0; i < a.size(); ++i) {
94 summer += a[i] * pow(Tstar, t[i]);
96 return data.
C * sqrt(HEOS.
T()) / summer;
98 throw NotImplementedError(
"TransportRoutines::viscosity_dilute_collision_integral_powers_of_T is only for pure and pseudo-pure");
104 HEOS.
components[0].transport.viscosity_higher_order.modified_Batschinski_Hildebrand;
110 for (
unsigned int i = 0; i < HO.
a.size(); ++i) {
111 S += HO.
a[i] * pow(delta, HO.
d1[i]) * pow(tau, HO.
t1[i]) * exp(HO.
gamma[i] * pow(delta, HO.
l[i]));
116 for (
unsigned int i = 0; i < HO.
f.size(); ++i) {
117 F += HO.
f[i] * pow(delta, HO.
d2[i]) * pow(tau, HO.
t2[i]);
122 for (
unsigned int i = 0; i < HO.
g.size(); ++i) {
123 summer_numer += HO.
g[i] * pow(tau, HO.
h[i]);
126 for (
unsigned int i = 0; i < HO.
p.size(); ++i) {
127 summer_denom += HO.
p[i] * pow(tau, HO.
q[i]);
132 return S + F * (1 / (delta0 - delta) - 1 / delta0);
134 throw NotImplementedError(
"TransportRoutines::viscosity_higher_order_modified_Batschinski_Hildebrand is only for pure and pseudo-pure");
142 const std::vector<CoolPropDbl>&b = data.
b, &t = data.
t;
149 for (
unsigned int i = 0; i < b.size(); ++i) {
150 summer += b[i] * pow(Tstar, t[i]);
153 B_eta = 6.02214129e23 * pow(sigma, 3) * B_eta_star;
156 throw NotImplementedError(
"TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend is only for pure and pseudo-pure");
165 const std::vector<CoolPropDbl>&n = data.
n, &d = data.
d, &t = data.
t;
171 for (
unsigned int i = 0; i < n.size(); ++i) {
172 summer += n[i] * pow(delta, d[i]) * pow(tau, t[i]);
176 throw NotImplementedError(
"TransportRoutines::viscosity_initial_density_dependence_empirical is only for pure and pseudo-pure");
180 static void visc_Helper(
double Tbar,
double rhobar,
double* mubar_0,
double* mubar_1) {
181 std::vector<std::vector<CoolPropDbl>> H(6, std::vector<CoolPropDbl>(7, 0));
186 *mubar_0 = 100.0 * sqrt(Tbar) / (1.67752 + 2.20462 / Tbar + 0.6366564 /
powInt(Tbar, 2) - 0.241605 /
powInt(Tbar, 3));
189 for (i = 0; i <= 5; i++) {
190 for (j = 0; j <= 6; j++) {
196 H[0][0] = 5.20094e-1;
197 H[1][0] = 8.50895e-2;
199 H[3][0] = -2.89555e-1;
201 H[0][1] = 2.22531e-1;
202 H[1][1] = 9.99115e-1;
205 H[5][1] = 1.20573e-1;
207 H[0][2] = -2.81378e-1;
208 H[1][2] = -9.06851e-1;
209 H[2][2] = -7.72479e-1;
210 H[3][2] = -4.89837e-1;
211 H[4][2] = -2.57040e-1;
213 H[0][3] = 1.61913e-1;
214 H[1][3] = 2.57399e-1;
216 H[0][4] = -3.25372e-2;
217 H[3][4] = 6.98452e-2;
219 H[4][5] = 8.72102e-3;
221 H[3][6] = -4.35673e-3;
222 H[5][6] = -5.93264e-4;
226 for (i = 0; i <= 5; i++) {
227 for (j = 0; j <= 6; j++) {
228 sum +=
powInt(1 / Tbar - 1, i) * (H[i][j] *
powInt(rhobar - 1, j));
231 *mubar_1 = exp(rhobar * sum);
234 double Tbar = HEOS.
T() / 643.847, rhobar = HEOS.
rhomass() / 358;
235 double A[] = {1.000000, 0.940695, 0.578377, -0.202044};
236 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};
237 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};
238 double Bij[] = {0.4864192, -0.2448372, -0.8702035, 0.8716056, -1.051126, 0.3458395, 0.3509007, 1.315436, 1.297752,
239 1.353448, -0.2847572, -1.037026, -1.287846, -0.02148229, 0.07013759, 0.4660127, 0.2292075, -0.4857462,
240 0.01641220, -0.02884911, 0.1607171, -0.009603846, -0.01163815, -0.008239587, 0.004559914, -0.003886659};
241 double mu0 = sqrt(Tbar) / (A[0] + A[1] / Tbar + A[2] /
POW2(Tbar) + A[3] /
POW3(Tbar));
243 for (
int i = 0; i < 26; ++i) {
244 summer += Bij[i] * pow(1 / Tbar - 1, I[i]) * pow(rhobar - 1, J[i]);
246 double mu1 = exp(rhobar * summer);
247 double mubar = mu0 * mu1;
248 return 55.2651e-6 * mubar;
251 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, Tstar, rhostar;
252 double delta, tau, mubar_0, mubar_1, mubar_2, drhodp, drhodp_R, DeltaChibar, zeta, w, L, Y, psi_D, Tbar, rhobar;
253 double drhobar_dpbar, drhobar_dpbar_R, R_Water;
258 Tbar = HEOS.
T() / Tstar;
259 rhobar = HEOS.
rhomass() / rhostar;
263 visc_Helper(Tbar, rhobar, &mubar_0, &mubar_1);
271 drhobar_dpbar = pstar / rhostar * drhodp;
275 / (R_Water * Tbar_R * Tstar
278 drhobar_dpbar_R = pstar / rhostar * drhodp_R;
280 DeltaChibar = rhobar * (drhobar_dpbar - drhobar_dpbar_R * Tbar_R / Tbar);
281 if (DeltaChibar < 0) DeltaChibar = 0;
282 zeta = zeta_0 * pow(DeltaChibar / LAMBDA_0, nu / gamma);
283 if (zeta < 0.3817016416) {
284 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));
286 psi_D = acos(pow(1 +
powInt(qd * zeta, 2), -1.0 / 2.0));
287 w = sqrt(std::abs((qc * zeta - 1) / (qc * zeta + 1))) * tan(psi_D / 2.0);
289 L = log((1 + w) / (1 - w));
291 L = 2 * atan(std::abs(w));
293 Y = 1.0 / 12.0 * sin(3 * psi_D) - 1 / (4 * qc * zeta) * sin(2 * psi_D)
294 + 1.0 /
powInt(qc * zeta, 2) * (1 - 5.0 / 4.0 *
powInt(qc * zeta, 2)) * sin(psi_D)
295 - 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);
297 mubar_2 = exp(x_mu * Y);
299 return (mubar_0 * mubar_1 * mubar_2) / 1e6;
303 CoolPropDbl c[] = {19.919216, -2.6557905, -135.904211, -7.9962719, -11.014795, -10.113817};
304 return 1e-6 * pow(
static_cast<double>(rhor), 2.0 / 3.0) * sqrt(Tr)
305 * ((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 CoolPropDbl c[] = {0, 6.43449673e-6, 4.56334068e-2, 2.32797868e-1, 9.58326120e-1, 1.27941189e-1, 3.63576595e-1};
311 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));
315 CoolPropDbl c[] = {-9.98945, 86.06260, 2.74872, 1.11130, -1.0, -134.1330, -352.473, 6.60989, 88.4174};
316 return 1e-6 * pow(rhor,
static_cast<CoolPropDbl>(2.0 / 3.0)) * sqrt(Tr)
317 * (c[0] * pow(rhor, 2) + c[1] * rhor / (c[2] + c[3] * Tr + c[4] * rhor)
318 + (c[5] * rhor + c[6] * pow(rhor, 2)) / (c[7] + c[8] * pow(rhor, 2)));
325 double c[] = {2.53402335 / 1e6, -9.724061002 / 1e6, 0.469437316, 158.5571631, 72.42916856 / 1e6,
326 10.60751253, 8.628373915, -6.61346441, -2.212724566};
327 return pow(rhor,
static_cast<CoolPropDbl>(2.0 / 3.0)) * sqrt(Tr)
328 * (c[0] / Tr + c[1] / (c[2] + Tr + c[3] * rhor * rhor)
329 + c[4] * (1 + rhor) / (c[5] + c[6] * Tr + c[7] * rhor + rhor * rhor + c[8] * rhor * Tr));
337 double c[] = {0, 22.15000 / 1e6, -15.00870 / 1e6, 3.71791 / 1e6, 77.72818 / 1e6, 9.73449, 9.51900, -6.34076, -2.51909};
338 return pow(rhor,
static_cast<CoolPropDbl>(2.0L / 3.0L)) * sqrt(Tr)
339 * (c[1] * rhor + c[2] * pow(rhor, 2) + c[3] * pow(rhor, 3)
340 + c[4] * rhor / (c[5] + c[6] * Tr + c[7] * rhor + rhor * rhor + c[8] * rhor * Tr));
344 double c1 = 0.360603235428487, c2 = 0.121550806591497, gamma = 8.06282737481277;
345 double Tt = HEOS.
Ttriple(), rho_tL = 1178.53;
346 double Tr = HEOS.
T() / Tt, rhor = HEOS.
rhomass() / rho_tL;
348 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);
350 double residual = eta_tL * (c1 * Tr * pow(rhor, 3) + (pow(rhor, 2) + pow(rhor, gamma)) / (Tr - c2));
360 double psi1 = exp(tau) - F.
c1;
361 double psi2 = exp(pow(tau, 2)) - F.
c2;
363 double ki = (F.
Ai[0] + F.
Ai[1] * psi1 + F.
Ai[2] * psi2) * tau;
365 double ka = (F.
Aa[0] + F.
Aa[1] * psi1 + F.
Aa[2] * psi2) * pow(tau, F.
Na);
366 double kr = (F.
Ar[0] + F.
Ar[1] * psi1 + F.
Ar[2] * psi2) * pow(tau, F.
Nr);
367 double kaa = (F.
Aaa[0] + F.
Aaa[1] * psi1 + F.
Aaa[2] * psi2) * pow(tau, F.
Naa);
372 krr = (F.
Arr[0] + F.
Arr[1] * psi1 + F.
Arr[2] * psi2) * pow(tau, F.
Nrr);
375 if (!F.
Aii.empty()) {
376 kii = (F.
Aii[0] + F.
Aii[1] * psi1 + F.
Aii[2] * psi2) * pow(tau, F.
Nii);
378 if (!F.
Arrr.empty() && !F.
Aaaa.empty()) {
379 krrr = (F.
Arrr[0] + F.
Arrr[1] * psi1 + F.
Arrr[2] * psi2) * pow(tau, F.
Nrrr);
380 kaaa = (F.
Aaaa[0] + F.
Aaaa[1] * psi1 + F.
Aaaa[2] * psi2) * pow(tau, F.
Naaa);
383 double p = HEOS.
p() / 1e5;
388 double deltapr = pr - pid;
390 double eta_f = ka * pa + kr * deltapr + ki * pid + kaa * pa * pa + kdrdr * deltapr * deltapr + krr * pr * pr + kii * pid * pid
391 + krrr * pr * pr * pr + kaaa * pa * pa * pa;
395 throw NotImplementedError(
"TransportRoutines::viscosity_higher_order_friction_theory is only for pure and pseudo-pure");
400 double eta_0, eta_0_slash, eta_E_slash, B, C, D, ln_eta, x;
418 B = -47.5295259 / x + 87.6799309 - 42.0741589 * x + 8.33128289 * x * x - 0.589252385 * x * x * x;
419 C = 547.309267 / x - 904.870586 + 431.404928 * x - 81.4504854 * x * x + 5.37008433 * x * x * x;
420 D = -1684.39324 / x + 3331.08630 - 1632.19172 * x + 308.804413 * x * x - 20.2936367 * x * x * x;
421 eta_0_slash = -0.135311743 / x + 1.00347841 + 1.20654649 * x - 0.149564551 * x * x + 0.012520841 * x * x * x;
422 eta_E_slash = rho * B + rho * rho * C + rho * rho * rho * D;
425 ln_eta = eta_0_slash + eta_E_slash;
427 return exp(ln_eta) / 10.0 / 1e6;
429 ln_eta = eta_0_slash + eta_E_slash;
430 eta_0 = 196 * pow(T,
static_cast<CoolPropDbl>(0.71938)) * exp(12.451 / T - 295.67 / T / T - 4.1249);
432 return (exp(ln_eta) + eta_0 - exp(eta_0_slash)) / 10.0 / 1e6;
440 N_A = 6.02214129e23, M = 32.04216,
451 CoolPropDbl b[9] = {-19.572881, 219.73999, -1015.3226, 2471.01251, -3375.1717, 2491.6597, -787.26086, 14.085455, -0.34664158};
452 CoolPropDbl t[9] = {0, -0.25, -0.5, -0.75, -1.0, -1.25, -1.5, -2.5, -5.5};
454 for (
unsigned int i = 0; i < 9; ++i) {
455 summer += b[i] * pow(Tstar, t[i]);
458 B_eta = N_A * pow(sigma0, 3) * B_eta_star;
461 C_eta_star = c[0] * pow(Tstar, 3) * exp(c[1] * pow(Tstar,
static_cast<CoolPropDbl>(-0.5)));
462 C_eta = pow(N_A * pow(sigma0, 3), 2) * C_eta_star;
465 CoolPropDbl eta_g = 1 + B_eta * rhomolar + C_eta * rhomolar * rhomolar;
466 CoolPropDbl a[13] = {1.16145, -0.14874, 0.52487, -0.77320, 2.16178, -2.43787, 0.95976e-3,
467 0.10225, -0.97346, 0.10657, -0.34528, -0.44557, -2.58055};
468 CoolPropDbl d[7] = {-1.181909, 0.5031030, -0.6268461, 0.5169312, -0.2351349, 5.3980235e-2, -4.9069617e-3};
469 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};
471 CoolPropDbl OMEGA_22_star_LJ = a[0] * pow(Tstar, a[1]) + a[2] * exp(a[3] * Tstar) + a[4] * exp(a[5] * Tstar);
472 CoolPropDbl OMEGA_22_star_delta = a[7] * pow(Tstar, a[8]) + a[9] * exp(a[10] * Tstar) + a[11] * exp(a[12] * Tstar);
473 CoolPropDbl OMEGA_22_star_SM = OMEGA_22_star_LJ * (1 + pow(delta, 2) / (1 + a[6] * pow(delta, 6)) * OMEGA_22_star_delta);
474 CoolPropDbl eta_0 = 2.66957e-26 * sqrt(M * T) / (pow(sigma0, 2) * OMEGA_22_star_SM);
477 for (
int i = 0; i < 7; ++i) {
478 summerd += d[i] / pow(Tr, i);
480 for (
int j = 1; j < 10; ++j) {
481 summerd += e[j] * pow(rhor, j);
487 CoolPropDbl g_sigma_HS = (1 - 0.5 * zeta) / pow(1 - zeta, 3);
488 CoolPropDbl eta_E = 1 / g_sigma_HS + 0.8 * b * rhomolar + 0.761 * g_sigma_HS * pow(b * rhomolar, 2);
491 return eta_0 * (f * eta_g + (1 - f) * eta_E);
496 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;
498 double a[] = {0.4425728, -0.5138403, 0.1547566, -0.02821844, 0.001578286};
499 double e_k = 243.91, sigma = 0.4278;
500 double Tstar = HEOS.
T() / e_k;
501 double logTstar = log(Tstar);
502 double Omega = exp(a[0] + a[1] * logTstar + a[2] * pow(logTstar, 2) + a[3] * pow(logTstar, 3) + a[4] * pow(logTstar, 4));
503 double eta_DG = 1.25 * 0.021357 * sqrt(molar_mass * HEOS.
T()) / (sigma * sigma * Omega);
505 double rhobar = HEOS.
rhomolar() / 1000;
506 double eta_L = C2 * (rhoL * rhoL) / (rhoL - rhobar) * sqrt(HEOS.
T()) * exp(rhobar / (rhoL - rhobar) * DeltaGstar / (Ru * HEOS.
T()));
508 double chi = rhobar - rhocbar;
509 double tau = HEOS.
T() - Tc;
511 double DELTAeta_c = 4 * DELTAeta_max / ((exp(chi) + exp(-chi)) * (exp(tau) + exp(-tau)));
513 return (pow((rhoL - rhobar) / rhoL, C1) * eta_DG + pow(rhobar / rhoL, C1) * eta_L + DELTAeta_c) / 1e6;
518 double D[] = {-2.05581e-3, 2.38762, 0, 10.4497, 15.9587};
519 double n[] = {10.3, 3.3, 25, 0.7, 0.4};
520 double E[] = {2.65651e-3, 0, 1.77616e-12, -18.2446, 0};
521 double k[] = {0.8, 0, 4.4};
522 double Tr = HEOS.
T() / 630.259, rhor = HEOS.
rhomolar() / 1000.0 / 2.6845;
524 double A0 = -1.4933, B0 = 473.2, C0 = -57033, T = HEOS.
T();
525 double ln_Seta = A0 + B0 / T + C0 / (T * T);
526 double eta0 = 0.22225 * sqrt(T) / exp(ln_Seta);
528 double A1 = 13.2814, B1 = -10862.4, C1 = 1664060, rho_molL = HEOS.
rhomolar() / 1000.0;
529 double eta1 = (A1 + B1 / T + C1 / (T * T)) * rho_molL;
531 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])
532 + (D[3] * rhor + E[3] * Tr) * pow(rhor, n[3]) + D[4] * pow(rhor, n[4]);
533 double DELTAeta = pow(rhor, 2.0 / 3.0) * sqrt(Tr) * f;
535 return (eta0 + eta1 + DELTAeta) / 1e6;
539 double D[] = {-0.268950, -0.0290018, 0, 14.7728, 17.1128};
540 double n[] = {6.8, 3.3, 22.0, 0.6, 0.4};
541 double E[] = {0.320971, 0, 1.72866e-10, -18.9852, 0};
542 double k[] = {0.3, 0, 3.2};
543 double Tr = HEOS.
T() / 616.89, rhor = HEOS.
rhomolar() / 1000.0 / 2.665;
545 double A0 = -1.4933, B0 = 473.2, C0 = -57033, T = HEOS.
T();
546 double ln_Seta = A0 + B0 / T + C0 / (T * T);
547 double eta0 = 0.22115 * sqrt(T) / exp(ln_Seta);
549 double A1 = 13.2814, B1 = -10862.4, C1 = 1664060, rho_molL = HEOS.
rhomolar() / 1000.0;
550 double eta1 = (A1 + B1 / T + C1 / (T * T)) * rho_molL;
552 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])
553 + (D[3] * rhor + E[3] * Tr) * pow(rhor, n[3]) + D[4] * pow(rhor, n[4]);
554 double DELTAeta = pow(rhor, 2.0 / 3.0) * sqrt(Tr) * f;
556 return (eta0 + eta1 + DELTAeta) / 1e6;
560 double Tr = HEOS.
T() / 616.168, rhor = HEOS.
rhomolar() / 1000.0 / 2.69392;
562 double A0 = -1.4933, B0 = 473.2, C0 = -57033, T = HEOS.
T();
563 double ln_Seta = A0 + B0 / T + C0 / (T * T);
564 double eta0 = 0.22005 * sqrt(T) / exp(ln_Seta);
566 double A1 = 13.2814, B1 = -10862.4, C1 = 1664060, rho_molL = HEOS.
rhomolar() / 1000.0;
567 double eta1 = (A1 + B1 / T + C1 / (T * T)) * rho_molL;
569 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)
570 - 0.004585 * pow(rhor, 11);
571 double sum2 = 15.337 * pow(rhor, 1.5) - 0.0004382 * pow(rhor, 11) + 0.00002307 * pow(rhor, 15);
572 double DELTAeta = pow(rhor, 2.0 / 3.0) * (sum1 + 1 / sqrt(Tr) * sum2);
574 return (eta0 + eta1 + DELTAeta) / 1e6;
579 0, -3.0328138281, 16.918880086, -37.189364917, 41.288861858, -24.615921140, 8.9488430959, -1.8739245042, 0.20966101390, -9.6570437074e-3};
580 double OMEGA_2_2 = 0, e_k = 245, Tstar;
582 Tstar = HEOS.
T() / e_k;
583 for (
int i = 1; i <= 9; i++) {
584 OMEGA_2_2 += C[i] * pow(Tstar, (i - 1) / 3.0 - 1);
587 return 12.0085 * sqrt(Tstar) * OMEGA_2_2 / 1e6;
592 CoolPropDbl S_eta = exp(-1.5093 + 364.87 / T - 39537 / pow(T, 2));
593 return 0.19592 * sqrt(T) / S_eta / 1e6;
598 double eta0, eta1, DELTAetar, den, Bstar;
602 1749.354893188350, -369.069300007128, 5423856.34887691, -2.21283852168356, -269503.247933569, 73145.021531826, 5.34368649509278};
605 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)) + a[6] * sqrt(T);
606 eta0 = 0.0010055 * sqrt(T) / den;
611 double r[] = {0, 1, 1, 2, 2, 2, 3, 3, 4, 4, 1, 1};
612 double s[] = {0, 0, 1, 0, 1, 1.5, 0, 2, 0, 1, 0, 1};
613 double g[] = {0, 0.47177003, -0.23950311, 0.39808301, -0.27343335, 0.35192260,
614 -0.21101308, -0.00478579, 0.07378129, -0.030435255, -0.30435286, 0.001215675};
616 double sum1 = 0, sum2 = 0, tau = 305.33 / HEOS.
T(), delta = HEOS.
rhomolar() / 6870;
618 for (
int i = 1; i <= 9; ++i) {
619 sum1 += g[i] * pow(delta, r[i]) * pow(tau, s[i]);
621 for (
int i = 10; i <= 11; ++i) {
622 sum2 += g[i] * pow(delta, r[i]) * pow(tau, s[i]);
624 return 15.977 * sum1 / (1 + sum2) / 1e6;
630 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};
631 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};
632 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};
633 double a3[] = {0, 1189.02000, 0.37283e-1, 3898.27000, 31.41780, 31.52670, -18.15070, 69.34660, -4.11661, 4.02528, -0.72663};
644 double mu_r = 131.3 * mu_D / sqrt(Vc_cm3mol * Tc);
646 for (
int i = 1; i <= 10; ++i) {
647 A[i] = a0[i] + a1[i] * acentric + a2[i] * pow(mu_r, 4) + a3[i] * kappa;
649 double F_c = 1 - 0.2756 * acentric + 0.059035 * pow(mu_r, 4) + kappa;
650 double epsilon_over_k = Tc / 1.2593;
652 double rho_molcm3 = HEOS.
rhomolar() / 1e6;
654 double Tstar = T / epsilon_over_k;
655 double Omega_2_2 = 1.16145 * pow(Tstar, -0.14874) + 0.52487 * exp(-0.77320 * Tstar) + 2.16178 * exp(-2.43787 * Tstar)
656 - 6.435e-4 * pow(Tstar, 0.14874) * sin(18.0323 * pow(Tstar, -0.76830) - 7.27371);
657 double eta0_P = 4.0785e-5 * sqrt(M_gmol * T) / (pow(Vc_cm3mol, 2.0 / 3.0) * Omega_2_2) * F_c;
659 double Y = rho_molcm3 * Vc_cm3mol / 6.0;
660 double G_1 = (1.0 - 0.5 * Y) / pow(1 - Y, 3);
661 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]);
662 double eta_k_P = eta0_P * (1 / G_2 + A[6] * Y);
664 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
665 * exp(A[8] + A[9] / Tstar + A[10] / pow(Tstar, 2));
667 return (eta_k_P + eta_p_P) / 10.0;
669 throw NotImplementedError(
"TransportRoutines::viscosity_Chung is only for pure and pseudo-pure");
679 for (std::size_t i = 0; i < data.
A.size(); ++i) {
680 summer1 += data.
A[i] * pow(Tr, data.
n[i]);
682 for (std::size_t i = 0; i < data.
B.size(); ++i) {
683 summer2 += data.
B[i] * pow(Tr, data.
m[i]);
686 return summer1 / summer2;
688 throw NotImplementedError(
"TransportRoutines::conductivity_dilute_ratio_polynomials is only for pure and pseudo-pure");
698 for (std::size_t i = 0; i < data.
B.size(); ++i) {
699 summer += data.
B[i] * pow(tau, data.
t[i]) * pow(delta, data.
d[i]);
703 throw NotImplementedError(
"TransportRoutines::conductivity_residual_polynomial is only for pure and pseudo-pure");
711 HEOS.
components[0].transport.conductivity_residual.polynomial_and_exponential;
714 for (std::size_t i = 0; i < data.
A.size(); ++i) {
715 summer += data.
A[i] * pow(tau, data.
t[i]) * pow(delta, data.
d[i]) * exp(-data.
gamma[i] * pow(delta, data.
l[i]));
719 throw NotImplementedError(
"TransportRoutines::conductivity_residual_polynomial_and_exponential is only for pure and pseudo-pure");
730 double k = data.
k, R0 = data.
R0, nu = data.
nu, gamma = data.
gamma, GAMMA = data.
GAMMA, zeta0 = data.
zeta0, qD = data.
qD,
735 cp, cv, delta, num, zeta, mu, pi =
M_PI, OMEGA_tilde, OMEGA_tilde0;
742 delta = HEOS.
delta();
745 double X = Pcrit / pow(rhoc, 2) * HEOS.
rhomolar() / dp_drho;
747 double tau_ref = Tc / Tref;
751 double Xref = Pcrit / pow(rhoc, 2) * HEOS.
rhomolar() / dp_drho_ref * Tref / HEOS.
T();
759 zeta = zeta0 * pow(num / GAMMA, nu / gamma);
765 OMEGA_tilde = 2.0 / pi * ((cp - cv) / cp * atan(zeta * qD) + cv / cp * (zeta * qD));
766 OMEGA_tilde0 = 2.0 / pi * (1.0 - exp(-1.0 / (1.0 / (qD * zeta) + 1.0 / 3.0 * (zeta * qD) * (zeta * qD) / delta / delta)));
768 double lambda = HEOS.
rhomolar() * cp * R0 * k * HEOS.
T() / (6 * pi * mu * zeta) * (OMEGA_tilde - OMEGA_tilde0);
771 throw NotImplementedError(
"TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers is only for pure and pseudo-pure");
776 double a13 = 0.486742e-2, a14 = -100, a15 = -7.08535;
777 return a13 * exp(a14 * pow(HEOS.
tau() - 1, 4) + a15 * pow(HEOS.
delta() - 1, 2));
782 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};
785 alpha = 1 - a[10] * acosh(1 + a[11] * pow(pow(1 - Tr, 2), a[12]));
788 CoolPropDbl numer = rhor * exp(-pow(rhor, a[1]) / a[1] - pow(a[2] * (Tr - 1), 2) - pow(a[3] * (rhor - 1), 2));
789 CoolPropDbl braced = (1 - 1 / Tr) + a[4] * pow(pow(rhor - 1, 2), 0.5 / a[5]);
790 CoolPropDbl denom = pow(pow(pow(braced, 2), a[6]) + pow(pow(a[7] * (rhor - alpha), 2), a[8]), a[9]);
791 return nc * numer / denom;
796 double e_k = 251.196, Tstar;
797 double b[] = {0.4226159, 0.6280115, -0.5387661, 0.6735941, 0, 0, -0.4362677, 0.2255388};
798 double c[] = {0, 2.387869e-2, 4.350794, -10.33404, 7.981590, -1.940558};
802 for (
int i = 1; i <= 5; i++)
803 summer += c[i] * pow(HEOS.
T() / 100.0, 2 - i);
804 double cint_k = 1.0 + exp(-183.5 / HEOS.
T()) * summer;
807 double r = sqrt(2.0 / 5.0 * cint_k);
813 Tstar = HEOS.
T() / e_k;
816 for (
int i = 0; i <= 7; i++)
817 summer += b[i] / pow(Tstar, i);
818 double Gstar_lambda = summer;
821 double lambda_0 = 0.475598e-3 * sqrt(HEOS.
T()) * (1 + r * r) / Gstar_lambda;
828 double tau = HEOS.
tau();
829 double l[]={0.0151874307, 0.0280674040, 0.0228564190, -0.00741624210};
831 double lambda_0 = pow(tau, -0.5)/(l[0] + l[1] * tau + l[2] * pow(tau, 2) + l[3] * pow(tau, 3));
833 return lambda_0/1000;
839 double tau = 305.33 / HEOS.
T(), Tstar = HEOS.
T() / e_k;
840 double fint = 1.7104147 - 0.6936482 / Tstar;
852 double summer = E.
A[0] * eta0_uPas;
853 for (std::size_t i = 1; i < E.
A.size(); ++i)
857 throw NotImplementedError(
"TransportRoutines::conductivity_dilute_eta0_and_poly is only for pure and pseudo-pure");
862 double Tbar = HEOS.
T() / 643.847, rhobar = HEOS.
rhomass() / 358;
863 double A[] = {1.00000, 37.3223, 22.5485, 13.0465, 0, -2.60735};
864 double lambda0 = A[0] + A[1] * Tbar + A[2] *
POW2(Tbar) + A[3] *
POW3(Tbar) + A[4] *
POW4(Tbar) + A[5] *
POW5(Tbar);
865 double Be = -2.506, B[] = {-167.310, 483.656, -191.039, 73.0358, -7.57467};
866 double DELTAlambda = B[0] * (1 - exp(Be * rhobar)) + B[1] * rhobar + B[2] *
POW2(rhobar) + B[3] *
POW3(rhobar) + B[4] *
POW4(rhobar);
867 double f_1 = exp(0.144847 * Tbar + -5.64493 *
POW2(Tbar));
868 double f_2 = exp(-2.80000 *
POW2(rhobar - 1)) - 0.080738543 * exp(-17.9430 *
POW2(rhobar - 0.125698));
869 double tau = Tbar / (std::abs(Tbar - 1.1) + 1.1);
870 double f_3 = 1 + exp(60 * (tau - 1) + 20);
871 double f_4 = 1 + exp(100 * (tau - 1) + 15);
872 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));
873 double DELTAlambda_L = -741.112 * pow(f_1, 1.2) * (1 - exp(-pow(rhobar / 2.5, 10)));
874 double lambdabar = lambda0 + DELTAlambda + DELTAlambda_c + DELTAlambda_L;
875 return lambdabar * 0.742128e-3;
880 double L[5][6] = {{1.60397357, -0.646013523, 0.111443906, 0.102997357, -0.0504123634, 0.00609859258},
881 {2.33771842, -2.78843778, 1.53616167, -0.463045512, 0.0832827019, -0.00719201245},
882 {2.19650529, -4.54580785, 3.55777244, -1.40944978, 0.275418278, -0.0205938816},
883 {-1.21051378, 1.60812989, -0.621178141, 0.0716373224, 0, 0},
884 {-2.7203370, 4.57586331, -3.18369245, 1.1168348, -0.19268305, 0.012913842}};
886 double lambdabar_0, lambdabar_1, lambdabar_2, rhobar, Tbar, sum;
887 double Tstar = 647.096, rhostar = 322, pstar = 22064000, lambdastar = 1e-3, mustar = 1e-6;
890 double R = 461.51805;
892 Tbar = HEOS.
T() / Tstar;
897 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));
900 for (i = 0; i <= 4; i++) {
901 for (j = 0; j <= 5; j++) {
902 sum += L[i][j] *
powInt(1.0 / Tbar - 1.0, i) *
powInt(rhobar - 1, j);
906 lambdabar_1 = exp(rhobar * sum);
908 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,
909 delta = HEOS.
delta();
912 double drhobar_dpbar = pstar / rhostar * drhodp;
913 double drhodp_Trbar = 1
914 / (R * Tr_bar * Tstar
917 double drhobar_dpbar_Trbar = pstar / rhostar * drhodp_Trbar;
918 double cp = HEOS.
cpmass();
919 double cv = HEOS.
cvmass();
920 double cpbar = cp / R;
921 double mubar = HEOS.
viscosity() / mustar;
922 double DELTAchibar_T = rhobar * (drhobar_dpbar - drhobar_dpbar_Trbar * Tr_bar / Tbar);
923 if (DELTAchibar_T < 0)
926 xi = xi_0 * pow(DELTAchibar_T / Lambda_0, nu / gamma);
927 double y = qd_bar * xi;
930 double kappa = cp / cv;
934 Z = 2 / (pi * y) * (((1 - 1 / kappa) * atan(y) + y / kappa) - (1 - exp(-1 / (1 / y + y * y / 3 / rhobar / rhobar))));
936 lambdabar_2 = GAMMA * rhobar * cpbar * Tbar / mubar * Z;
938 return (lambdabar_0 * lambdabar_1 + lambdabar_2) * lambdastar;
947 DeltaGstar = 2508.58,
950 DELTAlambda_max = 25,
955 double lambda_DG = B1 + B2 * T;
957 double rhobar = HEOS.
rhomolar() / 1000;
958 double lambda_L = C2 * (rhoL * rhoL) / (rhoL - rhobar) * sqrt(T) * exp(rhobar / (rhoL - rhobar) * DeltaGstar / (Ru * T));
960 double chi = rhobar - rhocbar;
963 double DELTAlambda_c = 4 * DELTAlambda_max / ((exp(chi) + exp(-chi)) * (exp(tau) + exp(-tau)));
965 return (pow((rhoL - rhobar) / rhoL, C1) * lambda_DG + pow(rhobar / rhoL, C1) * lambda_L + DELTAlambda_c) / 1e3;
977 double T = HEOS.
T(), Tc = 405.4, rhoc = 235, rho;
978 double LAMBDA = 1.2, nu = 0.63, gamma = 1.24, DELTA = 0.50, t, zeta_0_plus = 1.34e-10, a_zeta = 1, GAMMA_0_plus = 0.423e-8;
979 double pi = 3.141592654, a_chi, k_B = 1.3806504e-23, X_T, DELTA_lambda, dPdT, eta_B, DELTA_lambda_id, DELTA_lambda_i;
982 t = std::abs((T - Tc) / Tc);
983 a_chi = a_zeta / 0.7;
984 eta_B = (2.60 + 1.6 * t) * 1e-5;
985 dPdT = (2.18 - 0.12 / exp(17.8 * t)) * 1e5;
986 X_T = 0.61 * rhoc + 16.5 * log(t);
988 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
989 * GAMMA_0_plus * pow(t, -gamma) * (1 + a_chi * pow(t, DELTA));
990 DELTA_lambda_id = DELTA_lambda_i * exp(-36 * t * t);
991 if (rho < 0.6 * rhoc) {
992 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);
994 DELTA_lambda = DELTA_lambda_id * (X_T * X_T) / (X_T * X_T +
powInt(rho - 0.96 * rhoc, 2));
1004 double rhoc = 68.0, lambda_e, lambda_c, T = HEOS.
T(), rho = HEOS.
rhomass();
1005 double summer = 3.739232544 / T - 2.620316969e1 / T / T + 5.982252246e1 / T / T / T - 4.926397634e1 / T / T / T / T;
1006 double lambda_0 = 2.7870034e-3 * pow(T, 7.034007057e-1) * exp(summer);
1007 double c[] = {1.862970530e-4, -7.275964435e-7, -1.427549651e-4, 3.290833592e-5, -5.213335363e-8, 4.492659933e-8,
1008 -5.924416513e-9, 7.087321137e-6, -6.013335678e-6, 8.067145814e-7, 3.995125013e-7};
1010 lambda_e = (c[0] + c[1] * T + c[2] * pow(T, 1 / 3.0) + c[3] * pow(T, 2.0 / 3.0)) * rho
1011 + (c[4] + c[5] * pow(T, 1.0 / 3.0) + c[6] * pow(T, 2.0 / 3.0)) * rho * rho * rho
1012 + (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);
1017 if (3.5 < T && T < 12) {
1018 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;
1020 double DeltaT = std::abs(1 - T / Tc), DeltaRho = std::abs(1 - rho / rhoc_crit);
1025 double W = pow(DeltaT / 0.2, 2) + pow(DeltaRho / 0.25, 2);
1030 double x = pow(DeltaT / DeltaRho, 1 / beta);
1031 double h = E1 * (1 + x / x0) * pow(1 + E2 * pow(1 + x / x0, 2 / beta), (gamma - 1) / (2 * beta));
1042 * (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)
1043 + pow(beta, 2) * pow(E2 * pow((x + x0) / x0, 2 / beta) + 1, (1.0 / 2.0) * (2 * beta + gamma - 1) / beta))
1044 / (pow(beta, 2) * x0 * (E2 * pow((x + x0) / x0, 2 / beta) + 1));
1046 double RHS = pow(DeltaRho, delta - 1) * (delta * h - x / beta * dhdx);
1047 K_Tprime = 1 / (RHS * pow(rho / rhoc_crit, 2) * pc);
1048 K_Tbar = W * K_T + (1 - W) * K_Tprime;
1054 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));
1056 return lambda_0 + lambda_e + lambda_c;
1060 double delta = HEOS.
rhomolar() / 10139.0, tau = 190.55 / HEOS.
T();
1061 double lambda_dilute, lambda_residual, lambda_critical;
1066 0, -3.0328138281, 16.918880086, -37.189364917, 41.288861858, -24.615921140, 8.9488430959, -1.8739245042, 0.20966101390, -9.6570437074e-3};
1067 double OMEGA22_summer = 0;
1068 double t = HEOS.
T() / 174.0;
1069 for (
int i = 1; i <= 9; ++i) {
1070 OMEGA22_summer += C[i] * pow(t, (i - 1.0) / 3.0 - 1.0);
1072 double eta_dilute = 10.50 * sqrt(t) * OMEGA22_summer;
1073 double re[] = {0, 1, 1, 2, 2, 2, 3, 3, 4, 4, 1, 1};
1074 double se[] = {0, 0, 1, 0, 1, 1.5, 0, 2, 0, 1, 0, 1};
1075 double ge[] = {0, 0.41250137, -0.14390912, 0.10366993, 0.40287464, -0.24903524,
1076 -0.12953131, 0.06575776, 0.02566628, -0.03716526, -0.38798341, 0.03533815};
1077 double summer1 = 0, summer2 = 0;
1078 for (
int i = 1; i <= 9; ++i) {
1079 summer1 += ge[i] * pow(delta, re[i]) * pow(tau, se[i]);
1081 for (
int i = 10; i <= 11; ++i) {
1082 summer2 += ge[i] * pow(delta, re[i]) * pow(tau, se[i]);
1084 double eta_residual = 12.149 * summer1 / (1 + summer2);
1085 double eta = eta_residual + eta_dilute;
1088 double f_int = 1.458850 - 0.4377162 / t;
1089 lambda_dilute = 0.51828 * eta_dilute * (3.75 - f_int * (
POW2(HEOS.
tau()) * HEOS.
d2alpha0_dTau2() + 1.5));
1091 double rl[] = {0, 1, 3, 4, 4, 5, 5, 2};
1092 double sl[] = {0, 0, 0, 0, 1, 0, 1, 0};
1093 double jl[] = {0, 2.4149207, 0.55166331, -0.52837734, 0.073809553, 0.24465507, -0.047613626, 1.5554612};
1095 for (
int i = 1; i <= 6; ++i) {
1096 summer += jl[i] * pow(delta, rl[i]) * pow(tau, sl[i]);
1098 double delta_sigma_star = 1.0;
1102 lambda_residual = 6.29638 * (summer + jl[7] *
POW2(delta) / delta_sigma_star);
1104 double Tstar = 1 - 1 / tau;
1105 double rhostar = 1 - delta;
1106 double F_T = 2.646, F_rho = 2.678, F_A = -0.637;
1107 double F = exp(-F_T * sqrt(std::abs(Tstar)) - F_rho *
POW2(rhostar) - F_A * rhostar);
1109 if (std::abs(Tstar) < 0.03) {
1110 if (std::abs(rhostar) < 1e-16) {
1112 const double LAMBDA = 0.0801, gamma = 1.190;
1113 CHI_T_star = LAMBDA * pow(std::abs(Tstar), -gamma);
1114 }
else if (std::abs(rhostar) < 0.03) {
1116 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;
1117 double OMEGA = W * Tstar * pow(std::abs(rhostar), -1 / beta);
1119 if (Tstar < -pow(std::abs(rhostar), -1 / beta) / S) {
1120 theta = 1 + E * pow(1 + S * Tstar * pow(std::abs(rhostar), -1 / beta), 2 * beta);
1122 CHI_T_star = Q * pow(std::abs(rhostar), -a) * pow(theta, b) / (theta + OMEGA * (theta + R));
1133 * pow(CHI_T_star, 0.4681) * F;
1134 return (lambda_dilute + lambda_residual + lambda_critical) * 0.001;
1140 double resid = 9e30, resid_old = 9e30;
1153 r(0) = HEOS_Reference.
alphar() - alphar;
1163 Eigen::Vector2d v = J.colPivHouseholderQr().solve(-r);
1164 bool good_solution =
false;
1165 double T0_init = HEOS_Reference.
T(), rhomolar0_init = HEOS_Reference.
rhomolar();
1167 resid_old = sqrt(
POW2(r(0)) +
POW2(r(1)));
1168 for (
double frac = 1.0; frac > 0.001; frac /= 2) {
1171 double T_new = T0_init + frac * v(0);
1172 double rhomolar_new = rhomolar0_init + frac * v(1);
1176 if (resid > resid_old) {
1179 good_solution =
true;
1181 rhomolar0 = rhomolar_new;
1187 if (!good_solution) {
1192 throw ValueError(
format(
"conformal_state_solver took too many iterations; residual is %g; prior was %g", resid, resid_old));
1194 }
while (std::abs(resid) > 1e-9);
1207 for (std::size_t i = 0; i < ECS.
psi_a.size(); i++) {
1272 double psi_liq = 1 / (1 + exp(-100.0 * (x - 2)));
1275 const std::vector<double>&cL = data.
c_liq, cV = data.
c_vap;
1276 double f_liq = cL[0] + x * (cL[1] + x * (cL[2] + x * (cL[3])));
1277 double f_vap = cV[0] + x * (cV[1] + x * (cV[2] + x * (cV[3])));
1280 double etastar_ref = exp(psi_liq * f_liq + (1 - psi_liq) * f_vap);
1283 double etastar_fluid = 1 + data.
C * (etastar_ref - 1);
1285 return etastar_fluid * eta_dilute;
1293 R_kJkgK = R_u / M_kmol;
1300 for (std::size_t i = 0; i < ECS.
psi_a.size(); ++i) {
1306 for (std::size_t i = 0; i < ECS.
f_int_a.size(); ++i) {
1317 CoolPropDbl lambda_int = fint * eta_dilute * (cp0 - 2.5 * R) / 1e3;
1320 CoolPropDbl lambda_dilute = 15.0e-3 / 4.0 * R_kJkgK * eta_dilute;
1343 }
catch (std::exception& e) {
1344 throw ValueError(
format(
"Conformal state solver failed; error was %s", e.what()));
1364 CoolPropDbl lambda = lambda_int + lambda_dilute + lambda_resid * F_lambda + lambda_critical;